[TOC]
ARMA
Комбинация AR(p) и MA(q) представляет собой ARMA(p,q), авторегрессионную скользящую среднюю.
Формула выглядит следующим образом:
Формула говорит:
Значение текущего временного шага является константой плюс сумма авторегрессионного запаздывания и его множителя, плюс сумма запаздывания скользящего среднего и его множителя, плюс немного белого шума.
Он фиксирует влияние условий запаздывания и остатков и носит более общий характер. Определите порядок p, q и обновите коэффициенты уравнения в соответствии с непараметрическими оценками, такими как метод наименьших квадратов или оценка максимального правдоподобия.
Просмотрите процесс моделирования временных рядов:
- Проверка стационарности:
- Определить, является ли ряд стационарным
Если она не является стационарной, последовательность необходимо преобразовать (обычно с использованием разности);
- Определите, является ли стационарная последовательность белым шумом
Если стационарная последовательность представляет собой белый шум, условия моделирования не выполняются.
- Оценки модели:
- Определить значение p,q
Согласно историческим статьям, это можно определить по коэффициенту автокорреляции (ACF) и частичному коэффициенту автокорреляции (PACF), AR (p) представляет собой усечение p-порядка, MA (q) представляет собой усечение q-порядка;
- Информация
Если графики ACF и PACF не показывают четкого усечения, для оценки используется информационный критерий, как правило, с использованиемBIC,AIC
- Сочетание двух
- Модельный остаточный тест
- Распределяются ли остатки нормально со средним значением 0 и постоянной дисперсией (нормальностью)
- Тест на корреляцию остатков (корреляция)
ARIMA
Разница между авторегрессионным интегрированным скользящим средним (ARIMA) и ARMA заключается в том, что существует дополнительный параметр d, который преобразует нестационарную последовательность в стационарную последовательность, что означает, что разница d-порядка преобразуется в стационарную последовательность. Модель ARIMA — это просто модель ARMA в разном временном ряду.
Символы для моделей ARIMAARIMA(p, d, q) Выражать.
Например, в модели ARIMA(1,1,0) (1,1,0) означает, что имеется авторегрессионное запаздывание, данные являются первыми разностями и член скользящего среднего отсутствует.
- p
Авторегрессивная часть модели включает в модель влияние прошлых значений, то есть исторические значения оказывают влияние на будущее;
- d — интегральная часть модели.
Количество разностей, необходимых для того, чтобы сделать временной ряд стационарным. Например, если разница температур за последние три дня очень мала, завтрашняя температура может быть аналогична температуре предыдущих дней;
- q
Скользящая средняя часть модели, ошибка модели может быть линейной комбинацией значений ошибок, наблюдаемых в прошлые моменты времени.
SARIMA
SARIMA (сезонная авторегрессионная интегрированная модель скользящего среднего), сезонная авторегрессионная модель скользящего среднего с экзогенной регрессионной моделью, называемая季节性ARIMA. То есть на базе ARIMA добавляется сезонная часть. Сезонность относится к повторяющимся закономерностям в данных с фиксированной частотой: закономерностям, которые повторяются каждый день, каждые две недели, каждые четыре месяца и так далее.
Модель SARIMA может быть выражена как SARIMA（p，d，q）x（P，D，Q）s, формула удовлетворяет принципу умножения, первая половина представляет несезонную часть, последняя представляет сезонную часть, а s представляет сезонную частоту.
季节性成分可能捕捉长期模式，而非季节性成分调整了对短期变化的预测.
САРИМА бой
Во-первых, давайте последовательно пройдем процесс моделирования анализа временных рядов.
1. Тест на стационарность последовательности
используется здесь单位根检验.
Единичный корневой тест: Проверка единичного корня временного ряда является проверкой стационарности временного ряда.Если в нестационарном временном ряду есть единичный корень, единичный корень обычно можно исключить методом разности для получения стационарного ряд.
Единичный корень t-критерий:
- Нулевая гипотеза: существует единичный корень
- p
def test_stationarity(timeseries,
maxlag=None, regression=None, autolag=None,
window=None, plot=False, verbose=False):
'''
单位根检验
'''
# set defaults (from function page)
if regression is None:
regression = 'c'
if verbose:
print('Running Augmented Dickey-Fuller test with paramters:')
print('maxlag: {}'.format(maxlag))
print('regression: {}'.format(regression))
print('autolag: {}'.format(autolag))
if plot:
if window is None:
window = 4
#Determing rolling statistics
rolmean = timeseries.rolling(window=window, center=False).mean()
rolstd = timeseries.rolling(window=window, center=False).std()
#Plot rolling statistics:
orig = plt.plot(timeseries, color='blue', label='Original')
mean = plt.plot(rolmean, color='red', label='Rolling Mean ({})'.format(window))
std = plt.plot(rolstd, color='black', label='Rolling Std ({})'.format(window))
plt.legend(loc='best')
plt.title('Rolling Mean & Standard Deviation')
plt.show(block=False)
#Perform Augmented Dickey-Fuller test:
dftest = smt.adfuller(timeseries, maxlag=maxlag, regression=regression, autolag=autolag)
dfoutput = pd.Series(dftest[0:4], index=['Test Statistic',
'p-value',
'#Lags Used',
'Number of Observations Used',
])
for key,value in dftest[4].items():
dfoutput['Critical Value (%s)'%key] = value
if verbose:
print('Results of Augmented Dickey-Fuller Test:')
print(dfoutput)
return dfoutput
2, acf, pacf диаграмма
Нарисуйте исходную диаграмму последовательности, диаграммы ACF и PACF и грубо оцените тенденцию исторических данных и порядок p, q последовательности.
def tsplot(y, lags=None, title='', figsize=(14, 8)):
fig = plt.figure(figsize=figsize)
layout = (2, 2)
ts_ax = plt.subplot2grid(layout, (0, 0))
hist_ax = plt.subplot2grid(layout, (0, 1))
acf_ax = plt.subplot2grid(layout, (1, 0))
pacf_ax = plt.subplot2grid(layout, (1, 1))
y.plot(ax=ts_ax)
ts_ax.set_title(title)
y.plot(ax=hist_ax, kind='hist', bins=25)
hist_ax.set_title('Histogram')
smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)
smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)
[ax.set_xlim(0) for ax in [acf_ax, pacf_ax]]
sns.despine()
fig.tight_layout()
return ts_ax, acf_ax, pacf_ax
3. Остаточная статистика модели
Тест на нормальность стандартизованных остатков (критерий нормальности Жака-Бера):
- Нулевая гипотеза: это нормально
- p>альфа, принять нулевую гипотезу, остатки в норме
Тест на серийную корреляцию остатков (тест Льюнга-Бокса):
- Нулевая гипотеза: нет последовательной корреляции
- p>alpha , принять нулевую гипотезу, остаточный ряд не имеет корреляции
Тест остаточной серийной корреляции (тест Дарбина-Ватсона):
- Чем ближе статистическое значение к 2, тем лучше, и обычно между 1 и 3 проблем нет;
- Меньше 1 указывает, что остатки являются автокоррелированными.
def model_resid_stats(model_results,
het_method='breakvar',
norm_method='jarquebera',
sercor_method='ljungbox',
verbose=True,
):
(het_stat, het_p) = model_results.test_heteroskedasticity(het_method)[0]
# Jarque-Bera 正态性检验
norm_stat, norm_p, skew, kurtosis = model_results.test_normality(norm_method)[0]
# Ljung-Box检验 相关性检验
sercor_stat, sercor_p = model_results.test_serial_correlation(method=sercor_method)[0]
sercor_stat = sercor_stat[-1] # last number for the largest lag
sercor_p = sercor_p[-1] # last number for the largest lag
# Durbin-Watson检验 相关性检验
dw_stat = sm.stats.stattools.durbin_watson(model_results.filter_results.standardized_forecasts_error[0, model_results.loglikelihood_burn:])
# check whether roots are outside the unit circle (we want them to be);
# will be True when AR is not used (i.e., AR order = 0)
arroots_outside_unit_circle = np.all(np.abs(model_results.arroots) > 1)
# will be True when MA is not used (i.e., MA order = 0)
maroots_outside_unit_circle = np.all(np.abs(model_results.maroots) > 1)
if verbose:
print('Test heteroskedasticity of residuals ({}): stat={:.3f}, p={:.3f}'.format(het_method, het_stat, het_p));
print('\nTest normality of residuals ({}): stat={:.3f}, p={:.3f}'.format(norm_method, norm_stat, norm_p));
print('\nTest serial correlation of residuals ({}): stat={:.3f}, p={:.3f}'.format(sercor_method, sercor_stat, sercor_p));
print('\nDurbin-Watson test on residuals: d={:.2f}\n\t(NB: 2 means no serial correlation, 0=pos, 4=neg)'.format(dw_stat))
print('\nTest for all AR roots outside unit circle (>1): {}'.format(arroots_outside_unit_circle))
print('\nTest for all MA roots outside unit circle (>1): {}'.format(maroots_outside_unit_circle))
stat = {'het_method': het_method,
'het_stat': het_stat,
'het_p': het_p,
'norm_method': norm_method,
'norm_stat': norm_stat,
'norm_p': norm_p,
'skew': skew,
'kurtosis': kurtosis,
'sercor_method': sercor_method,
'sercor_stat': sercor_stat,
'sercor_p': sercor_p,
'dw_stat': dw_stat,
'arroots_outside_unit_circle': arroots_outside_unit_circle,
'maroots_outside_unit_circle': maroots_outside_unit_circle,
}
return stat
4. Поиск по сетке параметров модели
В модели SARIMA есть 6 параметров.Если вы выбираете его вручную, вы должны настроить облысение, поэтому вам нужно настроить поиск по сетке.На самом деле, вы также можете использовать байесовскую оценку для настройки параметров.Введена только сетка здесь.
Здесь мы в основном сосредоточимся на вызове функции SARIMAX и значении ее параметров.
- порядок: (p,d,q) соответствует ARIMA
- сезонный_порядок: (P,D,Q,s)
def model_gridsearch(ts,
p_min,
d_min,
q_min,
p_max,
d_max,
q_max,
sP_min,
sD_min,
sQ_min,
sP_max,
sD_max,
sQ_max,
trends,
s=None,
enforce_stationarity=True,
enforce_invertibility=True,
simple_differencing=False,
plot_diagnostics=False,
verbose=False,
filter_warnings=True,
):
'''Run grid search of SARIMAX models and save results.
'''
cols = ['p', 'd', 'q', 'sP', 'sD', 'sQ', 's', 'trend',
'enforce_stationarity', 'enforce_invertibility', 'simple_differencing',
'aic', 'bic',
'het_p', 'norm_p', 'sercor_p', 'dw_stat',
'arroots_gt_1', 'maroots_gt_1',
'datetime_run']
# Initialize a DataFrame to store the results
df_results = pd.DataFrame(columns=cols)
mod_num=0
for trend,p,d,q,sP,sD,sQ in itertools.product(trends,
range(p_min,p_max+1),
range(d_min,d_max+1),
range(q_min,q_max+1),
range(sP_min,sP_max+1),
range(sD_min,sD_max+1),
range(sQ_min,sQ_max+1),
):
# initialize to store results for this parameter set
this_model = pd.DataFrame(index=[mod_num], columns=cols)
if p==0 and d==0 and q==0:
continue
try:
model = sm.tsa.SARIMAX(ts,
trend=trend,
order=(p, d, q),
seasonal_order=(sP, sD, sQ, s),
enforce_stationarity=enforce_stationarity,
enforce_invertibility=enforce_invertibility,
simple_differencing=simple_differencing,
)
if filter_warnings is True:
with warnings.catch_warnings():
warnings.filterwarnings("ignore")
model_results = model.fit(disp=0)
else:
model_results = model.fit()
if verbose:
print(model_results.summary())
if plot_diagnostics:
model_results.plot_diagnostics();
stat = model_resid_stats(model_results,
verbose=verbose)
this_model.loc[mod_num, 'p'] = p
this_model.loc[mod_num, 'd'] = d
this_model.loc[mod_num, 'q'] = q
this_model.loc[mod_num, 'sP'] = sP
this_model.loc[mod_num, 'sD'] = sD
this_model.loc[mod_num, 'sQ'] = sQ
this_model.loc[mod_num, 's'] = s
this_model.loc[mod_num, 'trend'] = trend
this_model.loc[mod_num, 'enforce_stationarity'] = enforce_stationarity
this_model.loc[mod_num, 'enforce_invertibility'] = enforce_invertibility
this_model.loc[mod_num, 'simple_differencing'] = simple_differencing
this_model.loc[mod_num, 'aic'] = model_results.aic
this_model.loc[mod_num, 'bic'] = model_results.bic
# this_model.loc[mod_num, 'het_method'] = stat['het_method']
# this_model.loc[mod_num, 'het_stat'] = stat['het_stat']
this_model.loc[mod_num, 'het_p'] = stat['het_p']
# this_model.loc[mod_num, 'norm_method'] = stat['norm_method']
# this_model.loc[mod_num, 'norm_stat'] = stat['norm_stat']
this_model.loc[mod_num, 'norm_p'] = stat['norm_p']
# this_model.loc[mod_num, 'skew'] = stat['skew']
# this_model.loc[mod_num, 'kurtosis'] = stat['kurtosis']
# this_model.loc[mod_num, 'sercor_method'] = stat['sercor_method']
# this_model.loc[mod_num, 'sercor_stat'] = stat['sercor_stat']
this_model.loc[mod_num, 'sercor_p'] = stat['sercor_p']
this_model.loc[mod_num, 'dw_stat'] = stat['dw_stat']
this_model.loc[mod_num, 'arroots_gt_1'] = stat['arroots_outside_unit_circle']
this_model.loc[mod_num, 'maroots_gt_1'] = stat['maroots_outside_unit_circle']
this_model.loc[mod_num, 'datetime_run'] = pd.to_datetime('today').strftime('%Y-%m-%d %H:%M:%S')
df_results = df_results.append(this_model)
mod_num+=1
except:
continue
return df_results
5. Постройте модель
5.1 Импорт данных
Набор тестов и набор проверки здесь разделены, что отличается от моделирования машинного обучения в форме случайной выборки, потому что данные временных рядов упорядочены.
liquor = pd.read_csv('liquor.csv', header=0, index_col=0, parse_dates=[0])
n_sample = liquor.shape[0]
n_train=int(0.95*n_sample)+1
n_forecast=n_sample-n_train
# 拆分测试集序列和验证集序列
liquor_train = liquor.iloc[:n_train]['Value']
liquor_test = liquor.iloc[n_train:]['Value']
print(liquor_train.shape)
print(liquor_test.shape)
print("Training Series:", "\n", liquor_train.tail(), "\n")
print("Testing Series:", "\n", liquor_test.head())
5.2 Визуализация исходной последовательности, acf и pacf
tsplot(liquor_train, title='Liquor Sales (in millions of dollars)', lags=40);
Из исходной диаграммы последовательности видно, что последовательность не является стационарной последовательностью, а из диаграмм acf и pacf нет очевидного усечения, и невозможно судить о p и q.
5.3 Преобразование нестационарных рядов в стационарные ряды
# 检验平稳性
test_stationarity(liquor_train)
Тест на единичный корень, p>0,05, не может отвергнуть нулевую гипотезу (с единичным корнем), ряд нестационарен.
# 差分
test_stationarity(liquor_train.diff().dropna())
Разность первого порядка, p
5.4 Поиск сетки параметров модели
p_min = 0
d_min = 0
q_min = 0
p_max = 2
d_max = 1
q_max = 2
sP_min = 0
sD_min = 0
sQ_min = 0
sP_max = 1
sD_max = 1
sQ_max = 1
# 以一年为一个周期
s=12
# trends=['n', 'c']
trends=['n']
enforce_stationarity=True
enforce_invertibility=True
simple_differencing=False
plot_diagnostics=False
verbose=False
df_results = model_gridsearch(liquor['Value'],
p_min,
d_min,
q_min,
p_max,
d_max,
q_max,
sP_min,
sD_min,
sQ_min,
sP_max,
sD_max,
sQ_max,
trends,
s=s,
enforce_stationarity=enforce_stationarity,
enforce_invertibility=enforce_invertibility,
simple_differencing=simple_differencing,
plot_diagnostics=plot_diagnostics,
verbose=verbose,
)
5.5 Выбор модели и конструкция
df_results.sort_values(by='bic').head(10)
Здесь BIC выбран в качестве показателя оценки модели,选择最小的BIC对应的参数Выполняется моделирование, т.е. (p,d,q)=(2,1,0), (P,D,Q)s = (0,1,1)12.
Подставим указанные выше оптимальные параметры в модель:
mod = sm.tsa.statespace.SARIMAX(liquor_train, order=(2,1,0), seasonal_order=(0,1,1,12))
sarima_fit2 = mod.fit()
print(sarima_fit2.summary())
Давайте посмотрим на статистику полученной модели обучающего набора:
- coef : Коэффициент регрессии
- p>|z|: является ли коэффициент значимым
- JB: тест нормальности для остатков
- LB: корреляционный тест для остаточных последовательностей
Визуальный осмотр остатков модели:
- Случайность: белый это шум или нет
- Нормальность: это нормальное распределение
- Корреляция: низкая ли корреляция между остатками?
Только при выполнении вышеуказанных условий модель успешно устанавливается.
sarima_fit2.plot_diagnostics(figsize=(16, 12));
5.6 Прогноз
fig, ax1 = plt.subplots(nrows=1, ncols=1, figsize=(12, 8))
ax1.plot(liquor_train, label='In-sample data', linestyle='-')
# subtract 1 only to connect it to previous point in the graph
ax1.plot(liquor_test, label='Held-out data', linestyle='--')
# yes DatetimeIndex
pred_begin = liquor_train.index[sarima_fit2.loglikelihood_burn]
pred_end = liquor_test.index[-1]
pred = sarima_fit2.get_prediction(start=pred_begin.strftime('%Y-%m-%d'),
end=pred_end.strftime('%Y-%m-%d'))
pred_mean = pred.predicted_mean
pred_ci = pred.conf_int(alpha=0.05)
ax1.plot(pred_mean, 'r', alpha=.6, label='Predicted values')
ax1.fill_between(pred_ci.index,
pred_ci.iloc[:, 0],
pred_ci.iloc[:, 1], color='k', alpha=.2)
ax1.set_xlabel("Year")
ax1.set_ylabel("Liquor Sales")
ax1.legend(loc='best');
fig.tight_layout();
Увеличьте масштаб, эффект подгонки очень хороший~ (стабильно видно невооруженным глазом)
Ссылка на ссылку:
- статистика models.source forge.net/but eve/root…
- wiki.MiaoParisBay.com/wiki/%E5%8D…
- Woohoo.stats models.org/Dev/ аааа он…
- Woohoo.stats models.org/Dev/ аааа он…
- Woohoo.stats models.org/Dev/ аааа он…
- Woohoo.stats models.org/Dev/ аааа он…
- cloud.Tencent.com/developer/ ах…
- blog.CSDN.net/Снятие Муму…
Приглашаем обратить внимание на личный публичный номер:Distinct数说
