Прогнозирование временных рядов в Python: от хаоса в данных к работающей модели SARIMA

Привет! Умение "заглядывать" в будущее — очень важный навык при работе с данными. Как изменятся продажи в следующем квартале? Какой спрос на товар ожидать перед праздниками? Какова будет посещаемость сайта через месяц?

Ответить на эти вопросы помогает анализ временных рядов.

Анализ временных рядов (Time series analysis, TSA) – это метод изучения характеристик целевой переменной во времени, где время является независимой переменной. А временной ряд — это, по сути, любая последовательность данных, измеренных через равные промежутки времени. TSA используется в различных областях, например, для прогнозирования погодных показателей (температуры, осадков, давления и т.д.), финансовых рынков (цен на акции, курсов валют, процентных ставок и т.д.), продаж (спроса, объемов производства, запасов и т.д.).

И хорошая новость в том, что все необходимые инструменты для этого уже есть в экосистеме Python.

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

Наша цель: построить на Python модель, которая проанализирует исторические данные о пассажиропотоке и сможет спрогнозировать его значения на будущее.

Мы пройдем все ключевые этапы:

• Первичный анализ: научимся загружать данные и интерпретировать их с помощью визуализации.
• Декомпозиция: разберем ряд на его базовые компоненты — тренд, сезонность и остатки.
• Подготовка данных: узнаем, что такое стационарность и почему это критически важно для большинства моделей.
• Построение модели: соберем все знания воедино и обучим модель ARIMA делать точные предсказания.

К концу этого руководства у тебя будет не только работающий код, но и ясное понимание принципов, стоящих за ним. Приступаем! 🚀

Что скрывает наш временной ряд?

Как и в любом data-проекте, наш первый шаг — познакомиться с данными вживую. Прежде чем строить сложные модели, нужно понять, с чем мы имеем дело, каковы особенности и структура нашего материала.

Для нашей учебной задачи мы возьмем хрестоматийный датасет AirPassengers. Это данные о ежемесячном количестве пассажиров международных авиалиний с 1949 по 1960 год.

Почему он так хорош для старта?

1. Чистый сигнал. В данных кристально четко прослеживаются и тренд (с годами люди летали все больше), и сезонность (каждое лето — пик отпусков). На таком примере очень легко объяснять ключевые концепции.
2. Простота. Датасет состоит всего из двух колонок: Month и #Passengers. Никаких лишних сущностей, которые могли бы отвлечь нас от основной задачи — изучения временного ряда.

Загрузка и первичное знакомство с данными

Перейдем к коду. Первым делом импортируем pandas и загрузим наш CSV-файл.

Сразу при загрузке мы сделаем несколько важных вещей:

• index_col='Month': Укажем, что колонка Month должна стать индексом нашего DataFrame. Временные ряды почти всегда удобнее анализировать, когда время является индексом.
• parse_dates=['Month']: Дадим pandas команду не просто прочитать даты как строки, а сразу преобразовать их в специальный datetime формат. Это откроет нам доступ к множеству удобных методов для работы со временем.

import pandas as pd

# Загружаем данные с прямой ссылки
url = 'https://raw.githubusercontent.com/obulygin/content/refs/heads/main/Airline%20Passengers/airline-passengers.csv'
df = pd.read_csv(url, index_col='Month', parse_dates=['Month'])

# Переименуем колонку для удобства, убрав спецсимволы
df = df.rename(columns={'#Passengers': 'Passengers'})

# Посмотрим на первые 5 строк
print(df.head())

Вывод будет таким:

Passengers
Month                 
1949-01-01         112
1949-02-01         118
1949-03-01         132
1949-04-01         129
1949-05-01         121

Отлично. Теперь проверим, все ли pandas понял правильно, и посмотрим на общую информацию о нашем датафрейме с помощью метода .info().

# Проверяем типы данных и наличие пропусков
df.info()

Результат:

<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 144 entries, 1949-01-01 to 1960-12-01
Data columns (total 1 columns):
 #   Column      Non-Null Count  Dtype
---  ------      --------------  -----
 0   Passengers  144 non-null    int64
dtypes: int64(1)
memory usage: 2.2 KB

Что мы видим?

1. DatetimeIndex: Индекс нашего датафрейма имеет правильный тип. pandas распознал его как даты.
2. 144 entries: У нас 144 записи, что соответствует 12 годам ежемесячных данных (12 * 12 = 144).
3. Non-Null Count: В колонке Passengers 144 непустых значения. Это значит, что в данных нет пропусков, и нам не придется тратить время на их обработку.

Данные загружены, очищены и готовы к первому визуальному осмотру.

Визуальный осмотр

Числа в таблице редко бывают красноречивыми. Чтобы по-настоящему "почувствовать" данные и уловить их характер, нужно построить график. Для этого воспользуемся вечной классикой — библиотекой matplotlib.

Поскольку наш индекс уже имеет тип datetime, pandas и matplotlib легко "договорятся" между собой, чтобы построить красивый и информативный график.

import matplotlib.pyplot as plt

# Задаем размер графика
plt.figure(figsize=(14, 7))

# Строим график
plt.plot(df.index, df['Passengers'])

# Добавляем заголовок и подписи осей
plt.title('Динамика пассажиропотока (1949-1960)')
plt.xlabel('Год')
plt.ylabel('Количество пассажиров')

# Включаем сетку для лучшей читаемости
plt.grid(True)

# Отображаем график
plt.show()

А вот теперь начинается самое интересное — анализ. Что этот график нам рассказывает?

1. Линия неуклонно ползет вверх. Это значит, что с течением времени общее количество авиапассажиров росло.
2. Внутри каждого года мы видим повторяющийся "узор": спад в начале года, подъем к лету и снова спад к зиме. Эти предсказуемые, циклические колебания, привязанные к определенному времени года (в нашем случае — к сезону отпусков), называются сезонностью.
3. Обрати внимание на важную деталь: летние "пики" в конце 50-х гораздо выше, чем пики в начале 50-х. А зимние "провалы" — глубже. Колебания не просто существуют, их размах (амплитуда) увеличивается вместе с ростом общего тренда.

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

Компоненты временного ряда

Мы только что на глаз определили ключевые компоненты нашего временного ряда. Теперь пришло время формализовать эти знания. Понимание внутренней структуры данных — это фундамент, на котором строятся точные прогнозы.

Тренд, сезонность и шум

Большинство классических временных рядов можно представить как комбинацию трех компонентов:

1. Тренд (Trend, T): Это то самое долгосрочное, основное направление движения ряда. В нашем случае — это общий рост числа пассажиров год от года. Тренд может быть восходящим, нисходящим или отсутствовать (горизонтальный).
2. Сезонность (Seasonality, S): Это те самые циклические, предсказуемые колебания вокруг тренда. Важно, что они происходят с фиксированной периодичностью: день, неделя, год. Наши летние пики и зимние спады — классическая годовая сезонность.
3. Остатки (Residuals или Шум, R): Это все, что осталось после вычитания тренда и сезонности. В идеальном мире это случайные, непредсказуемые колебания без какой-либо видимой структуры.

Важно отметить, что не все временные ряды содержат все эти компоненты.

Существует два основных способа, которыми эти компоненты могут "собираться" в итоговый ряд:

• Аддитивная модель: Y(t) = T(t) + S(t) + R(t)
  Подходит для рядов, где амплитуда сезонных колебаний постоянна и не зависит от общего уровня тренда.
• Мультипликативная модель: Y(t) = T(t) * S(t) * R(t)
  Подходит для рядов, где амплитуда сезонных колебаний растет или убывает вместе с трендом. Наше наблюдение из прошлой главы о растущем размахе колебаний — прямое указание на этот тип модели.

Понимание типа модели (аддитивная или мультипликативная) критически важно для правильного анализа. К счастью, в Python есть инструменты, которые помогут нам не только подтвердить наши догадки, но и наглядно разделить ряд на эти три составляющие.

Декомпозиция временного ряда

Пора вернуться к коду. Процесс разложения временного ряда на его компоненты называется декомпозицией. В библиотеке statsmodels, которая является стандартом де-факто для статистического анализа в Python, есть для этого удобная функция seasonal_decompose.

Основные аргументы этой функции:

• series: временной ряд, который необходимо декомпозировать (в нашем случае df['Passengers']).
• model: вид декомпозиции, может быть 'additive' (аддитивная) или 'multiplicative' (мультипликативная).
• period: период сезонности, например, 12 для ежемесячных данных или 4 для ежеквартальных.

Функция возвращает объект seasonal_decompose, который содержит следующие атрибуты:

• trend: тренд временного ряда.
• seasonal: сезонная компонента.
• resid: нерегулярная компонента / остатки.
• observed: исходный временной ряд.

from statsmodels.tsa.seasonal import seasonal_decompose

# Выполняем декомпозицию
# period=12, так как у нас годовая сезонность с 12 месяцами
result = seasonal_decompose(df['Passengers'], model='multiplicative', period=12)

# Визуализируем результат
fig = result.plot()
fig.set_size_inches(14, 10)
plt.show()

Давай разберем, что мы получили на этих четырех графиках:

1. Observed: Это просто наш исходный ряд, с которым мы уже знакомы.
2. Trend: А вот и наш тренд, выделенный в чистом виде! Функция сгладила сезонные колебания, и теперь мы видим только основное направление движения ряда.
3. Seasonal: Здесь показан сам сезонный паттерн. Обрати внимание, что по оси Y значения колеблются вокруг 1.0. Это и есть наш мультипликативный сезонный коэффициент. Например, значение 1.2 в летний месяц означает, что пассажиропоток в этом месяце в среднем на 20% выше, чем годовой тренд. Значение 0.85 в зимний месяц — на 15% ниже.
4. Resid: Те самые остатки, или шум. В идеале, здесь не должно быть никакой структуры — просто случайные точки. Наш график выглядит достаточно хаотичным, что является хорошим знаком.

Мы успешно извлекли из ряда его составные части. Теперь, когда мы понимаем его анатомию, мы готовы перейти к следующему, возможно, самому важному и нетривиальному концепту в анализе временных рядов — к стационарности.

Стационарность временных рядов

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

Что такое стационарный ряд?

Если говорить просто, временной ряд называется стационарным, если его статистические свойства — среднее, дисперсия и автоковариация — не изменяются во времени.

Давай расшифруем на пальцах:

• Постоянное среднее: У ряда нет тренда. Он колеблется вокруг одного и того же уровня.
• Постоянная дисперсия: "Размах" колебаний ряда примерно одинаков на всем его протяжении. Он не становится более "диким" или более "спокойным" со временем.
• Постоянная автоковариация: Связь между значениями ряда Y(t) и Y(t+k) зависит только от лага k (расстояния между точками), а не от того, в какой части ряда мы находимся — в начале или в конце.

Наш ряд AirPassengers — это эталонный пример нестационарного ряда.

1. У него есть ярко выраженный восходящий тренд (среднее растет).
2. У него есть растущая дисперсия (размах колебаний увеличивается).

Почему это так важно? Представь, что ты пытаешься построить модель на данных о росте ребенка. Прогнозировать его рост в 20 лет на основе данных от 1 до 5 лет — гиблое дело, потому что "правила игры" (скорость роста) постоянно меняются. Модели временных рядов работают так же: им нужна стабильность. Они могут выучить закономерности только в "стабильной" среде, где правила игры не меняются.

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

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

Тесты Дики-Фуллера (ADF) и KPSS

Полагаться только на визуальный анализ рискованно — иногда тренд или сезонность могут быть не так очевидны. Чтобы принять объективное решение, используют статистические тесты. Самый популярный из них — расширенный тест Дики-Фуллера (Augmented Dickey-Fuller Test, ADF).

Суть теста проста. Он проверяет одну-единственную гипотезу:

• Нулевая гипотеза (H₀): Ряд нестационарен (в нем присутствует "единичный корень").
• Альтернативная гипотеза (H₁): Ряд стационарен.

Наша задача — попытаться опровергнуть нулевую гипотезу. Мы делаем это с помощью p-value, которое возвращает тест:

• Если p-value > 0.05: У нас нет достаточных оснований, чтобы отвергнуть H₀. Мы принимаем, что ряд нестационарен.
• Если p-value <= 0.05: Мы отвергаем H₀ и с уверенностью можем считать наш ряд стационарным.

Давай применим тест к нашему исходному ряду. Для этого в statsmodels есть готовая функция adfuller.

from statsmodels.tsa.stattools import adfuller

# Проводим ADF-тест
result_adf = adfuller(df['Passengers'])

# Выводим результат
print('ADF Statistic:', result_adf[0])
print('p-value:', result_adf[1])

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

ADF Statistic: 0.8153688792060482
p-value: 0.991880243437641

Наш p-value равен 0.99. Это значительно больше порогового значения 0.05.

Вывод: Тест ADF с уверенностью говорит нам: "Ребята, у вас нет никаких оснований считать этот ряд стационарным". Наша визуальная догадка полностью подтвердилась.

Для полной уверенности часто используют еще один тест — KPSS (Kwiatkowski-Phillips-Schmidt-Shin). KPSS-тест основывается на анализе суммы квадратов отклонений накопленных частичных сумм от линейного тренда. Он работает «наоборот»:

• Нулевая гипотеза (H₀): Ряд стационарен.
• Альтернативная гипотеза (H₁): Ряд нестационарен.

Здесь мы, наоборот, хотим получить p-value меньше 0.05, чтобы отвергнуть гипотезу о стационарности."

Применим и его.

from statsmodels.tsa.stattools import kpss

# Проводим KPSS-тест
result_kpss = kpss(df['Passengers'], regression='c')

# Выводим результат
print('\nKPSS Test:')
print(f'KPSS Statistic: {result_kpss[0]}')
print(f'p-value: {result_kpss[1]}')

Результат:

KPSS Test:
KPSS Statistic: 1.6513122354165206
p-value: 0.01

InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is smaller than the p-value returned.

p-value здесь равен 0.01 (и даже меньше), что ниже порога 0.05. Это значит, мы отвергаем нулевую гипотезу KPSS о том, что ряд стационарен.

Финальный вердикт: Оба теста сошлись во мнении. ADF не смог отвергнуть гипотезу о нестационарности, а KPSS уверенно отверг гипотезу о стационарности. Диагноз окончательный и обжалованию не подлежит: наш ряд нестационарен.

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

Преобразование нестационарных данных в стационарные

Итак, у нашего ряда две проблемы:

1. Нестабильная дисперсия (размах колебаний растет со временем).
2. Явный тренд (среднее значение ряда растет со временем).

Будем решать их последовательно.

Шаг 1: Стабилизируем дисперсию с помощью логарифмирования

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

import numpy as np

# Применяем натуральный логарифм
df['Passengers_log'] = np.log(df['Passengers'])

# Визуализируем результат
plt.figure(figsize=(14, 7))
plt.plot(df.index, df['Passengers_log'])
plt.title('Пассажиропоток после логарифмирования')
plt.xlabel('Год')
plt.ylabel('Log(Пассажиров)')
plt.grid(True)
plt.show()

Посмотрите, как изменился график! Тренд никуда не делся, но "пики" и "спады" теперь имеют примерно одинаковую амплитуду. Мы победили растущую дисперсию! Но ряд все еще нестационарен из-за тренда, давай убедимся в этом с помощью тестов.

# Проверяем логарифмированный ряд
result_adf = adfuller(df['Passengers_log'])
print(f'ADF p-value: {result_adf[1]}')
# ADF p-value: 0.42236677477039125

result_kpss = kpss(df['Passengers_log'], regression='c')
print(f'KPSS p-value: {result_kpss[1]}')
# KPSS p-value: 0.01

ADF все еще не может отвергнуть гипотезу о нестационарности (p > 0.05), а KPSS уверенно отвергает гипотезу о стационарности (p < 0.05). Вердикт: проблему дисперсии мы решили, но тренд все еще делает ряд нестационарным.

Шаг 2: Убираем тренд с помощью дифференцирования

Теперь наша главная цель — избавиться от тренда. Самый мощный инструмент для этого — дифференцирование, то есть переход от анализа абсолютных значений к анализу их разностей. Попробуем несколько стратегий.

Стратегия А: Простое дифференцирование

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

diff(Y_t) = Y_t - Y_{t-1}

В pandas для этого есть удобный метод .diff(). Применим его к нашему логарифмированному ряду.

# Берем первую разность от логарифмированного ряда
log_diff = df['Passengers_log'].diff().dropna()

# Почему .dropna()?
# При вычислении разности для самого первого элемента (1949-01-01)
# нет предыдущего значения, поэтому результат — NaN (Not a Number).
# Статистические тесты не работают с пропусками, поэтому мы обязаны удалить эту строку.

Проверим, что получилось:

result_adf = adfuller(log_diff)
print(f'ADF p-value: {result_adf[1]}')
# ADF p-value: 0.07112054815085875

result_kpss = kpss(log_diff, regression='c')
print(f'KPSS p-value: {result_kpss[1]}')
# KPSS p-value: 0.1

Попробуем другой подход.

Стратегия Б: Сезонное дифференцирование

Логичный следующий шаг — атаковать сезонность напрямую. Для этого применим сезонное дифференцирование: будем вычитать значение не за прошлый месяц, а за аналогичный месяц прошлого года (т.е. с лагом 12).

# Берем сезонную разность от логарифмированного ряда
log_seasonal_diff = df['Passengers_log'].diff(12).dropna()

# Теперь мы теряем первые 12 месяцев (весь 1949 год),
# так как для них не существует данных за 12 месяцев до этого.
# Их также необходимо удалить.

И снова проверка:

result_adf = adfuller(log_seasonal_diff)
print(f'ADF p-value: {result_adf[1]}')
# ADF p-value: 0.07239567181769484

result_kpss = kpss(log_seasonal_diff, regression='c')
print(f'KPSS p-value: {result_kpss[1]}')
# KPSS p-value: 0.1

Опять почти, но не совсем! Мы убрали сезонные пики, но, видимо, какой-то остаточный тренд все еще присутствует в данных.

Стратегия В: Комбинированный удар!

Раз ни один метод не сработал поодиночке, объединим их силу! Применим и сезонное, и простое дифференцирование последовательно. Это стандартный прием для рядов с сильным трендом и сильной сезонностью.

# Сначала берем сезонную разность, затем от результата — простую
final_diff = df['Passengers_log'].diff(12).diff().dropna()

Давайте посмотрим на график этого финального ряда и проверим его тестами.

# Визуализация
plt.figure(figsize=(14, 7))
plt.plot(final_diff.index, final_diff)
plt.title('Ряд после всех преобразований')
plt.xlabel('Год')
plt.ylabel('Двойная разность логарифмов')
plt.grid(True)
plt.show()

# Финальная проверка
result_adf = adfuller(final_diff)
print(f'ADF p-value: {result_adf[1]}')
# ADF p-value: 0.001215446476048047

result_kpss = kpss(final_diff, regression='c')
print(f'KPSS p-value: {result_kpss[1]}')
# KPSS p-value: 0.1

Бинго!

• ADF: p-value значительно меньше 0.05. Мы уверенно отвергаем гипотезу о нестационарности.
• KPSS: p-value все еще 0.1, что означает, что мы не можем отвергнуть его нулевую гипотезу о стационарности.

Оба теста наконец-то сошлись: наш ряд стационарен.

Мы прошли самый сложный и самый важный этап подготовки данных. Теперь, имея на руках "послушный" стационарный ряд, мы наконец-то готовы строить нашу машину времени — модель ARIMA.

Строим модель SARIMA

Мы вплотную подошли к финалу. Все предыдущие главы были тщательной подготовкой к этому моменту. Теперь мы будем собирать все наши знания в единую модель. Что за модель?

Autoregressive Integrated Moving Average (ARIMA) – модель авторегрессионной интегрированной скользящей средней.

• AR ==> Использует прошлые значения для прогнозирования будущего.
• MA ==> Использует прошлые значения ошибок в данном ряду для прогнозирования будущего.
• I ==> Использует разности наблюдений, чтобы сделать данные стационарными.
• AR + I + MA = ARIMA

Модели ARIMA определяются тремя параметрами:

• p (AR - Autoregressive): Порядок авторегрессии. Сколько прошлых значений нашего ряда мы используем для предсказания текущего.
• d (I - Integrated): Порядок простого дифференцирования. Сколько раз мы применяли .diff() для удаления тренда.
• q (MA - Moving Average): Порядок скользящего среднего. Сколько прошлых ошибок прогноза мы используем для предсказания текущего значения.

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

Расшифровываем SARIMA

Это две ARIMA-модели в одной: одна для несезонной части, другая — для сезонной. То есть добавляются сезонные параметры (P, D, Q)m:

• P (Seasonal AR): Порядок сезонной авторегрессии. То же, что и p, но для сезонных лагов (например, смотрим на значение 12 месяцев назад, 24 месяца назад и т.д.).
• D (Seasonal I): Порядок сезонного дифференцирования. Сколько раз мы применяли .diff(m) для удаления сезонности.
• Q (Seasonal MA): Порядок сезонного скользящего среднего. То же, что и q, но для ошибок на сезонных лагах.
• m: Период сезонности. В нашем случае это 12, так как у нас годовая сезонность на ежемесячных данных.

Чтобы найти эти недостающие параметры, нам понадобятся еще два инструмента — автокорреляционная и частично автокорреляционная функции.

Ищем p, q, P и Q с помощью ACF и PACF

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

ACF (Autocorrelation Function) — функция автокорреляции. Измеряет степень сходства между значением временного ряда и его прошлыми значениями с разными временными лагами. Другими словами, ACF измеряет корреляцию между временными отклонениями ряда и его отставаниями на различные лаги.

Корреляция, близкая к 1 или -1, указывает на сильную зависимость между значениями на различных временных точках, в то время как корреляция, близкая к 0, указывает на отсутствие зависимости.

На практике ACF позволяет определить наличие сезонных или циклических паттернов во временном ряду.

PACF (Partial Autocorrelation Function) — функция частичной автокорреляции. Показывает "чистую", прямую корреляцию между Y_t и Y_{t-k}, убирая влияние всех промежуточных точек.

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

Есть классическое правило (хоть и не всегда идеально работающее), которое помогает определить параметры по графикам ACF и PACF:

• Параметры q и Q (MA) определяются по ACF. Ищем последний значимый лаг перед тем, как график резко обрывается и уходит в незначимую зону (синий доверительный интервал).
• Параметры p и P (AR) определяются по PACF. Логика та же — ищем последний значимый лаг перед обрывом.

Давайте построим эти графики для нашего финального, стационарного ряда (final_diff).

import statsmodels.api as sm

# Наш стационарный ряд из прошлой главы
final_diff = df['Passengers_log'].diff(12).diff().dropna()

# Строим графики
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 10))

# График ACF
sm.graphics.tsa.plot_acf(final_diff, lags=40, ax=ax1)
ax1.set_title('Автокорреляционная функция (ACF)')

# График PACF
sm.graphics.tsa.plot_pacf(final_diff, lags=40, ax=ax2)
ax2.set_title('Частичная автокорреляционная функция (PACF)')

plt.tight_layout()
plt.show()

Давайте внимательно их изучим:

1. Несезонные параметры (p, q):

• ACF: Видим один значимый пик на лаге 1, после чего корреляция резко падает. Это сильный кандидат на q = 1.
• PACF: Также видим один значимый пик на лаге 1. Это кандидат на p = 1.

2. Сезонные параметры (P, Q):

• ACF: Есть очень заметный значимый пик на лаге 12. Это наш главный сезонный лаг. Это сильный кандидат на Q = 1.
• PACF: Наблюдаем ту же картину — значимый пик на лаге 12 (и еще несколько вокруг него, 11 и 13, что типично для сезонных эффектов). Это кандидат на P = 1.

Выбор параметров по графикам ACF/PACF — это отчасти искусство, отчасти наука. Наша гипотеза выглядит очень хорошо, но окончательный вердикт вынесет сама модель. Теперь у нас есть все необходимое, чтобы наконец-то ее обучить.

Обучаем, проверяем и уточняем модель

Теперь у нас есть все, чтобы собрать и обучить нашу модель. Мы будем использовать класс SARIMAX из statsmodels. Обрати внимание на X в конце — он означает, что модель может работать и с экзогенными (внешними) переменными, но мы эту возможность использовать не будем.

Мы начнем с гипотезы, которую получили после анализа графиков ACF/PACF: SARIMAX(1, 1, 1)(1, 1, 1, 12).

Процесс состоит из трех шагов:

1. Создать экземпляр модели, передав ей наши данные и подобранные параметры. Важно: модель обучается на исходных, недифференцированных данных (df['Passengers_log']), а параметры d, D и т.д. сообщают ей, какие преобразования нужно выполнить "под капотом".
2. Обучить модель с помощью метода .fit().
3. Сделать прогноз с помощью метода .predict() или .forecast().

Давай напишем код.

# Создаем и обучаем модель SARIMA
# Обрати внимание, что мы используем логарифмированные данные
model = sm.tsa.SARIMAX(
    df['Passengers_log'],
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 12)
)

results = model.fit()

# Выведем сводную информацию о модели
print(results.summary())

После выполнения этого кода ты увидишь большую сводную таблицу с результатами.

                                     SARIMAX Results                                      
==========================================================================================
Dep. Variable:                     Passengers_log   No. Observations:                  131
Model:             SARIMAX(1, 1, 1)x(1, 1, 1, 12)   Log Likelihood                 218.619
Date:                            Tue, 01 Jul 2025   AIC                           -427.239
Time:                                    18:30:03   BIC                           -413.385
Sample:                                02-01-1950   HQIC                          -421.614
                                     - 12-01-1960                                         
Covariance Type:                              opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.2673      0.204      1.311      0.190      -0.132       0.667
ma.L1         -0.6519      0.175     -3.718      0.000      -0.996      -0.308
ar.S.L12      -0.0680      0.202     -0.337      0.736      -0.464       0.328
ma.S.L12      -0.5257      0.225     -2.333      0.020      -0.967      -0.084
sigma2         0.0014      0.000      8.022      0.000       0.001       0.002
===================================================================================
Ljung-Box (L1) (Q):                   0.12   Jarque-Bera (JB):                 2.89
Prob(Q):                              0.73   Prob(JB):                         0.24
Heteroskedasticity (H):               0.51   Skew:                            -0.04
Prob(H) (two-sided):                  0.04   Kurtosis:                         3.76
===================================================================================

А теперь — самый важный этап анализа. Смотрим на таблицу и столбец P>|z| (p-value коэффициентов).

• ar.L1 (p) -> 0.190 ( > 0.05) — незначим!
• ma.L1 (q) -> 0.000 ( < 0.05) — значим.
• ar.S.L12 (P) -> 0.736 ( > 0.05) — незначим!
• ma.S.L12 (Q) -> 0.020 ( < 0.05) — значим.

Модель говорит нам, что компоненты авторегрессии (AR) не вносят значимого вклада. В такой ситуации правильный шаг — упростить модель, оставив только значимые параметры.

Упрощаем модель: SARIMA(0, 1, 1)(0, 1, 1, 12)

Теперь мы обучим новую, более простую модель и посмотрим на ее результаты.

# Создаем и обучаем ФИНАЛЬНУЮ, УТОЧНЕННУЮ модель
model_final = sm.tsa.SARIMAX(
    df['Passengers_log'],
    order=(0, 1, 1),
    seasonal_order=(0, 1, 1, 12)
)
results_final = model_final.fit()

# Смотрим на новую сводку
print(results_final.summary())

                                     SARIMAX Results                                      
==========================================================================================
Dep. Variable:                     Passengers_log   No. Observations:                  131
Model:             SARIMAX(0, 1, 1)x(0, 1, 1, 12)   Log Likelihood                 218.067
Date:                            Tue, 01 Jul 2025   AIC                           -430.133
Time:                                    18:35:58   BIC                           -421.821
Sample:                                02-01-1950   HQIC                          -426.758
                                     - 12-01-1960                                         
Covariance Type:                              opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ma.L1         -0.4051      0.072     -5.609      0.000      -0.547      -0.264
ma.S.L12      -0.5552      0.114     -4.882      0.000      -0.778      -0.332
sigma2         0.0014      0.000      8.307      0.000       0.001       0.002
===================================================================================
Ljung-Box (L1) (Q):                   0.02   Jarque-Bera (JB):                 1.28
Prob(Q):                              0.88   Prob(JB):                         0.53
Heteroskedasticity (H):               0.49   Skew:                             0.04
Prob(H) (two-sided):                  0.03   Kurtosis:                         3.51
===================================================================================

Теперь таблица с коэффициентами выглядит гораздо лучше — все p-value в ней ниже 0.05, что говорит о значимости всех оставшихся параметров. Мы получили более простую и более "здоровую" модель.

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

Оцениваем полученную модель

Как понять, что наша модель действительно хороша? Один из важных критериев: ее остатки (residuals), то есть ошибки прогноза на исторических данных, должны быть максимально похожи на случайный "белый шум". Если в остатках есть какая-то структура (тренд, сезонность, корреляция), это значит, что модель не смогла извлечь всю полезную информацию из данных.

В statsmodels есть встроенная функция plot_diagnostics, которая создает 4 графика для проверки модели.

# Строим диагностические графики
results_final.plot_diagnostics(figsize=(15, 12))
plt.show()

Давай разберем, что мы видим и что это значит:

1. Вверху слева (Standardized residual): График стандартизированных остатков. Он должен выглядеть как хаотичные колебания вокруг нуля, без каких-либо явных паттернов. Наш график именно таков — это хороший знак.
2. Вверху справа (Histogram plus estimated density): Гистограмма распределения остатков. Зеленая кривая (KDE) показывает это распределение, а оранжевая (N(0,1)) — стандартное нормальное распределение. Они должны быть максимально похожи. У нас они почти совпадают, что подтверждает гипотезу о нормальности ошибок.
3. Внизу слева (Normal Q-Q): Еще одна, более строгая проверка на нормальность. Синие точки должны лежать на красной диагональной линии. Небольшие отклонения по краям допустимы. У нас почти идеальное совпадение.
4. Внизу справа (Correlogram): Это ACF-график для остатков. Он показывает, коррелируют ли ошибки модели друг с другом. В идеале все столбцы (кроме нулевого лага) должны находиться внутри синего доверительного интервала. Это означает отсутствие значимой автокорреляции. Наш график полностью удовлетворяет этому требованию.

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

Делаем прогноз

Модель обучена, и ее параметры выглядят адекватно. Теперь сделаем прогноз на будущее! Давай предскажем значения на следующие 3 года (36 месяцев).

# Делаем прогноз на 36 месяцев вперед
forecast_steps = 36
forecast = results_final.get_forecast(steps=forecast_steps)

# Получаем предсказанные значения и доверительные интервалы
forecast_mean = forecast.predicted_mean
confidence_intervals = forecast.conf_int()

Прогноз получен! Но есть один нюанс: модель работала с логарифмированными данными. Чтобы вернуть прогноз к исходным масштабам (количеству пассажиров), нам нужно применить обратное преобразование — экспоненту (np.exp).

# Возвращаем прогноз к исходному масштабу
forecast_mean_exp = np.exp(forecast_mean)
confidence_intervals_exp = np.exp(confidence_intervals)

Теперь у нас есть все для финальной визуализации. Построим график, на котором будут и наши исторические данные, и прогноз на будущее.

# Визуализируем результат
plt.figure(figsize=(14, 8))

# Исходные данные
plt.plot(df.index, df['Passengers'], label='Исторические данные')

# Прогноз
plt.plot(forecast_mean_exp.index, forecast_mean_exp, color='red', label='Прогноз SARIMA')

# Доверительный интервал
plt.fill_between(
    confidence_intervals_exp.index,
    confidence_intervals_exp.iloc[:, 0],
    confidence_intervals_exp.iloc[:, 1],
    color='pink',
    alpha=0.7,
    label='Доверительный интервал 95%'
)

plt.title('Прогноз пассажиропотока на 3 года')
plt.xlabel('Год')
plt.ylabel('Количество пассажиров')
plt.legend(loc='upper left')
plt.grid(True)
plt.show()

Посмотри на этот график! Наша модель не просто продолжила общий тренд на рост, но и идеально уловила сезонный паттерн с его летними пиками и зимними спадами.

Заключение: куда двигаться дальше?

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

Давай кратко вспомним наш путь:

1. Разведка: Мы загрузили данные и с помощью простого графика обнаружили в них тренд и растущую сезонность.
2. Анатомия: Мы разложили ряд на компоненты, чтобы лучше понять его структуру.
3. Стабилизация: С помощью логарифмирования и двойного дифференцирования мы превратили нестационарный ряд в стационарный — пригодный для анализа.
4. Моделирование: Мы подобрали параметры для модели SARIMA, обучили ее, проанализировали результаты, улучшили модель и, наконец, построили точный прогноз.

Ты не просто скопировал код. Ты научился думать как аналитик: ставить гипотезы, проверять их с помощью тестов, интерпретировать неоднозначные результаты и принимать решения на основе данных. Это и есть самый ценный навык.

Наша модель — это только начало

Модель SARIMA, которую мы построили, — это мощный, но классический инструмент. Мир анализа временных рядов огромен, и если тебе понравился этот процесс, вот несколько направлений для дальнейшего развития:

• Автоматический подбор параметров: Вместо ручного анализа ACF/PACF можно использовать инструменты вроде auto_arima из библиотеки pmdarima, которые делают это автоматически.
• Модели на основе машинного обучения: Можно решать задачи прогнозирования с помощью градиентного бустинга, создавая лаговые признаки вручную.
• Нейросетевые подходы: Для очень сложных рядов часто используют рекуррентные нейронные сети (LSTM, GRU) или современные архитектуры типа N-BEATS и Transformer.
• Модель Prophet отлично справляется с рядами, имеющими несколько видов сезонности и пропущенные данные.

Анализ временных рядов — это не просто технический навык, это способ мышления. И я надеюсь, эта статья стала для тебя хорошей отправной точкой в этом.