В первой главе исследования был описан набор данных с временными рядами о погоде, который мы будем использовать для выполнения задачи прогнозирования температуры, а также были приведены шаги по его предварительной подготовке.

В данной главе мы рассмотрим процессы авторегрессии-проинтегрированного скользящего среднего по методологии АРПСС (в англоязычной терминологии - ARIMA). Разберёмся, почему процесс АРПСС позволяет получить широкий класс стационарных и нестационарных моделей, которые адекватно описывают многие встречающиеся на практике временные ряды. А затем применим эту методологию с целью нахождения подходящего подкласса моделей из общего семейства моделей АРПСС для адекватного прогнозирования будущих значений температуры.

1. Описание набора климатических данных и его предварительная подготовка.

2. Статистический анализ данных, включая оценку сезонной компоненты, для разработки моделей из семейства ARIMA с целью прогнозирования температуры.

3. Разработка моделей глубокого обучения Keras с использованием слоя LSTM для аналогичной задачи.

4. Сравнение и оценка эффективности результатов прогнозирования.

5. Рассмотрение подхода прогнозирования смещённой температуры с помощью модели глубокого обучения Keras с использованием слоёв LSTM.

Методы машинного обучения и глубокого обучения в целом часто рассматриваются как универсальные средства для решения различных задач, включая прогнозное моделирование, обработку естественного языка и распознавание изображений. А искусственный интеллект так и вовсе воспринимается как волшебные бубенцы, звон которых обеспечивает моментальное решение всевозможных проблем. Однако в направлении прогнозирования временных рядов, как специфического направления прогнозного моделирования, ситуация обстоит не так однозначно. Далеко не так однозначно.

Ещё в 2018 году было опубликовано исследование [ссылка], в котором оценивалась и сравнивалась эффективность как традиционных статистических методов прогнозирования, так и современных методов машинного обучения, включая рекуррентные нейронные сети, для решения задач прогнозирования временных рядов. В основу экспериментов было заложено выполнение прогнозирования как на один шаг вперёд, так и на несколько интервалов (горизонтов) разной протяжённости. Авторы исследования констатируют: «Сравнив поствыборочную точность популярных методов ML [Machine Learning. — Примечание] с точностью восьми традиционных статистических методов [среди которых ARIMA. — Примечание], мы обнаружили, что последние доминируют по обоим показателям точности [sMAPE и MASE. — Примечание] и для всех рассмотренных горизонтов прогнозирования. Кроме того, мы отметили, что их вычислительные требования значительно выше, чем у статистических методов.»

Также авторы обозначили проблемы использования прогнозов, выполненных с помощью моделей машинного обучения, которые вызывают множество сомнений в их качестве и интерпретации: «Получение цифр из "чёрного ящика" неприемлемо для практиков, которым необходимо знать, как возникают прогнозы и как на них можно повлиять или скорректировать, чтобы получить работоспособные предсказания».

Итоговый посыл приводимого исследования можно заключить в необходимость использования традиционных статистических методов в качестве базового решения проблемы прогнозирования временных рядов: «Проблема академической литературы по ML‑прогнозированию заключается в том, что большинство опубликованных исследований предоставляют прогнозы и заявляют об удовлетворительной точности, не сравнивая их с простыми статистическими методами или даже упрощёнными эталонами [яркий тому пример — kaggle. — Примечание]. Это порождает ожидания, что методы ML обеспечивают точные прогнозы, но без каких‑либо эмпирических доказательств того, что это так.»

Исследование очень интересное; рекомендую ознакомиться с ним. Опираясь на его выводы применительно к нашей задаче прогнозирования температуры мы зададим очень высокую планку, применяя в качестве базового решения статистическую модель SARIMAX, и постараемся допрыгнуть до её результатов с помощью модели глубокого обучения Keras на базе LSTM. И попытаемся выяснить — так ли уж всё плохо с глубоким обучением в решении задач прогнозирования временных рядов.

Итак, начинаем.

В начале проясним следующее.

Модель Бокса-Дженкинса, или ARIMA (от англ. AutoRegressive Integrated Moving Average, или Авторегрессионная проинтегрированная скользящая средняя (АРПСС)) — это в первую очередь методология, а не одна определённая модель. Давайте разбираться.

ARIMA состоит из трёх основных компонентов:

• AutoRegressive (AR): Эта часть модели отвечает за зависимость текущего значения временного ряда от его предыдущих значений. Это означает, что текущее значение ряда может быть описано как линейная комбинация его прошлых значений.

• Integrated (I): Эта часть модели связана с преобразованием временного ряда для достижения стационарности. Если временной ряд не стационарен (например, имеет тенденцию или сезонные колебания), его можно сделать стационарным путём преобразования — в основном с помощью операций взятия последовательной разности (дифференцирования) или логарифмирования [ссылка].

• Moving Average (MA): Эта часть модели учитывает влияние случайных ошибок (или «шумов») на текущее значение временного ряда. Это означает, что текущее значение зависит от предыдущих ошибок прогноза. Другими словами, модель скользящего среднего предполагает, что в ошибках модели в предшествующие периоды сосредоточена информация обо всей предыстории временного ряда.

Каждый вышерассмотренный компонент выражается в буквенном обозначении:

• p - порядок авторегрессии;

• d - порядок разности;

• q - порядок скользящего среднего.

Параметры p, d и q могут принимать различные значения, что приводит к созданию множества реализаций ARIMA (=различных моделей). Ниже приведены примеры конфигураций:

• AR (2, 0, 0) или AR (p=2):

  модель только с процессом авторегрессии с двумя предыдущими значениями (лагами).

• MA (0, 0, 1) или MA (q=1):

  модель только с процессом скользящего среднего с одним лагом.

• ARMA (1, 0, 1) или ARMA (p=1, q=1):

  модель со смешанным процессом авторегрессии-скользящего среднего: с одним лагом для авторегрессии и одним лагом для скользящего среднего.

• ARIMA (2, 1, 2)

  более сложная модель с комбинацией всех трёх компонентов: авторегрессия с двумя лагами, разность первого порядка и скользящее среднее с двумя лагами.

Конкретная конфигурация подбирается в ходе анализа характеристик рассматриваемого целевого временного ряда. Например, если рассматриваемый временной ряд по статистическим тестам уже является стационарным, то выполнять дифференцирование (или логарифмирование), как правило, не нужно: в этом случае параметр d принимается равным 0 (об этапах определения компонентов p, d и q будет рассказано далее). В этой связи стационарность временного ряда необходима для лучшего, более точного определения параметров ARMA (p, q).

Методология заключается в том, что построение модели ARIMA по реализации случайного процесса (если будущие значения временного ряда могут быть описаны только с помощью распределения вероятностей) её авторы — Дж.Бокс и Г.Дженкинс — предложили разбить на несколько этапов:

1. Идентифицировать порядок разности d, то есть добиться стационарности ряда, взяв достаточное количество последовательных разностей [из книги «Анализ временных рядов: прогноз и управление /Дж.Бокс и Г.Дженкинс. Выпуск 1. — Издательство «Мир», Москва 1974. — 194 с.»]. Для определения значения d может быть применён эвристический критерий.

2. Построить для полученного стационарного временного ряда АКФ [Автокорреляционную функцию. — Примечание] и ЧАКФ [Частную автокорреляционную функцию. — Примечание] [для определения тесноты статистической связи между запаздывающими значениями (лагами) временного ряда. — Примечание]. Исследуя характер их поведения, выдвинуть гипотезы о значениях параметров p и q, то есть подбирается модель ARMA (p, q). На данном этапе формируется базовый набор моделей, включающий одну, две или даже большее количество моделей. [из книги «Анализ временных рядов и прогнозирование: учебник /В.Н. Афанасьев, М.М. Юзбашев. Изд. 2-е, перераб. и доп. — М.: Финансы и статистика, 2010. — 290 с.»]

3. Оценить параметры модели и выполнить «диагностическую проверку», сопоставляя результаты модели с фактическими данными для выявления дефектов подгонки и диагностирования их причин [из книги «Анализ временных рядов: прогноз и управление /Дж.Бокс и Г.Дженкинс. Выпуск 1. — Издательство „Мир“, Москва 1974. — 36, 232, 313 с.»]. Используются графические методы и статистические тесты для анализа остатков модели (разницы между фактическими и предсказанными значениями), чтобы убедиться, что они ведут себя как белый шум и что модель адекватно описывает данные. При обнаружении проблем модель корректируется.

4. Прогнозирование. После определения точной модели и её подгонки к фактическим данным, её можно использовать для прогнозирования будущих значений целевого временного ряда [из книги «Анализ временных рядов: прогноз и управление /Дж.Бокс и Г.Дженкинс. Выпуск 1. — Издательство „Мир“, Москва 1974. — 144 с.»]. Прогнозы, основанные на содержащейся в модели информации, могут быть как краткосрочными, так и долгосрочными.

Таким образом, методология Бокса‑Дженкинса заключается в анализе временного ряда для определения оптимальных параметров модели прогнозирования (p, d, q) процессов ARIMA. В свою очередь, модель ARIMA — это модель из семейства (или классов) моделей ARIMA, в которой к временному ряду применялась операция взятия разности (I) хотя бы один раз, чтобы сделать его стационарным, и объединены условия процессов авторегрессии (AR) и скользящей средней (MA).

Также добавим, что существует расширение модели ARIMA под названием SARIMA. Это модель с сезонной составляющей (S), что добавляет ещё больше возможностей для моделирования временных рядов.

В завершении определим что такое «эндогенные» и «экзогенные» переменные, представив определения из книги [Практический анализ временных рядов: прогнозирование со статистикой и машинное обучение.: Пер. с англ. — СПб.: ООО «Диалектика», 2021. — 238 с.]: «С точки зрения статистических исследований переменные модели, которые оказывают влияние друг на друга, называются эндогенными — их значения описываются поведением модели. В противоположность им экзогенные переменные не находят объяснения в модели — они не описываются предположениями модели, а потому их значения принимаются такими, какими они есть, и не подвержены динамическим изменениям.»

Теперь, когда мы хорошенько прокачали знания в методологии ARIMA, предлагаю вернуться к нашим временным рядам и проанализировать их с помощью статистических методов, чтобы подобрать подходящие параметры p, d и q.

Сначала импортируем необходимые модули и атрибуты, которые мы будем использовать во второй главе исследования, а также изменим некоторые параметры метода pyplot библиотеки matplotlib по умолчанию на пользовательские:

from matplotlib import pyplot
from numpy import arange
from numpy.random import seed
from pandas import read_parquet, Series
from sklearn.metrics import mean_absolute_percentage_error
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tools.sm_exceptions import InterpolationWarning
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller, kpss
from termcolor import colored
from warnings import simplefilter

simplefilter('ignore', InterpolationWarning)

pyplot.rcParams["figure.autolayout"]=True
pyplot.rcParams["figure.figsize"]=(14, 7)

seed(1)

С помощью функции read_parquet() библиотеки pandas импортируем подготовленный в первой главе исследования набор данных, передав ей название файла вместе с его расширением. Поскольку имеется нерешённый вопрос (open issue на github) о потери информации с частотой данных при их сохранении в формат parquet, то дополнительно с помощью метода asfreq('h') переустановим почасовую частоту.

df = read_parquet('jena_climate_2009_2016_hour_grouped_data.parquet').asfreq('h')

Создадим отдельную переменную с целевым временным рядом для прогнозирования с типом объекта pandas Series. Напомню: объект Series «представляет собой одномерный однородный маркированный массив, содержащий значения и индекс» [Pandas в действии / Борис Пасхавер. — СПб.: Питер, 2023. — 89 С.].

temperature = df['T (degC)']

Итак, первое, что нам нужно сделать согласно вышеописанной методологии – это провести проверку имеющихся временных рядов на стационарность (определить параметр d). Сделать это можно визуально по графикам (эмпирически), но надёжнее всего прибегнуть к специально созданным для этого статистическим тестам. На практике часто встречается применение связки из двух тестов - расширенного теста Дики-Фуллера (в англоязычной терминологии — ADF‑test) и теста Квятковского‑Филлипса‑Шмидта‑Шина (KPSS‑test) [ссылка]. Выполнение сразу двух тестов гарантирует более тщательную проверку временного ряда на стационарность. Цель первого — проверить наличие единичного корня; вместе с тем этот тест имеет ряд существенных недостатков [ссылка], которых лишён тест КФШШ, целью которого является проверка наличия в данных тенденции и сезонности.

При этом методологический подход теста КПСС полностью отличается от подхода расширенного теста Дики-Фуллера, главное различие которого заключается в перестановке нулевой и альтернативной гипотез:

Обычно для обоих тестов уровень значимости принимается равным 0,05. Сравнивая с ним рассчитанные в результате выполнения тестов значения p определяют, является ли временной ряд стационарным согласно каждому тесту или нет.

Итак, выполним расширенный тест Дики-Фуллера и извлечём из его результатов только p-значение (result[1]), которое затем будет сравниваться с уровнем значимости для отклонения либо нулевой, либо альтернативной гипотезы. Для удобного восприятия вывода информации определившееся стационарное состояние временных рядов подсветим зелёным цветом, а нестационарное — красным. При этом, если в ходе тестирования обнаружатся нестационарные временные ряды, то выведем их названия.

def adf_test_columns_stationarity(dataframe):
    """Проверка временных рядов на стационарность по тесту Дики-Фуллера"""
    msg = "Выполнение расширенного теста Дики-Фуллера:"
    print('\n' + msg)
    print('-' * len(msg))
    adf_non_stationary_columns = []

    for i in range(len(dataframe.columns)):
        result = adfuller(dataframe[dataframe.columns[i]])
        if result[1] < 0.05:
            h_text = colored('стационарен', 'grey', 'on_green')
            print(f"Ряд {df.columns[i]} {h_text}: p-value = {result[1]:.5f}")
        else:
            h_text = colored('нестационарен', 'grey', 'on_red')
            print(f"Ряд {df.columns[i]} {h_text}: p-value = {result[1]:.5f} > 0.05")
            adf_non_stationary_columns.append(dataframe.columns[i])

    if len(adf_non_stationary_columns) > 0:
        msg = "Нестационарные временные ряды:"
        print('\n' + msg)
        print('-' * len(msg))
        print(* adf_non_stationary_columns, sep='\n')


adf_test_columns_stationarity(dataframe=df)

Результаты тестирования представлены ниже:

Как видно, по расширенному тесту Дики-Фуллера все временные ряды определены как стационарные. Следом запустим проверку временных рядов на стационарность по тесту КПСС:

def kpss_test_columns_stationarity(dataframe):
    """Проверка временных рядов на стационарность по тесту КФШШ"""
    msg = "Выполнение теста Квятковского-Филлипса-Шмидта-Шина:"
    print('\n' + msg)
    print('-' * len(msg))
    kpss_non_stationary_columns = []

    for i in range(len(dataframe.columns)):
        result = kpss(dataframe[dataframe.columns[i]])
        if result[1] > 0.05:
            h_text = colored('стационарен', 'grey', 'on_green')
            print(f"Ряд {df.columns[i]} {h_text}: p-value = {result[1]:.5f}")
        else:
            h_text = colored('нестационарен', 'grey', 'on_red')
            print(f"Ряд {df.columns[i]} {h_text}: p-value = {result[1]:.5f} < 0.05")
            kpss_non_stationary_columns.append(dataframe.columns[i])

    if len(kpss_non_stationary_columns) > 0:
        msg = "Нестационарные временные ряды:"
        print('\n' + msg)
        print('-' * len(msg))
        print(* kpss_non_stationary_columns, sep='\n')


kpss_test_columns_stationarity(dataframe=df)

Результаты тестирования представлены ниже:

В противовес расширенному тесту Дики-Фуллера тест КПСС выявил нестационарность у большей части временных рядов, включая целевой. Случай, когда расширенный тест Дики-Фуллера указывает на стационарность временного ряда, а тест КПСС на нестационарность коллеги из statsmodels [ссылка] охарактеризовали как «ряд является разностно-стационарным: чтобы сделать ряд стационарным, необходимо использовать дифференцирование». В нашем случае можно сделать вывод, что временной ряд не является строго стационарным из-за присутствия сезонности/тенденции.

Обращаю внимание, что по умолчанию тест КПСС библиотеки stastmodels выполняет проверку стационарности временного ряда по константе (regression='c'), то есть по среднему значению ряда. Если нулевая гипотеза о стационарности отвергается, то это может указывать на наличие сезонных колебаний. При regression='ct' тест будет проверять стационарность временного ряда относительно линейной тенденции [ссылка], благодаря чему можно предположить о наличие во временном ряду тенденции. В данном исследовании мы будем использовать тест КППС с настройками по умолчанию и подтверждение его альтернативной гипотезы будет означать, что ряд является нестационарным по причине присутствия сезонных колебаний.

Но давайте ещё раз посмотрим на результаты теста КПСС. Ничего не напоминает? Нет, я не про Коммунистическую Партию Советского Союза (как можно убедиться выше, понимание того, как работает тест КПСС не требует партбилета). («ха» три раза) Вспомните подготовленный график распределения временных рядов во времени из первой главы нашего исследования. Тест КПСС напрямую подтвердил наши догадки о периодичности данных у большинства временных рядов.

Глядя на получившиеся результаты, с нескрываемым удовольствием сообщаю, что наш с вами случай — самый интересный в виду распределения временных рядов по состоянию стационарности. Как и целевой временной ряд — температура, — так и большая часть остальных временных рядов, которые мы планировали использовать в качестве экзогенных, являются нестационарными по сезонности/тренду. В этой связи нам предстоит проделать множество экспериментов, по результатам которым мы попробуем ответить на некоторые вопросы, которые до сих пор не определены среди статистиков/исследователей данных как разрешённые. Например: в случае использования временных рядов в качестве экзогенных требуется ли для них предварительно выполнять операцию взятия разности, если по результатам тестов они оказались нестационарными? Забавную странность этому вопросу придаёт тот факт, что в одних книгах/статьях, предоставляемых авторами вопросов в качестве доводов на stats.stackexchange.com утверждается одно, а в других статьях по прогнозированию — обратное. Другой пример: подразумевает ли установка параметра d > 0 в параметрах подгонки модели SARIMAX автоматическое интегрирование не только эндогенного регрессора, но и экзогенных? Вопросы интересные и для их решения мы с вами задействуем мощное и проверенное средство — эксперимент, который, как известно, в нашем исследовательском деле всему голова: выполним серию экспериментов, а затем проанализируем получившиеся результаты.

Итак, в виду того, что наш случай не просто самый интересный, но и затратный на количество увлекательных экспериментов, то, как и решение любой сложной задачи, будет правильным разбить его на несколько шагов и идти по пути от самого простого решения к сложному. В этой связи для решения задачи прогнозирования температуры мы сначала рассмотрим только один временной ряд — собственно температуру, которая будет выступать в качестве эндогенной, то есть, как мы помним, описываемой моделью на основе исторической информации. В англоязычной терминологии этому соответствует понятие univariate time series forecasting (прогнозирование одномерных временных рядов). Как только мы подберём для этой модели оптимальные параметры p, d и q и получим хорошие результаты, мы затем попробуем улучшить их — сначала попытаемся добавить сезонную компоненту, а затем и экзогенные переменные и выясним, смогут ли они улучшить результат.

За дело!

Для удобного восприятия значений, полученных в ходе выполнения статистических тестов, с помощью следующей функции мы сократим количество знаков дробной части числа до num, которое по умолчанию равняется 3, а целые числа оставим без изменений:

def round_series_floats(series, num=3):
    """Округление вещественных чисел до num-знаков после точки"""
    return series.apply(lambda x: x if isinstance(x, int) else round(x, num))

Итак, выполним углублённую проверку временного ряда с температурой на стационарность и выведем все результаты обоих тестов. Дополнительно построим график, на котором изобразим оригинальный временной ряд вместе со скользящей средней и стандартным отклонением по годовому периоду, равному 8760:

def plot_movingAverage(timeseries, adf_p_value, kpss_p_value, window=24*365):
    """График для определения тенденции"""
    movingAverage = timeseries.rolling(window=window).mean().dropna()
    movingSTD = timeseries.rolling(window=window).std().dropna()

    x_start, x_end = movingAverage.index[0], movingAverage.index[-1]
    y_start, y_end = movingAverage.iloc[0], movingAverage.iloc[-1]

    pyplot.figure()
    pyplot.plot(timeseries, color='lightgrey', label='Оригинальный ряд')
    pyplot.plot(movingAverage, color='coral', label='Скользящее среднее')
    pyplot.plot(movingSTD, color='black', label='Стандартное отклонение')

    pyplot.scatter(x_start, y_start, color='lightsalmon', marker='o')
    pyplot.scatter(x_end, y_end, color='lightsalmon', marker='o')

    pyplot.text(x_start, y_start, f'{y_start:.3f}', ha='right')
    pyplot.text(x_end, y_end, f'{y_end:.3f}', ha='left')

    pyplot.plot([x_start, x_end], [y_start, y_end], color='blue',
                dashes=[4, 3], label='Тенденция')

    text_1 = f"ADF p-значение: {adf_p_value:.5f}"
    text_2 = f"KPSS p-значение: {kpss_p_value:.5f}"
    pyplot.text(0.845, 0.2, s=text_1 + '\n' + text_2,
                ha='left', va='bottom', transform=pyplot.gca().transAxes,
                bbox=dict(facecolor='blue', alpha=0.2))

    pyplot.xlabel('Дата наблюдения', fontsize=12)
    pyplot.ylabel('Температура (град. C)', fontsize=12)
    pyplot.legend(loc='lower right')
    pyplot.title(f"\"{timeseries.name}\": скользящее среднее и стандартное отклонение (окно = {window})")
    pyplot.grid(True)
    pyplot.show()

def test_stationarity(timeseries, plot=True):
    """Выполнение углублённой проверки временного ряда на стационарность"""
    print('\n'f"Выполнение расширенной проверки ВР \"{timeseries.name}\" на стационарность:")
    adf_test = adfuller(timeseries)
    adf_result = Series(adf_test[0:4],
                        index=['Тестовая статистика',
                                'p-значение',
                                'Использованные лаги',
                                'Кол-во использованных наблюдений']
                        )
    for key, value in adf_test[4].items():
        adf_result[f"Критическое значение {key}"] = value
    adf_result = round_series_floats(adf_result)

    kpss_test = kpss(timeseries)
    kpss_result = Series(kpss_test[0:3], index=['Тестовая статистика',
                                                'p-значение',
                                                'Использованные лаги'])
    for key, value in kpss_test[3].items():
        kpss_result[f"Критическое значение {key}"] = value
    kpss_result = round_series_floats(kpss_result)

    msg = "Результаты теста Дики-Фуллера (ADF-тест):"
    print('\n' + msg)
    print('-' * len(msg))
    for index, value in adf_result.items():
        sep = ' ' * (40 - len(index))
        print(f'{index}: {sep} {value}')

    msg = "Результаты теста Квятковского-Филлипса-Шмидта-Шина (KPSS-тест):"
    print('\n' + msg)
    print('-' * len(msg))
    for index, value in kpss_result.items():
        sep = ' ' * (40 - len(index))
        print(f'{index}: {sep} {value}')

    if plot:
        plot_movingAverage(timeseries=timeseries,
                           adf_p_value=adf_result['p-значение'],
                           kpss_p_value=kpss_result['p-значение'])

        
test_stationarity(timeseries=temperature)

Кодможет показаться сложными для восприятия, однако с его помощью мы выведем всю необходимую информацию, которая позволит сделать правильные выводы о стационарности временного ряда. Функция test_stationarity() выполнит вышеописанные статистические тесты и выведет их полные результаты в удобном для чтения виде. Благодаря использованию функции plot_movingAverage() будет построен график, который поможет наглядно убедиться в наличии положительной тенденции в данных; на графике будут изображены как оригинальный временной ряд в качестве фона, так и скользящее среднее со стандартным отклонением с окном, равным 8760. При этом мы обозначим на графике начальную и конечную точки скользящей средней, координаты которых определим по её данным, а также для удобства обозначим рядом с ними их соответствующие значения температуры. По этим точкам проведём отрезок, по которому будем судить о тенденции в данных. Дополнительно обозначим на графике p‑значения обоих тестов.

Согласитесь, не так уж и сложно.

Результаты выполнения функций приведены ниже.

Если не изображать дополнительный отрезок вместе с сеткой графика, то можно было бы предположить, что как тенденция, так и стандартное отклонение постоянны во времени (= стационарны). Однако, как мы можем видеть, это не так. По графику видна еле заметная тенденция увеличения температуры со временем.

Итак, с графиком разобрались: мы убедились, что в данных временного ряда с температурой имеется тенденция к её увеличению. Следом посмотрим на результаты статистических тестов. Для удобства объединим их вместе с анализом в общую таблицу.

Таким образом, по вышеизложенным результатам и их анализу следует, что временной ряд «Температура» является стационарным по расширенному тесту Дики‑Фуллера и нестационарным по тесту КПСС. С помощью дополнительно построенного графика с оригинальным временным рядом вместе со скользящей средней и стандартным отклонением по периоду, равному 8760, а также отрезку, проведённому по координатам конечных точек скользящего среднего, мы убедились, что временной ряд имеет тенденцию к увеличению температуры. В этой связи мы принимаем, что временной ряд является нестационарным как по сезонности, так и по тенденции и, следовательно, необходимо провести операцию взятия последовательной разности для достижения его стационарности.

Давайте выполним дифференцирование [ссылка] температуры на один порядок (period=1) и проведём тест на стационарность с помощью следующей функции, основа которой вам уже хорошо знакома:

def series_p_value_verification(series, periods=None):
    """Проверка временных рядов на стационарность по двум тестам"""
    msg = "Выполнение расширенного теста Дики-Фуллера:"
    print('\n' + msg)
    print('-' * len(msg))

    if periods:
        name = series.name
        series = series.diff(periods=periods).dropna()
        series.rename(name + ' diff=' + str(periods), inplace=True)

    adf_result = adfuller(series)
    if adf_result[1] < 0.05:
        h_text = colored('стационарен', 'grey', 'on_green')
        print(f"Ряд {series.name} {h_text}: p-value = {adf_result[1]:.5f}")
    else:
        h_text = colored('не стационарен', 'grey', 'on_red')
        print(f"Ряд {series.name} {h_text}: p-value = {adf_result[1]:.5f} > 0.05")

    msg = "Выполнение теста Квятковского-Филлипса-Шмидта-Шина:"
    print('\n' + msg)
    print('-' * len(msg))

    kpss_result = kpss(series)
    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")


series_p_value_verification(series=temperature, periods=1)

Результаты представлены ниже.

По результатам видно, что одного порядка разности оказалось достаточно для достижения стационарности и его дальнейшее дифференцирование не требуется. Таким образом, параметр d для целевого временного ряда принимаем равным 1. Напомню, что вручную выполнять дифференцирование эндогенного временного ряда перед обучением модели не нужно: модель ARIMA самостоятельно выполнит операцию взятия последовательной разности согласно установленному параметру d. При этом надо иметь в виду, что «на практике d обычно равно 0, 1 или максимум 2», как отмечают в своей книге Дж. Бокс и Г.Дженкинс.

Для определения параметров p и q воспользуемся следующей функцией, которая построит график стационарного временного ряда с разностью в один порядок, графики автокорреляционной функции (АКФ) и частной автокорреляционной функции (ЧАКФ) с временным запаздыванием, равным 120, а также АКФ с количеством лагов, равным количеству наблюдений временного ряда для лучшего представления информации о долгосрочной связи между уровнями (лагами) стационарного временного ряда. С помощью построенных графиков АКФ мы также проанализируем наши данные на сезонность.

def plot_target_analysis(series, periods=None, lags=None):
    """График с анализом характеристик целевого признака для определения параметров p, q"""
    if periods:
        series = series.diff(periods=periods).dropna()

    fig, axs = pyplot.subplots(2, 3)

    axs[0, 0].plot(series)
    axs[0, 0].set_title(f"Ряд с разностью {periods}-го порядка")
    plot_acf(series, ax=axs[0, 1], lags=lags)
    axs[0, 1].set_title(f"АКФ ({periods}-й порядок)")
    plot_pacf(series, ax=axs[0, 2], lags=lags)
    axs[0, 2].set_title(f"ЧАКФ ({periods}-й порядок)")

    for ax in axs[1, :]:
        ax.remove()

    axs[1, 0] = fig.add_subplot(212)

    plot_acf(series, ax=axs[1, 0], lags=arange(len(series)), marker='')
    axs[1, 0].set_title(f"График автокорреляции временного ряда {periods}-го порядка")

    pyplot.show()


plot_target_analysis(series=temperature, periods=1, lags=120)

Хмм… Дорогие друзья, а мы приехали!

Помимо годовой сезонности, которая представлена на нижнем графике АКФ в виде синусоиды, затухающей с увеличением количества лагов (далее мы выведем этот график отдельно для лучшего понимания), обнаружилась и сильная дневная сезонность, равная 24, которая отчётливо видна на верхнем графике АКФ. Сложность в том, что модель SARIMA напрямую не поддерживает более одной сезонности и считается, что её выполнение в этом случае будет не самой лучшей идеей. Отмечу, что для подобных случаев с несколькими сезонностями в данных были разработаны специальные модели, одна из которых TBATS. TBATS — это аббревиатура по заложенным в ней компонентам — Trigonometric, Box‑Cox transformation, ARMA errors, Trend and Seasonal components, что по‑русски можно перевести как «Тригонометрические функции — преобразование Бокса‑Кокса — ошибки АРСС — тенденция — сезонные компоненты». Эта модель заслуживает отдельного рассмотрения не только из‑за возможности учёта нескольких сезонных компонент, но и возможности применения преобразования Бокса‑Кокса к данным временного ряда для получения их нормального распределения [ссылка], о котором было упомянуто в конце первой главы нашего исследования (нормальное распределение измеряемых характеристик является важным моментом в статистике не только из‑за простоты анализа, но и лёгкости интерпретации результатов тестов). Однако и у неё есть свои недостатки, главный из которых — медленное выполнение вычислений.

Строгоговоря, для достижения наилучших результатов при наличии нескольких сезонностей во временном ряду хорошей практикой будет использование нескольких стратегий и подходов. Решение в лоб — разделить данные по годам и моделировать каждую сезонность отдельно, а затем суммировать результаты. Наши данные позволяют это сделать. Но мы оставим эту задачу последователям, а сами пойдём по сложному и чреватому на разочарование пути — и попытаемся выполнить модель SARIMA для задачи прогнозирования температуры с учётом двух сезонностей.

Но об этом, дорогие друзья, мы поговорим с вами в следующей, завершающей части второй главы нашего исследования.