В прошлой части мы с вами остановились на том, что обнаружили у временного ряда с температурой две сезонности и, несмотря на это, решили двигаться дальше в выполнении сезонной модели САРПСС по методологии АРПСС. В этой части второй главы мы с вами продолжим применение методологии для поиска оптимальных параметров модели, которая будет адекватно описывать целевой временной ряд с температурой.
Напоминаю, что в нашем ряду проявляются сильные сезонные зависимости и имеется тенденция к увеличению температуры со временем. Давайте выполним дополнительную проверку и разложим оригинальный временной ряд на три компоненты - годовую сезонность, тенденцию и случайную компоненту. Для этого применим метод statsmodels seasonal_decompose(), установив период равный 8760-ти шагам.
По выполненным в предыдущих частях графикам с температурой можно предположить, что амплитуда колебаний годовой сезонности остаётся относительно постоянной с увеличением тенденции. Данное предположение оправдывает построение аддитивной модели. Хотя Дж.Бокс и Г.Дженкинс из своего опыта советовали «в качестве начального приближения» строить мультипликативную модель, метод seasonal_decompose() не даст нам это сделать по причине присутствия множества нулевых значений во временном ряду. Конечно, было бы правильным потратить n-ое время на преобразования данных и построить мультипликативную модель, тем более что концы у каждой годовой дуги разнятся. Однако предлагаю выполнить аддитивную модель: будем работать с тем, что есть.
from statsmodels.tsa.seasonal import seasonal_decompose
def show_seasonal_decompose_plot(series, period):
"""Декомпозиция временного ряда на компоненты"""
decomposition = seasonal_decompose(series, period=period, model='additive')
fig, axs = pyplot.subplots(4, 1)
axs[0].plot(series, color='blue', alpha=0.8)
axs[0].set_title('Оригинальный временной ряд')
axs[1].plot(decomposition.seasonal, color='coral')
axs[1].set_title('Сезонная компонента')
axs[2].plot(decomposition.trend, color='green')
axs[2].set_title('Тенденция')
axs[3].plot(decomposition.resid, color='blue', alpha=0.3)
axs[3].set_title('Случайная компонента')
pyplot.suptitle(f"Декомпозиция временного ряда (период={period})")
pyplot.show()
show_seasonal_decompose_plot(series=temperature, period=365*24)
Вдобавок построим отдельный график АКФ, чтобы подтвердить годовую сезонность в 8760 шагов:
def show_autocorrelation_plot(series, name_column='Температура'):
"""График автокорреляции временного ряда"""
plot_acf(series, lags=arange(len(series)), marker='')
pyplot.xlabel('Шаг запаздывания', fontsize=12)
pyplot.ylabel('Автокорреляция', fontsize=12)
pyplot.title(f"График автокорреляции временного ряда "{name_column}"")
pyplot.grid(True)
pyplot.show()
show_autocorrelation_plot(series=temperature)
Итак, сразу бросается в глаза синусообразная форма автокорреляционной функции, которая затухает с увеличением количества шагов запаздывания. Годовой период равняется 8760 шагам и, как хорошо видно по оригинальному временному ряду и по сезонной компоненте разложенного ряда, имеет повторяющуюся форму в виде дуги, которую можно попробовать смоделировать с помощью математической функции. Дневной период не имеет чётко выраженной структуры (далее на графике с температурой, разграниченной на части, мы с вами убедимся в этом): несмотря на то, что в дневных данных есть тенденция к чередованию значений «вверх-вниз», будущие значения не могут быть точно описаны математической функцией. В этой связи можно предположить, что целевой временной ряд, состоящий из случайных элементов, может быть частично детерминированным по годовой сезонности.
Также в литературе по прогнозированию временных рядов часто отмечается, что синусообразная и «бесконечная» (именно такой термин используют Дж.Бокс и Г.Дженкинс) форма автокорреляционной функции может указывать на необходимость дополнительного взятия разности (параметр d). Однако мы должны помнить, что «конечную разность берут столько раз, сколько необходимо, чтобы обеспечить стационарность» [194 с.], поэтому полностью доверяться поведению АКФ с сильной дневной зависимостью в данных с более 70 тыс. строк не стоит. В дальнейшем мы выполним анализ остатков прогнозной модели, включая проверку в них автокорреляции, чтобы выяснить, имеет ли смысл увеличить порядок разности с 1 до 2 или нет.
Для более точного определения параметров p и q временного ряда с разностью 1-го порядка повторно выполним команду из первой части, но уменьшим количество шагов временного запаздывания до 25:
plot_target_analysis(temperature, lags=25)
Хотя ЧАКФ и обрывается после первого шага, далее она представлена в форме плавно затухающей синусоиды и бесконечна. Форма затухающей синусоиды у обеих функций - АКФ и ЧАКФ - говорит о том, что мы имеем дело со смешанным процессом авторегрессии - скользящего среднего (АРСС).
Отдельно отмечу, что Дж.Бокс и Г.Дженкинс в своей книге разбирают случаи поведения АКФ и ЧАКФ, поэтому для общего понимания рекомендую ознакомиться с таблицей 3.2 «Сводка характеристик процессов авторегрессии, скользящего среднего и смешанных (АРСС) процессов» (97 с. книги 1-го выпуска).
По верхнему графику АКФ видно, что большая часть шагов запаздывания значительно удалена от доверительного интервала, который в данном случае вовсе не заметен, поэтому мы не можем сразу определить какое значение q взять за основу. Вообще говоря, в книгах по прогнозированию временных рядов отмечается - в частности в полюбившейся нам книге Дж.Бокса и Г.Дженкинса - что «адекватное описание наблюдаемых временных рядов достигается при помощи моделей авторегресии, скользящего среднего или комбинированной модели, в которых p и q не больше, а часто и меньше 2» [26 с.]. При этом коллеги эконометристы оставляют важное дополнение: «если нет сезонной компоненты» [из книги «Эконометрика. Начальный курс: Учеб.» - 5-е изд., испр. / Магнус Я.Р., Катышев П.К., Пересецкий А.А. - М.: Дело, 2001. - 400]. Учитывать их опыт необходимо, но ограничивать себя этими рамками, имея дело с «большими» временными рядами с сильными сезонными особенностями, было бы ошибкой.
Нужно понимать, что наши данные далеко не игрушечные: они имеют сильное влияние сезонной компоненты (дневная и годовая сезонности) и их объём значителен (к сведению: набор временных рядов конкурса М3 из приведённого ранее исследования 2018 года представлен временными рядами с тремя тысячами наблюдений [ссылка]). Поэтому для начала предлагаю принять q равным 24. Дополнительно возьмём значения 24 и 48.
Обращаю внимание, что нулевые шаги запаздывания как в АКФ, так и в ЧАКФ неинформативны в плане определения зависимостей между наблюдениями в разные моменты времени, поскольку представляют собой корреляцию временного ряда сами с собой в тот же самый момент времени. И они всегда равны 1. Поэтому их не берут в расчёт.
На графике с ЧАКФ видно, что первое запаздывание представлено в виде обособленной пики, что является хорошим знаком в плане наличия статистически значимой автокорреляции на этом шаге, поэтому начальное значение для p определим равным 1. Однако и здесь наблюдается удалённость большей части запаздывающих шагов от доверительного интервала; к тому же на графиках с ЧАКФ (lags=120 и 25) также выявилась сезонность в шагах запаздывания. При этом видно, что наиболее предпочтителен для дальнейшего использования и анализа 23-й шаг запаздывания, поскольку он наиболее удалён от нулевой линии, чем соседние шаги. Но и здесь мы не будем себя сдерживать: для учёта краткосрочных влияний процесса авторегрессии помимо значений 1 и 23 возьмём также 24. Вдобавок, ради интереса прихватим значение 46, чтобы посмотреть, способно ли дальнейшее увеличение количества шагов запаздывания повлиять на результат.
Важно отметить, что принимая такие значения запаздываний, как 23 и 24, в качестве основы мы не только пытаемся с их помощью учесть краткосрочные (дневные) зависимости процессов АРСС, но и обыграть ограничение модели SARIMA/SARIMAX в одну сезонность, которую отнесём к долгосрочной (годовой) зависимости с использованием параметра m=8760. При этом определение моделей с «избыточными» параметрами (значения 46, 48) оправдано не только из-за наличия в данных двух сезонностей. Выполнение диагностической проверки с использованием метода введения избыточного числа параметров позволит выявить («угадать») неадекватные свойства модели: «Идентифицировав модель и считая её правильной [в нашем случае мы предполагаем, что это ARIMA[23, 1, 24]. - Примечание], мы затем подгоняем более сложную модель. Это ставит под угрозу идентифицированную модель, потому что более сложная модель содержит дополнительные параметры, с помощью которых можно ликвидировать возможные отклонения. Нужно тщательно продумать вопрос о том, как следует усложнить модель.» [из книги «Анализ временных рядов: прогноз и управление /Дж.Бокс и Г.Дженкинс. Выпуск 1. - Издательство «Мир», Москва 1974. - 314 с.»]
Итак, с учётом результатов статистических тестов для определения стационарности и анализа вышерассмотренных графиков АКФ и ЧАКФ примем следующие порядки (p, d, q) для модели АРПСС:
arima_orders = [
(1, 1, 24),
(23, 1, 1),
(23, 1, 24),
(24, 1, 24),
(46, 1, 48)
]Заметим, что в начале мы начинаем с порядков, смешанный процесс которых ограничен одним «сильным» процессом — 24 для процесса скользящего среднего и 23 для процесса авторегрессии. Это поможет лучше понять «природу прогнозирования» временного ряда с температурой.
В дальнейших расчётах мы будем использовать модель SARIMAX библиотеки statsmodels, но поскольку на начальном этапе мы не учитываем сезонную компоненту и экзогенные регрессоры, то для несезонных порядков выполним модель ARIMA той же библиотеки.
Сезонный порядок (P, D, Q, m) значительно увеличит время расчётов и потребление памяти, поэтому в целях экономии времени рассмотрение всех временных рядов я выношу за рамки данного исследования. В ходе дальнейших шагов мы рассмотрим только следующие из них, которые были выбраны выборочно («на глаз») с учётом разнообразия состояния (стационарное/нестационарное), без корреляционного анализа или анализа связи по тесту Грейнджера:
colls = [
'p (mbar)',
'T (degC)',
'VPmax (mbar)',
'sh (g/kg)',
'rho (g/m**3)',
'wv (m/s)'
]При этом в расчётах буду задействованы не все годы, а только первые четыре. Нижеследующий график наглядно продемонстрирует разделение данных. Как можно догадаться, принимая эти ограничения мы хорошенько подготавливаемся к борьбе с вычислительной затратностью модели SARIMAX применительно к годовой сезонности в 8760 шагов.
def show_plot_split_temperature():
"""График температуры, разграниченной по частям"""
fig, axs = pyplot.subplots(2, 1)
axs[0].plot(temperature['2009-01-01':'2011-12-31'], label='Обучающая часть')
axs[0].plot(temperature['2012-01-01':'2012-01-07'], label='Валидационная часть')
axs[0].plot(temperature['2012-01-08':'2012-01-14'], label='Тестовая часть')
axs[0].set_ylabel('Температура (град. C)')
axs[0].legend(loc='best')
axs[1].plot(temperature['2009-01-01':'2011-12-31'][-168:], label='Обучающая часть')
axs[1].plot(temperature['2012-01-01':'2012-01-07'], label='Валидационная часть')
axs[1].plot(temperature['2012-01-08':'2012-01-14'], label='Тестовая часть')
axs[1].set_xlabel('Дата наблюдения')
axs[1].set_ylabel('Температура (град. C)')
axs[1].legend(loc='best')
pyplot.suptitle('График температуры, разграниченной по частям')
pyplot.show()
show_plot_split_temperature()
Период [2009-01-01 : 2011-12-31] включительно будет использоваться для обучения моделей, следующий за ним период [2012-01-01 : 2012-01-07] (то есть 7 дней) - для их валидации. Тестовая часть, охватывающая период [2012-01-08 : 2012-01-14] (следующие 7 дней после валидационной части), будет использоваться для сравнения с прогнозом окончательно определившейся модели по результатам отбора.
Я дополнительно выполнил проверку на стационарность отобранных данных с учётом их сокращения согласно периоду [2009-01-01 : 2011-12-31]: отличий от ранее выполненных тестов - расширенного теста Дики-Фуллера и теста КПСС - не выявлено.
for column in colls:
print('n'f"Проверка стационарности временного ряда {column}")
series_p_value_verification(df[column].loc[:'2011-12-31'])Целью данного исследования в части выполнения моделей является выполнение прогнозирования почасовых значений температуры на следующие 7 дней, что с учётом почасового формата данных составит горизонт прогнозирования, равный 168 наблюдениям.
Оценивать модели мы будем с помощью двух основных подходов к оценке модели в прогнозировании временных рядов - предсказание в пределах выборки и прогнозирование вне выборки (в англоязычной терминологии in-sample prediction и out-of-sample forecasting соответственно). Чаще всего к данным, которые находятся в пределах выборки, относят как обучающие, так и валидационные данные, которые затем объединяются друг с другом для окончательного обучения модели, и уже после этого выполняется прогноз на основе последнего известного значения, который сравнивается с тестовыми данными (с вневыборочными данными).
Данный подход к оцениванию модели хорошо объясняет разницу между понятиями «предсказать» и «спрогнозировать» применительно к временным рядам: предсказание — это оценка известных данных (обучающих и валидационных), а прогноз — это оценка новых данных (тестовых, как в нашем случае), которые не использовались при обучении модели.
Это важно запомнить, поскольку по этой же логике выполнено несколько методов моделей ARIMA/SARIMAX библиотеки statsmodels. Методы predict() и get_prediction() применяются для выполнения внутривыборочных предсказаний, forecast() и get_forecast() — для выполнения вневыборочных прогнозов на основе последнего известного значения. При этом методы predict() и forecast() возвращают только предсказанные/спрогнозированные значения, а методы get_prediction() и get_forecast() позволяют получить «дополнительные результаты» [ссылка] - как предсказанные/спрогнозированные значения (mean), так и граничные значения интервала предсказания/прогнозирования для этих значений, но никак не доверительные интервалы, как, по всей видимости, ошибочно написано в документации statsmodels к этим методам.
Доверительный интервал обычно используется для оценки надёжности параметров модели (далее мы убедимся в этом по итоговой таблице с результатами обученной модели), в то время как интервал прогнозирования охватывает диапазон возможных значений, которые могут быть получены при применении модели по отношению к новым данным. Поскольку метод get_forecast() предоставляет прогнозы для будущих наблюдений, то логично, что интервал, который он возвращает, должен отражать неопределённость этих прогнозов, а не оценок параметров модели.
Возвращаясь к нашим данным, я предлагаю не выполнять предсказания обученной модели на обучающих данных, поскольку это будет неинформативно в плане обобщения: модели будут показывать очень высокий результат на обучающих данных и жуткий на валидационных данных, что связано с такой проблемой, как переобучение модели. В этой связи о качестве моделей мы будем судить как по информационному критерию Акаике AIC [ссылка на оригинальную статью] из результатов обученной модели, так и по ошибкам на валидационных данных с помощью метрики «Средняя абсолютная процентная ошибка» (в англоязычной терминологии - mean absolute percentage error, mape). Чем они меньше, тем лучше модель описывает данные. Тестовые данные будут использоваться только для сравнения с прогнозом, выполненного определившеюся моделью по результатам отбора.
Итак, с помощью нижеприведённого кода мы для каждого порядка из списка arima_orders последовательно выполним следующие действия: обучим модель, выполним предсказания по меткам валидационных данных на 168 шагов и рассчитаем ошибку mape, а результаты (AIC, BIC, preds_mean и mape) сохраним в словарь results, который как результат функции build_arima_model() вложим в общий для всех порядков (=несезонных моделей) словарь sum_results. По ключам и элементам этого общего словаря мы будем оценивать эффективность выполнения предсказывания рассматриваемых моделей.
train_data = temperature['2009-01-01':'2011-12-31'].copy()
val_data = temperature['2012-01-01':'2012-01-07'].copy()
arima_orders = [
(1, 1, 24),
(23, 1, 1),
(23, 1, 24),
(24, 1, 24),
(46, 1, 48)
]
sum_results = {}
def build_arima_model(order):
"""Обучение модели ARIMA и сохранение её результатов"""
results = {}
arima_model = ARIMA(endog=train_data,
order=order,
trend=None,
enforce_stationarity=False,
enforce_invertibility=False)
arima_model = arima_model.fit()
preds = arima_model.get_prediction(start=val_data.index.min(),
end=val_data.index.max())
mape = mean_absolute_percentage_error(y_true=val_data,
y_pred=preds.predicted_mean)
results['AIC'] = arima_model.aic
results['BIC'] = arima_model.bic
results['preds_mean'] = preds.predicted_mean
results['MAPE'] = mape
return results
for order in arima_orders:
key = f"ARIMA{order}"
print(f"Выполняется обучение модели {key}")
sum_results[key] = build_arima_model(order)Завершив выполнение расчётов, объединим результаты и выведем их на общем графике с помощью следующей функции:
def show_plot_prediction(model_dict):
"""График сравнения истинных значений с предсказанными"""
pyplot.figure()
pyplot.plot(val_data, label='Валидационные данные')
for key, model_data in model_dict.items():
predictions = model_data['preds_mean']
label_order = f"{key} "
label_mape = f"[mape: {model_data['MAPE']*100:.2f}%] "
label_aic = f"[AIC: {int(model_data['AIC'])}]"
pyplot.plot(predictions, label=(label_order + label_mape + label_aic))
pyplot.xlabel('Дата наблюдения')
pyplot.ylabel('Температура (град C)')
pyplot.title('Предсказания на валидационных данных')
pyplot.grid(True)
pyplot.legend(loc='best')
pyplot.show()
show_plot_prediction(model_dict=sum_results)
По результатам видно, что модель ARIMA(1, 1, 24) показала наихудший результат из всех по критерию AIC, а предсказания этой модели представлены в виде практически прямой линии с минимальными колебаниями, которые наблюдаются только в начале предсказания.
ARIMA(23, 1, 1) показала значительно лучший результат по AIC, но по mape она оказалась немного хуже. Предсказания этой модели уже не выглядят как постоянная прямая, а демонстрируют незначительные колебания.
Совмещение процесса авторегрессии со значением 23 и процесса скользящего среднего со значением 24 в модели ARIMA(23, 1, 24) значительно повысило эффективность как по AIC, так и по mape. Предсказания модели стали более вариативными, учитывая подъёмы и спады в данных.
ARIMA(24, 1, 24) показало немного лучший результат по AIC, чем ARIMA(23, 1, 24), но хуже по mape. Предсказания обеих моделей практически идентичны.
Дальнейшее увеличение количества временных запаздываний, представленное в модели ARIMA(46, 1, 48) лишь незначительно улучшило AIC, но при этом ошибка mape ухудшилась по сравнению с моделями ARIMA(23, 1, 24) и ARIMA(24, 1, 24).
Таким образом, заключаем, что модель ARIMA(23, 1, 24) оказалась наиболее подходящей для дальнейшего использования в силу баланса между минимизацией критериев оценки модели AIC и mape.
Важно понимать, что использование модели ARIMA(23, 1, 24) фактически учитывает только последние 24 значения для предсказания текущего значения временного ряда. Это означает, что модель будет использовать последние сутки как основу для предсказания будущего часового значения температуры. Однако, несмотря на то что предсказания модели учитывают колебания, они недостаточно точно отражают валидационные данные и представлены в виде волнообразной линии. В этой связи можно предположить, что модель, не учитывающая годовую сезонную компоненту, не смогла воспроизвести валидационные данные в должной мере и, следовательно, необходимо пересмотреть параметры модели с учётом сильного влияния годовой сезонности на предсказания.
Но сперва посмотрим на результаты обученной модели ARIMA(23, 1, 24), воспользовавшись удобными методами библиотеки statsmodels — summary() и plot.diagnostics() и проанализируем их. (Для этого я отдельно обучу модель ARIMA(23, 1, 24) и между делом зафиксирую, что для её выполнения потребуется 14-14.2 ГБ памяти согласно мониторингу TPU.)
Итак, выведем итоговую информацию обученной модели с помощью метода summary(). Для улучшения читаемости воспользуемся функцией print(), чтобы отформатировать вывод и придать ему структурированный вид.
print(arima_model.summary())Поскольку итоговая таблица на загляденье удобна и информативна, то ниже привожу описание её основных разделов, каждый из которых содержит важную информацию о модели и её параметрах. Для улучшения восприятия разделы таблицы выделены разным цветом.
Основные разделы, на которые можно подразделить вывод, следующие:
В итоговой таблице «поджабным» цветом подсвечены коэффициенты, p-значение которых меньше уровня значимости 0.05. Следовательно, для них мы отклоняем нулевую гипотезу, принимая, что эти коэффициенты статистически значимы. В связи с тем, что для процесса АР наиболее значимыми оказались шаги запаздывания 1 и 23, а для процесса СС - 22, 23 и 24, то это может означать, что во временном ряду с температурой имеется не только сильная сезонная структура, но и то, что текущее значение температуры зависит от предшествующих ему значений днём ранее. При этом надо признать, что наша модель ARIMA(23, 1, 24) смогла выделить эти значимые характеристики временного ряда и можно предположить, что она адекватно описывает данные (вспомните, мы приняли такие значения, чтобы учесть краткосрочные зависимости). Однако не будем спешить с выводами и посмотрим на результаты дополнительно выполненных статистических тестов, представленных в итоговой таблице.
В дополнение к итоговой таблице выполним графическую диагностику остатков модели на нормальность распределения и автокорреляцию.
arima_model.plot_diagnostics()
pyplot.show()
По графику стандартизированных остатков (z-оценка) видно, что имеется множество как положительных, так и отрицательных значений, которые значительно отклоняются от нуля как среднего значения. Эти значения модель принимает за так называемые выбросы.
На графике «квантиль-квантиль» представлено сравнение квантиля распределения остатков модели с квантилем теоретического нормального распределения в виде красной линии. Видно, что половине точек до лампочки эта «красная линия»: если в середине графика точки располагаются вдоль линии нормального распределения, то в крайних точках они сильно загибаются (имеются так называемые тяжёлые хвосты). Данное поведение говорит о том, что остатки модели избыточно рассеяны относительно среднего значения по сравнению с нормальным распределением и в них много выбросов [ссылка]. В целом, график «квантиль-квантиль» ясно показывает, что точки квантилей не лежат на теоретической нормальной кривой и, следовательно, данные распределены ненормально.
Это подтверждает гистограмма остатков: вытянутая пика и немного узкое распределение остатков говорят об отклонении их распределения от теоретического нормального. При этом на графике с гистограммой не видно, но по результатам вышерассмотренных тестов мы должны помнить, что имеется как отрицательная асимметрия, так и неоднородная дисперсия остатков.
Коррелограмма (график с АКФ) демонстрирует «быстрый спад» автокорреляционной функции и отсутствие значительной автокорреляции остатков; на графике нет значительных пик, которые бы указывали на необходимость пересмотра параметров модели (например, увеличения порядка разности параметра d или шагов процесса скользящего среднего).
Впечатляет, правда?! Всего две строчки кода, не считая pyplot’а — и мы получили довольно полную диагностику модели! Единственное, чем хотелось бы расширить графическую диагностику — графиком с ЧАКФ и коробчатой диаграммой для остатков модели. Давайте выполним их отдельно.
def dop_plot_diagnostics(lags=None):
"""ЧАКФ и коробчатая диаграмма остатков модели"""
fig, axs = pyplot.subplots(2, 1)
plot_pacf(residuals, lags=lags, ax=axs[0])
axs[0].set_title('Частная автокорреляционная функция')
axs[1].boxplot(residuals, vert=False)
axs[1].set_title('Коробчатая диаграмма')
pyplot.suptitle('Дополнительная графическая диагностика')
pyplot.show()
dop_plot_diagnostics(lags=48)
Пускай и не так заметно, но на графике с ЧАКФ выделяется 24-й шаг запаздывания, что указывает на необходимость пересмотра параметров модели применительно к увеличению шагов процесса авторегрессии. Согласно коробчатой диаграмме в остатках модели присутствует множество значительно удалённых точек (значений), преимущественно отрицательных, обуславливающих отрицательную асимметрию. Принимать эти значения за выбросы или нет, держа в уме тот аспект, что анализируется временной ряд с температурой, — это задача, требующая дополнительного изучения и применения специальных критериев. В данном исследовании мы не будем преобразовывать эти данные, а оставим как есть, считая, что они содержат важную информацию о динамике изменения температуры со временем.
Таким образом, на основании вышеприведённых элементов оценки параметров модели и диагностической проверки её остатков можно сделать вывод, что модель ARIMA(23, 1, 24), сформулированная в качестве отправной точки, неэффективно подогнана к обучающим данным. Есть признаки автокорреляции и гетероскедастичности в остатках, а также их несоответствие нормальному распределению. Поэтому остатки модели не ведут себя как белый шум. Следовательно, необходимо пересмотреть модель.
Что ж, если приведение распределения данных к нормальному мы вынесли за рамки нашего исследования и работаем с тем, что есть, то давайте попробуем избавиться от автокорреляции в остатках. Для этого воспользуемся подсказкой графика с ЧАКФ остатков предыдущей модели и увеличим количество шагов запаздывания процесса авторегрессии до 24. Затем обучим модель ARIMA(24, 1, 24) (потребуется 14.3 ГБ памяти согласно мониторингу TPU) и проверим, лишена ли эта модель автокорреляции в остатках или нет.
По итоговой таблице видно, что коэффициент 24-го шага запаздывания процесса авторегрессии имеет низкое значение стандартной ошибки и статистически значим. Его применение в новой модели немного улучшило все три информационных критерия, а полученное p-значение для теста Льюнга-Бокса (0.26), которое больше уровня значимости 0.05, позволяет сделать вывод, что в остатках модели отсутствует автокорреляция. Таким образом, с допущениями мы можем говорить, что модель ARIMA(24, 1, 24) адекватно описывает данные. Однако необходимо помнить о результатах выполнения предсказаний на валидационных данных, которые показывали, что предыдущая модель плохо справлялась с их воспроизводством. Поэтому нам следует сконцентрироваться на внесении глобальных изменений в модель, а не заниматься тонкой настройкой параметров модели после диагностической проверки.
Ничего не поделаешь, давайте приступим к определению сезонного порядка (P, D, Q, m) для модели SARIMAX... Но уже в самом начале нас поджидает неприятность в виде длительной отрисовки графиков АКФ и ЧАКФ с количеством шагов запаздывания, равным 8760*2 и более. (Прошу поберечь помидоры: хоть я и держал интригу до последнего, планируя представить хотя бы код для выполнения расчётов с данным годовым периодом, но длительность отрисовки (более 3-х часов ожидания) вконец «добила» эту затею.) Поэтому я настойчиво предлагаю не издеваться над компьютером, выполняя расчёты для сезонного параметра m=8760, и ограничиться теоретическими правилами (алгоритмом) нахождения параметров, на изучение и систематизацию которых я вложил немало времени.
Не берусь останавливать отчаянных последователей, которые захотят испытать свой компьютер; только предупрежу, что m=200 и более считается непосильной задачей для любого компьютера [ссылка] из-за обработки и временного хранения больших массивов данных (форма состояния пространства и применение фильтра Калмана для оценивания модели, матрицы ковариаций и пр.), требующих значительных ресурсов. Эта проблема известна, и разработчики библиотеки statsmodels предлагают решения для сокращения потребления памяти [ссылка], допуская различные ограничения, чтобы ускорить обработку. Но даже с этими настройками наш с вами сезонный годовой период в 8760 шагов не по зубам массовым компьютерам. Однако не стоит опускать руки: дальше мы применим общеизвестный метод, который позволяет моделировать сезонность любой длины и их любое количество, но об этом сильно далее и в своё время.
Итак, алгоритм определения порядков модели SARIMAX(p, d, q)(P, D, Q, m) сводится к следующему. Сначала начинают с определения параметров сезонного порядка. И первым делом, как по методологии АПРСС, рассмотренной в первой части главы 2, идентифицируют сезонный порядок разности D, чтобы устранить сезонность в данных. Поскольку у временного ряда температуры имеется явная сезонность, то выполняя дифференцирование временного ряда на период, равным годовой сезонности (8760), мы тем самым определяем D=1. Далее определим несезонный порядок разности d. Поскольку D>0, то выполним тест КПСС на дифференцированном (период=8760) временном ряду. Так как годовая сезонность устраняется с помощью параметра разности D, то тест КПСС проверяет стационарность временного ряда относительно линейной тенденции (regression='ct') и, если она имеется, устраняем её с помощью параметра d=1. При этом часто отмечается, что сумма порядков разности d и D не должна превышать значение 2.
Давайте запустим тест КПСС для определения значения параметра d:
def seasonal_trend_kpss_verification(series, periods=None):
"""Проверка ряда с сезонной разностью на стационарность вокруг тенденции"""
if periods:
name = series.name
series = series.diff(periods=periods).dropna()
series.rename(name + ' diff=' + str(periods), inplace=True)
msg = "Тест КФШШ для проверки стационарности вокруг тенденции (ct):"
print('n' + msg)
print('-' * len(msg))
kpss_result = kpss(series, regression='ct')
if kpss_result[1] > 0.05:
h_text = colored('стационарен', 'grey', 'on_green')
print(f"Ряд {series.name} {h_text}: p-value = {kpss_result[1]:.5f}")
else:
h_text = colored('нестационарен', 'grey', 'on_red')
print(f"Ряд {series.name} {h_text}: p-value = {kpss_result[1]:.5f} < 0.05")
seasonal_trend_kpss_verification(temperature, periods=8760)Результат представлен ниже.
По результатам теста видно, что предыдущее взятие сезонной разности не устранило тенденцию в данных. Необходимо устранить её с помощью параметра разности d, установив ему значение 1. Следовательно, мы принимаем параметры разности равными D=1 и d=1.
Для определения оставшихся параметров воспользуемся следующей функцией, которая построит график временного ряда температуры с разностью в годовую сезонность (период=8760), а также графики с АКФ и ЧАКФ с временным запаздыванием, равным 240 шагам. Такое количество шагов не требует длительного времени построения графиков и позволит понять «природу» как сезонных порядков P, Q, так и несезонных p, q. При этом я дополнительно отмечу на графике возможное поведение частной автокорреляционной функции на шагах 8760, 17520 для понимания того, как определяются сезонные параметры P и Q.
def plot_seasonal_order_analysis(series, periods=8760, lags=None):
"""График для определения параметров (P, Q), (p, q)"""
fig, axs = pyplot.subplots(3, 1)
series = series.diff(periods=periods).dropna()
axs[0].plot(series)
axs[0].set_title(f"Ряд с разностью {periods}-го порядка")
plot_acf(series, ax=axs[1], lags=lags)
axs[1].set_title(f"АКФ ({periods}-й порядок)")
plot_pacf(series, ax=axs[2], lags=lags)
axs[2].set_title(f"ЧАКФ ({periods}-й порядок)")
pyplot.show()
plot_seasonal_order_analysis(temperature, lags=240)
Итак, согласно этим графикам и таблице 3.2 (из книги Дж.Бокса и Г.Дженкинса) мы можем охарактеризовать ряд с разностью в годовой порядок как процесс авторегрессии: автокорреляционная функция бесконечна и представлена в виде затухающей экспоненты, а частная автокорреляционная функция конечна и обрывается после 2-х шагов запаздывания с появлением всплесков на шагах 8760 и 17520. Следовательно, сезонный параметр авторегрессии P может иметь значения 2 и 1 (по всплескам). Сезонный параметр скользящего среднего Q принимаем равным 0: никакие предыдущие сезонные остатки в модели не учитываются.
На том же графике с ЧАКФ видно, что два начальных шага представлены в виде обособленной пики: эти шаги (2 и 1) мы примем в качестве значений параметра авторегрессии p. Поскольку ранее мы обнаружили, что процесс скользящего среднего оказывает влияние на результат предсказаний, то примем q равным 1.
Таким образом, по результатам анализа вышерассмотренных графиков АКФ и ЧАКФ для временного ряда температуры с сезонной годовой разницей в 8760 шагов мы определили следующие возможные конфигурации модели SARIMAX:
(2-1, 1, 1)(2-1, 1, 0, 8760)Интерпретация данной модели следующая. Текущее часовое значение температуры зависит от двух (либо одного - в зависимости от установленного значения p) предыдущих часов, учитывая один предыдущий остаток (ошибку), и от двух (либо одного - в зависимости от установленного значения P) значений температуры в аналогичные часы из позапрошлого (или прошлого) года.
Не будь мы ограничены в вычислительных ресурсах, то проверить как сильно влияет годовая сезонность на прогнозирование не составило бы труда: «нажми на кнопку — получишь результат».
И тем не менее существуют известные методы, среди которых анализ Фурье, являющийся частью спектрального анализа, которые позволяют обойти ограничение длины сезонной компоненты. Это направление очень интересное и заслуживает отдельного рассмотрения в третьей, завершающей части второй главы. А пока разработаем модель SARIMAX с «адекватным» сезонным параметром применительно к вычислительным ресурсам, значение которого установим равным дневной сезонности (m=24). Тем самым мы учтём кратковременные зависимости с помощью сезонного параметра. А годовую сезонность попробуем учесть с помощью имеющихся экзогенных переменных.
Чтобы усилить интерес, предлагаю подбор оптимальной конфигурации модели SARIMAX выполнить с помощью мощного инструмента под названием авто-АРПСС. Наиболее известной python-библиотекой, предоставляющей возможность автоматического перебора параметров модели АРПСС, является pmdarima [ссылка]. Однако в этом исследовании я предлагаю воспользоваться библиотекой под названием statsforecast, которая, как заявляют её разработчики, в 20 раза быстрее, чем pmdarima (скомпилирована для выполнения высокопроизводительных расчётов с использованием numba [ссылка]).
Методы с автоматическим перебором параметров АРПСС обеих библиотек заимствованы из пакета forecast языка программирования R. Этот пакет, в свою очередь, опирается на систематизацию более ранних исследований и решений, направленных на автоматизацию процесса перебора с использованием программного обеспечения TRAMO, SEATS, Autobox и других подобных инструментов. Об этом справедливо упоминают авторы статьи «Автоматическое прогнозирование временных рядов: пакет forecast для языка R» (2008) [ссылка].
Рекомендую ознакомиться с этой работой, чтобы уверенно ориентироваться в разделах с инструментом auto-ARIMA различных python-библиотек. Я лишь кратко отмечу отличия от «ручного» метода определения параметров, который мы с вами применяли ранее для определения несезонных и сезонных порядков модели АРПСС/САРПСС.
Итак, если выбор как несезонных, так и сезонных параметров процессов АР и СС (p, q, P и Q) осуществляется путём минимизации информационного критерия Акаике (AIC) (либо информационного критерия Байеса, BIC), то определение параметров разности выполняется с помощью различных статистических тестов. В отличие от python, в языке R эти тесты давно реализованы и широко используются. В частности, auto.arima() для определения параметра сезонной разницы D использует функцию nsdiffs(), которая охватывает несколько тестов [ссылка]. После того, как параметр D установлен, с помощью другой функции — ndiffs() - определяют параметр d. Данная функция использует тест КПСС по умолчанию [ссылка] и применяется к ряду либо с сезонной разницей (если D=1), либо к исходному (если D=0).
С помощью следующего кода мы запустим автоматический перебор параметров для поиска наилучшей модели SARIMA с длиной сезона в 24 шага (season_length=24):
from statsforecast.models import AutoARIMA
train_data = temperature['2009-01-01':'2011-12-31'].copy()
val_data = temperature['2012-01-01':'2012-01-07'].copy()
model = AutoARIMA(season_length=24, seasonal=True, trace=True)
model = model.fit(y=train_data.to_numpy())Результаты представлены ниже.
ARIMA(2,0,2)(1,1,1)[24] :inf
ARIMA(0,0,0)(0,1,0)[24] :139845.93882574368
ARIMA(1,0,0)(1,1,0)[24] :62793.959059700705
ARIMA(0,0,1)(0,1,1)[24] :107428.13870839613
ARIMA(1,0,0)(0,1,0)[24] :68380.80916728784
ARIMA(1,0,0)(2,1,0)[24] :60514.528469509656
ARIMA(1,0,0)(2,1,1)[24] :inf
ARIMA(1,0,0)(1,1,1)[24] :inf
ARIMA(0,0,0)(2,1,0)[24] :138470.59089873798
ARIMA(2,0,0)(2,1,0)[24] :54144.152318747736
ARIMA(2,0,0)(1,1,0)[24] :56682.81810435849
ARIMA(2,0,0)(2,1,1)[24] :inf
ARIMA(2,0,0)(1,1,1)[24] :inf
ARIMA(3,0,0)(2,1,0)[24] :inf
ARIMA(2,0,1)(2,1,0)[24] :54143.24675833847
ARIMA(2,0,1)(1,1,0)[24] :56682.01262490523
ARIMA(2,0,1)(2,1,1)[24] :inf
ARIMA(2,0,1)(1,1,1)[24] :inf
ARIMA(1,0,1)(2,1,0)[24] :55008.69589543988
ARIMA(3,0,1)(2,1,0)[24] :54102.41679797225
ARIMA(3,0,1)(1,1,0)[24] :56643.717496704645
ARIMA(3,0,1)(2,1,1)[24] :inf
ARIMA(3,0,1)(1,1,1)[24] :inf
ARIMA(4,0,1)(2,1,0)[24] :54055.97547555753
ARIMA(4,0,1)(1,1,0)[24] :56595.9640338018
ARIMA(4,0,1)(2,1,1)[24] :inf
ARIMA(4,0,1)(1,1,1)[24] :inf
ARIMA(4,0,0)(2,1,0)[24] :54056.11095450652
ARIMA(5,0,1)(2,1,0)[24] :53941.59162502506
ARIMA(5,0,1)(1,1,0)[24] :56441.562388514416
ARIMA(5,0,1)(2,1,1)[24] :inf
ARIMA(5,0,1)(1,1,1)[24] :139861.94553157385
ARIMA(5,0,0)(2,1,0)[24] :inf
ARIMA(5,0,2)(2,1,0)[24] :inf
ARIMA(4,0,2)(2,1,0)[24] :54057.93584322888
Now re-fitting the best model(s) without approximations...
ARIMA(5,0,1)(2,1,0)[24] :53941.59162502506 В результате автоматического перебора параметров с помощью AutoARIMA лучшей моделью по минимизации AIC оказалась модель SARIMA(5, 0, 1)(2, 1, 0, 24). Согласно документации [ссылка] параметры авторегрессии p и P достигли своих максимальных значений, установленных по умолчанию, — max_p=5 и max_P=2 соответственно. Данное обстоятельство должно подвергнуть критике подобранную модель: в наших данных действительно имеется сильный авторегрессионный процесс, но ограничивается ли он теми максимальными значениями по авто-АРПСС?
Постройте для временного ряда температуры графики с АКФ и ЧАКФ с сезонной разностью m=24, а затем проверьте, можно ли включить большее число шагов запаздывания в процессы авторегрессии. Возможно, кто-нибудь из вас пойдёт дальше — и обучит модель SARIMA подобранной (возможно, более подходящей) конфигурации с помощью библиотеки statsmodels на TPU google.colab (потребуется более 20 ГБ памяти) и проанализирует её результаты. Делитесь ими в комментариях, с интересом обсудим!
А мы тем временем попробуем учесть годовую сезонность с помощью имеющихся экзогенных переменных.
Напоминаю, экзогенные переменные используются для учёта внешних воздействий (внешних факторов) с целью более точного и адекватного описания поведения эндогенного временного ряда. В моделях ARIMAX/SARIMAX экзогенные переменные рассматриваются как ковариаты, которые добавляются в модель наряду с эндогенной переменной. В отличие от эндогенной переменной к ним не применяют операторы разности, и они не подвергаются преобразованию, а используются в исходном виде. В этом можно убедиться по формулам моделей ARIMAX/SARIMAX [ссылка].
Предлагаю провести интересный эксперимент, в котором одна и та же сезонная модель будет тестироваться с разной экзогенной переменной. В первом случае будем использовать температуру в исходном виде (без преобразований), а во втором — температуру с разностью первого порядка (как мы обнаружили ранее, временной ряд с температурой нестационарен по тесту КПСС). Выполняя предсказания на валидационных данных, мы проверим, какая их этих двух экзогенных переменных способна улучшить результат и, следовательно, наиболее значима. На основании полученных результатов сделаем вывод о целесообразности применения оператора разности к нестационарной экзогенной переменной.
С помощью нижеприведённого кода мы сначала выделим обучающую и валидационную части для эндогенного параметра модели, определим отдельную переменную с дифференцированной температурой, а затем на их основе сформируем экзогенный параметр, который будет состоять из двух вышерассмотренных элементов (столбцов). Далее для каждого экзогенного столбца обучим одну и ту же модель SARIMAX, выполним предсказания по меткам валидационных данных на 168 шагов с учётом валидационной части экзогенной переменной, рассчитаем ошибку mape и сохраним результаты в словарь results. Этот словарь как результат функции build_sarimax_model() вложим в общий для обеих сезонных моделей словарь sarimax_sum_results, по ключам и элементам которого с помощью функции show_sarimax_temp_plot_prediction() будем оценивать эффективность выполнения предсказания моделей.
Так как при использовании экзогенных переменных в модели SARIMAX библиотеки statsmodels возникает ошибка при обучении (экзогенные параметры фактически игнорируются), мы будем применять модель ARIMA. С учётом использования сезонного порядка и экзогенных переменных эта модель будет эквивалентна SARIMAX.
#from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.arima.model import ARIMA
horizont = 168
endog_train = temperature['2009-01-01':'2011-12-31'].copy()
endog_val = temperature['2012-01-01':'2012-01-07'].copy()
temperature_diff = temperature.diff().fillna(0)
exog_train = endog_train.to_frame().copy()
exog_val = endog_val.to_frame().copy()
exog_train[f"{exog_train.columns[0]}_diff"] = temperature_diff['2009-01-01':'2011-12-31'].copy()
exog_val[f"{exog_val.columns[0]}_diff"] = temperature_diff['2012-01-01':'2012-01-07'].copy()
sarimax_sum_results = {}
def build_sarimax_model(colls):
"""
Обучение модели SARIMAX с экзогенными переменными и
сохранение её результатов
"""
order=(5, 0, 1)
seasonal_order=(2, 1, 0, 24)
results = {}
sarimax_model = ARIMA(endog=endog_train,
exog=exog_train[colls],
order=order,
seasonal_order=seasonal_order,
trend=None,
enforce_stationarity=False,
enforce_invertibility=False)
sarimax_model = sarimax_model.fit()
preds = sarimax_model.get_prediction(start=endog_val.index.min(),
end=endog_val.index.max(),
exog=exog_val[colls])
mape = mean_absolute_percentage_error(y_true=endog_val,
y_pred=preds.predicted_mean)
results['AIC'] = sarimax_model.aic
results['BIC'] = sarimax_model.bic
results['residuals'] = sarimax_model.resid
results['preds_mean'] = preds.predicted_mean
results['MAPE'] = mape
return results
for exoq_col in exog_train:
key = f"'{exoq_col}'"
print('n' f"Обучение модели SARIMAX с экзоген.переменной {key}")
sarimax_sum_results['SARIMAX[24] ' + key] = build_sarimax_model(exoq_col)После завершения расчётов выведем результаты обеих моделей на общем графике с помощью следующей функции:
def show_sarimax_temp_plot_prediction(model_dict):
"""График сравнения истинных значений с предсказанными"""
pyplot.figure()
pyplot.plot(endog_val, label='Валидационные данные')
for key, model_data in model_dict.items():
predictions = model_data['preds_mean']
label_model = f"{key} "
label_mape = f"[mape: {model_data['MAPE']*100:.2f}%] "
label_aic = f"[AIC: {int(model_data['AIC'])}]"
pyplot.plot(predictions, dashes=[4, 3],
label=(label_model + label_mape + label_aic))
pyplot.xlabel('Дата наблюдения')
pyplot.ylabel('Температура (град C)')
pyplot.title('Предсказания на валидационных данных: сравнение exoq и exoq_diff')
pyplot.grid(True)
pyplot.legend(loc='best')
pyplot.show()
show_sarimax_temp_plot_prediction(model_dict=sarimax_sum_results)
График показывает, что модель SARIMAX с нестационарной экзогенной переменной идеально точно воспроизвела валидационную выборку, ошибка mape равна нулю. Напротив, модель SARIMAX с предварительно продифференцированной экзогенной переменной показала значительно худшие результаты, ошибка mape подскочила до 56%.
Таким образом, предсказания на валидационных данных показали, что предварительная обработка экзогенной переменной привела к снижению точности. Это подтверждает, что экзогенные переменные следует использовать в исходном виде.
Давайте рассмотрим коэффициенты обеих экзогенных переменных, которые мы получим благодаря дополнительно выполненной модели SARIMAX.
Оба коэффициента имеют p-значение, равное 0, что подтверждает их статистическую значимость при любых стандартных уровнях значимости. Коэффициент первого параметра, выделенный зелёным цветом, равен 1, что указывает на его положительное влияние на модель, которое подтверждается высокой точностью и уверенностью в значениях, поскольку доверительный интервал совпадает в верхней и нижней границах и равен 1. Как отмечается в литературе по прогнозированию, это делает его интерпретацию «интуитивно понятной и предпочтительной».
Коэффициент второго параметра обладает отрицательным «околонулевым» значением (если быть точным, его значение в стандартной форме равно -0.00000188), что говорит о его незначительном отрицательном влиянии, которое подтверждается менее точным доверительным интервалом с отрицательными границами, что менее удобно для интерпретации.
Какому из этих коэффициентов вы бы отдали предпочтение? Вероятнее всего первому из-за его положительного влияния и понятной интерпретации. Таким образом, рассмотрение коэффициентов также подтверждает, что экзогенный фактор лучше оставлять без изменений.
А теперь давайте ещё раз посмотрим на график, на котором модель SARIMAX с нестационарной экзогенной переменной идеально точно воспроизвела валидационную выборку. Использование экзогенных параметров в прогнозных моделях накладывает важное ограничение, связанное с необходимостью предоставления соответствующей экзогенной информации во время выполнения предсказаний или прогнозирования. Другими словами, применение экзогенных данных подразумевает, что для выполнения предсказаний или прогнозов с помощью обученной модели также необходимо использовать соответствующие экзогенные данные (валидационные или тестовые). Это важно учитывать.
А если у нас нет валидационных экзогенных данных и мы «просто хотим выполнить прогноз»? Для решения этой проблемы наиболее распространённым подходом является использование последних n-строк обучающих экзогенных данных. Однако нужно учитывать, что в этом случае нарушается временная структура (последние наблюдения могут сильно отличаться от предыдущих) и это может привести к снижению точности предсказаний (прогнозов).
В этом легко убедиться благодаря следующему эксперименту. Обратите внимание, мы будем использовать несколько временных рядов в качестве экзогенных параметров (переменная colls), а их валидационная часть будет соответствовать последним 168 наблюдениям обучающей части (переменная exog_val).
from statsmodels.tsa.arima.model import ARIMA
horizont = 168
colls = [
'p (mbar)',
'T (degC)',
'VPmax (mbar)',
'sh (g/kg)',
'rho (g/m**3)',
'wv (m/s)'
]
endog_train = temperature['2009-01-01':'2011-12-31'].copy()
endog_val = temperature['2012-01-01':'2012-01-07'].copy()
exog_train = df[colls].loc['2009-01-01':'2011-12-31'].copy()
exog_val = exog_train.iloc[-horizont:, :].copy()
order=(5, 0, 1)
seasonal_order=(2, 1, 0, 24)
sarimax_model = ARIMA(endog=endog_train,
exog=exog_train,
order=order,
seasonal_order=seasonal_order,
trend=None,
enforce_stationarity=False,
enforce_invertibility=False)
sarimax_model = sarimax_model.fit()
preds = sarimax_model.get_prediction(start=endog_val.index.min(),
end=endog_val.index.max(),
exog=exog_val)
mape = mean_absolute_percentage_error(y_true=endog_val,
y_pred=preds.predicted_mean)
def show_sarimax_plot_exog_train_prediction():
"""График сравнения истинных значений с предсказанными"""
pyplot.figure()
pyplot.plot(endog_val, label='Валидационные данные')
label_mape = f" [mape: {mape*100:.2f}%]"
pyplot.plot(preds.predicted_mean, dashes=[4, 3],
label='Предсказанные данные '+ label_mape)
pyplot.xlabel('Дата наблюдения')
pyplot.ylabel('Температура (град C)')
pyplot.title('Предсказания на валидационных данных с exog_train[-history:, :]')
pyplot.grid(True)
pyplot.legend(loc='best')
pyplot.show()
show_sarimax_plot_exog_train_prediction()
График показывает, что применение данного подхода позволило отразить годовые изменения температуры, однако точность предсказаний существенно снизилась: ошибка mape на валидационных данных увеличилась до 78%.
Существуют ли методы, которые позволили бы нам улучшить точность предсказаний? Да, такие методы существуют, и одним из них является преобразование Фурье. В третьей, завершающей, части мы с вами, дорогие друзья, перейдём к более глубокому анализу целевого временного ряда - и применим некоторые методы спектрального анализа, чтобы точнее предсказать будущие значения температуры, минимизируя влияние шумов и случайных колебаний.
Автор: OLZ1
