Machine Learning Guru

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

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

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

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

Строгая стационарность (strict stationarity) — совместное распределение вероятностей любого набора наблюдений не зависит от времени. Другими словами, если сдвинуть временной ряд вперед или назад во времени, его вероятностные свойства останутся неизменными.

Слабая стационарность (weak stationarity) или стационарность в широком смысле (covariance stationarity) — среднее значение, дисперсия и автоковариация ряда не зависят от времени.

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

• Постоянное среднее значение: E[X_t] = μ для всех t;
• Постоянная дисперсия: Var(X_t) = σ² для всех t;
• Ковариация между X_t и X_{t+h} зависит только от временного лага h, а не от времени t: Cov(X_t, X_{t+h}) = γ(h) для всех t и h.

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

Почему стационарность так важна?

Стационарность имеет фундаментальное значение для анализа временных рядов по нескольким причинам:

1. Применимость статистических методов: Большинство классических методов анализа временных рядов, включая авторегрессионные модели и модели скользящего среднего, предполагают стационарность исходных данных.
2. Прогнозируемость: Стационарные ряды легче прогнозировать, поскольку их статистические свойства остаются постоянными с течением времени.
3. Интерпретируемость результатов: Результаты анализа стационарных рядов обычно более надежны и легче интерпретируются.
4. Предотвращение ложных корреляций: Нестационарность может привести к ложным корреляциям между переменными, что известно как проблема «ложной регрессии» (spurious regression).

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

Типы нестационарности

Нестационарность временных рядов может проявляться в различных формах. Понимание конкретного типа нестационарности критически важно для выбора правильного метода преобразования данных.

Тренд-нестационарность

Тренд-нестационарность возникает, когда временной ряд имеет детерминированную тенденцию роста или снижения. Математически это можно представить как:

X_t = f(t) + ε_t

где f(t) — детерминированная функция времени (тренд), а ε_t — стационарный случайный процесс.

Тренд может быть:

• Линейным: f(t) = α + βt;
• Полиномиальным: f(t) = α + β₁t + β₂t² + … + βₙtⁿ;
• Экспоненциальным: f(t) = αeᵝᵗ;
• Логарифмическим: f(t) = α + βln(t).

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

Разностно-стационарность (unit root non-stationarity)

Разностно-стационарные ряды характеризуются наличием единичного корня в характеристическом уравнении соответствующего авторегрессионного процесса. Простейшим примером является случайное блуждание:

X_t = X_{t-1} + ε_t

где ε_t — белый шум.

Такие ряды можно сделать стационарными путем взятия разностей. Если ряд становится стационарным после взятия d-й разности, говорят, что он интегрирован порядка d, что обозначается как I(d).

Финансовые временные ряды, такие как цены акций или валютные курсы, часто являются разностно-стационарными.

Сезонная нестационарность

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

• Годовым циклом (например, продажи зимней одежды выше в зимние месяцы);
• Недельным циклом (например, трафик веб-сайтов обычно ниже в выходные дни);
• Суточным циклом (например, потребление электроэнергии выше в дневное время).

Математически сезонный ряд можно представить как:

X_t = s_t + ε_t

где s_t — сезонная компонента с периодом S (s_t = s_{t+S}), а ε_t — стационарный случайный процесс.

Гетероскедастичность

Гетероскедастичность означает непостоянство дисперсии временного ряда. Даже если среднее значение ряда остается постоянным, изменение дисперсии со временем нарушает условие стационарности.

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

Структурные изменения

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

• Изменениями в экономической политике;
• Технологическими прорывами;
• Природными катастрофами;
• Геополитическими событиями.

Математически ряд со структурным изменением в точке τ можно представить как:

X_t = {
f₁(t) + ε_t, если t < τ
f₂(t) + ε_t, если t ≥ τ
}

где f₁(t) и f₂(t) — различные функции времени, а ε_t — стационарный случайный процесс.

Понимание типа нестационарности критически важно для выбора подходящего метода преобразования данных. В следующем разделе я рассмотрю способы проверки временных рядов на стационарность.

Методы проверки стационарности

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

Визуальный анализ

Визуальный анализ — это первый шаг при исследовании временного ряда. Он позволяет быстро выявить очевидные признаки нестационарности.

График временного ряда

Простой график временного ряда может многое рассказать о его свойствах:

• Наличие тренда (восходящего или нисходящего);
• Сезонные колебания;
• Изменение дисперсии со временем;
• Структурные сдвиги.

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

Приведу пример визуализации стационарного и нестационарного рядов на Python:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller

# Генерация данных
np.random.seed(42)
# Стационарный ряд - AR(1) процесс
n = 1000
phi = 0.7
stationary = np.zeros(n)
stationary[0] = 0
for t in range(1, n):
    stationary[t] = phi * stationary[t-1] + np.random.normal(0, 1)

# Нестационарный ряд - случайное блуждание
nonstationary = np.zeros(n)
nonstationary[0] = 0
for t in range(1, n):
    nonstationary[t] = nonstationary[t-1] + np.random.normal(0, 1)

# Визуализация
plt.figure(figsize=(14, 7))

plt.subplot(2, 1, 1)
plt.plot(stationary)
plt.title('Стационарный ряд (AR(1) процесс)')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(nonstationary)
plt.title('Нестационарный ряд (случайное блуждание)')
plt.grid(True)

plt.tight_layout()
plt.show()

Рис. 1: Стационарный ряд и нестационарный ряд

Автокорреляционная функция (ACF)

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

# Построение ACF для стационарного и нестационарного рядов
plt.figure(figsize=(14, 7))

plt.subplot(2, 1, 1)
sm.graphics.tsa.plot_acf(stationary, lags=50, alpha=0.05, title='ACF для стационарного ряда')

plt.subplot(2, 1, 2)
sm.graphics.tsa.plot_acf(nonstationary, lags=50, alpha=0.05, title='ACF для нестационарного ряда')

plt.tight_layout()
plt.show()

Рис. 2: График автокорреляции для стационарного временного ряда

Рис. 3: График автокорреляции для нестационарного временного ряда

Частная автокорреляционная функция (PACF)

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

# Построение PACF для стационарного и нестационарного рядов
plt.figure(figsize=(14, 7))

plt.subplot(2, 1, 1)
sm.graphics.tsa.plot_pacf(stationary, lags=50, alpha=0.05, title='PACF для стационарного ряда')

plt.subplot(2, 1, 2)
sm.graphics.tsa.plot_pacf(nonstationary, lags=50, alpha=0.05, title='PACF для нестационарного ряда')

plt.tight_layout()
plt.show()

Рис. 4: Частная автокорреляция стационарного временного ряда

Рис. 5: Частная автокорреляция нестационарного временного ряда

Статистические тесты на стационарность

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

Расширенный тест Дикки-Фуллера (ADF)

ADF-тест проверяет нулевую гипотезу о наличии единичного корня (нестационарности) против альтернативной гипотезы о стационарности ряда. Основное уравнение теста:

Δy_t = α + βt + γy_{t-1} + δ₁Δy_{t-1} + … + δₚΔy_{t-p} + ε_t

где Δy_t = y_t — y_{t-1}.

Если γ значимо меньше нуля, нулевая гипотеза отвергается, и ряд считается стационарным.

# Проведение ADF-теста
def perform_adf_test(series, title):
    result = adfuller(series)
    print(f"ADF-тест для {title}:")
    print(f"ADF статистика: {result[0]:.4f}")
    print(f"p-значение: {result[1]:.4f}")
    print(f"Критические значения:")
    for key, value in result[4].items():
        print(f"  {key}: {value:.4f}")
    print("Вывод:", end=" ")
    if result[1] <= 0.05:
        print("Временной ряд стационарен (Отвергаем H0)")
    else:
        print("Временной ряд нестационарен (Не отвергаем H0)")
    print()

perform_adf_test(stationary, "стационарного ряда")
perform_adf_test(nonstationary, "нестационарного ряда")

ADF-тест для стационарного ряда:
ADF статистика: -13.3900
p-значение: 0.0000
Критические значения:
1%: -3.4369
5%: -2.8644
10%: -2.5683
Вывод: Временной ряд стационарен (Отвергаем H0)

ADF-тест для нестационарного ряда:
ADF статистика: -1.7852
p-значение: 0.3878
Критические значения:
1%: -3.4369
5%: -2.8644
10%: -2.5683
Вывод: Временной ряд нестационарен (Не отвергаем H0)

Тест Квятковского-Филлипса-Шмидта-Шина (KPSS)

KPSS-тест, в отличие от ADF-теста, проверяет нулевую гипотезу о стационарности против альтернативной гипотезы о нестационарности. Это делает его полезным дополнением к ADF-тесту, поскольку они проверяют противоположные гипотезы.

from statsmodels.tsa.stattools import kpss

def perform_kpss_test(series, title):
    result = kpss(series)
    print(f"KPSS-тест для {title}:")
    print(f"KPSS статистика: {result[0]:.4f}")
    print(f"p-значение: {result[1]:.4f}")
    print(f"Критические значения:")
    for key, value in result[3].items():
        print(f"  {key}: {value:.4f}")
    print("Вывод:", end=" ")
    if result[1] <= 0.05:
        print("Временной ряд нестационарен (Отвергаем H0)")
    else:
        print("Временной ряд стационарен (Не отвергаем H0)")
    print()

perform_kpss_test(stationary, "стационарного ряда")
perform_kpss_test(nonstationary, "нестационарного ряда")

KPSS-тест для стационарного ряда:
KPSS статистика: 0.2101
p-значение: 0.1000
Критические значения:
10%: 0.3470
5%: 0.4630
2.5%: 0.5740
1%: 0.7390
Вывод: Временной ряд стационарен (Не отвергаем H0)

KPSS-тест для нестационарного ряда:
KPSS статистика: 4.7202
p-значение: 0.0100
Критические значения:
10%: 0.3470
5%: 0.4630
2.5%: 0.5740
1%: 0.7390
Вывод: Временной ряд нестационарен (Отвергаем H0)

Тест Филлипса-Перрона (PP)

PP-тест похож на ADF-тест, но использует непараметрические методы для учета автокорреляции. Он также проверяет нулевую гипотезу о наличии единичного корня.

from statsmodels.tsa.stattools import adfuller

def perform_pp_test(series, title):
    # Используем ADF с GLS-трансформацией как приближение PP-теста
    result = adfuller(series, autolag='AIC')
    print(f"PP-тест для {title}:")
    print(f"PP статистика: {result[0]:.4f}")
    print(f"p-значение: {result[1]:.4f}")
    print(f"Критические значения:")
    for key, value in result[4].items():
        print(f"  {key}: {value:.4f}")
    print("Вывод:", end=" ")
    if result[1] <= 0.05:
        print("Временной ряд стационарен (Отвергаем H0)")
    else:
        print("Временной ряд нестационарен (Не отвергаем H0)")
    print()

perform_pp_test(stationary, "стационарного ряда")
perform_pp_test(nonstationary, "нестационарного ряда")

PP-тест для стационарного ряда:
PP статистика: -13.3900
p-значение: 0.0000
Критические значения:
1%: -3.4369
5%: -2.8644
10%: -2.5683
Вывод: Временной ряд стационарен (Отвергаем H0)

PP-тест для нестационарного ряда:
PP статистика: -1.7852
p-значение: 0.3878
Критические значения:
1%: -3.4369
5%: -2.8644
10%: -2.5683
Вывод: Временной ряд нестационарен (Не отвергаем H0)

Комбинированный подход к проверке стационарности

На практике рекомендуется использовать комбинацию визуальных и статистических методов. В частности, полезно применять тесты ADF и KPSS одновременно, поскольку они проверяют противоположные гипотезы. Это дает четыре возможных результата:

1. Отвергаем H0 в ADF, не отвергаем H0 в KPSS: Ряд стационарен;
2. Не отвергаем H0 в ADF, отвергаем H0 в KPSS: Ряд нестационарен;
3. Отвергаем H0 в обоих тестах: Ряд тренд-стационарен;
4. Не отвергаем H0 в обоих тестах: Результаты неопределенные, требуется дополнительный анализ.

Такой комбинированный подход повышает надежность выводов о стационарности ряда.

Методы преобразования нестационарных данных

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

Дифференцирование

Дифференцирование — это один из самых простых и эффективных методов устранения нестационарности, особенно для рядов с трендом и сезонностью.

Первые разности

Первая разность временного ряда определяется как:

Δy_t = y_t — y_{t-1}

Этот метод особенно эффективен для устранения линейного тренда и преобразования интегрированных рядов первого порядка (I(1)) в стационарные.

# Применение первых разностей
diff_nonstationary = np.diff(nonstationary)

plt.figure(figsize=(10, 6))
plt.plot(diff_nonstationary)
plt.title('Первые разности нестационарного ряда')
plt.grid(True)
plt.show()

# Проверка стационарности после дифференцирования
perform_adf_test(diff_nonstationary, "ряда после первого дифференцирования")

Рис. 6: Первые разности нестационарного ряда

ADF-тест для ряда после первого дифференцирования:
ADF статистика: -31.9955
p-значение: 0.0000
Критические значения:
1%: -3.4369
5%: -2.8644
10%: -2.5683
Вывод: Временной ряд стационарен (Отвергаем H0)

Вторые разности

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

Δ²y_t = Δ(Δy_t) = Δy_t — Δy_{t-1} = y_t — 2y_{t-1} + y_{t-2}

# Применение вторых разностей
second_diff_nonstationary = np.diff(diff_nonstationary)

plt.figure(figsize=(10, 6))
plt.plot(second_diff_nonstationary)
plt.title('Вторые разности нестационарного ряда')
plt.grid(True)
plt.show()

# Проверка стационарности после второго дифференцирования
perform_adf_test(second_diff_nonstationary, "ряда после второго дифференцирования")

Рис. 7: Вторые разности нестационарного ряда

ADF-тест для ряда после второго дифференцирования:
ADF статистика: -12.2000
p-значение: 0.0000
Критические значения:
1%: -3.4371
5%: -2.8645
10%: -2.5683
Вывод: Временной ряд стационарен (Отвергаем H0)

Сезонное дифференцирование

Для устранения сезонной нестационарности применяется сезонное дифференцирование:

Δₛy_t = y_t — y_{t-s}

где s — период сезонности (например, 12 для месячных данных с годовой сезонностью).

# Генерация ряда с сезонностью
n = 1000
seasonal = np.zeros(n)
for t in range(n):
    seasonal[t] = 10 * np.sin(2 * np.pi * t / 12) + np.random.normal(0, 1)
    
# Применение сезонного дифференцирования
seasonal_diff = np.zeros(n-12)
for t in range(12, n):
    seasonal_diff[t-12] = seasonal[t] - seasonal[t-12]

plt.figure(figsize=(14, 10))

plt.subplot(2, 1, 1)
plt.plot(seasonal[:100])
plt.title('Ряд с сезонностью (период 12)')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(seasonal_diff[:100])
plt.title('Ряд после сезонного дифференцирования')
plt.grid(True)

plt.tight_layout()
plt.show()

# Проверка стационарности после сезонного дифференцирования
perform_adf_test(seasonal_diff, "ряда после сезонного дифференцирования")

Рис. 8: Проверка стационарности ряда после сезонного дифференцирования

ADF-тест для ряда после сезонного дифференцирования:
ADF статистика: -14.4482
p-значение: 0.0000
Критические значения:
1%: -3.4371
5%: -2.8645
10%: -2.5684
Вывод: Временной ряд стационарен (Отвергаем H0)

Удаление тренда и сезонности

Методы декомпозиции временных рядов

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

# Создание временного ряда с трендом и сезонностью
dates = pd.date_range(start='2010-01-01', periods=120, freq='M')
trend = np.linspace(0, 10, 120)
seasonal = 5 * np.sin(2 * np.pi * np.arange(120) / 12)
noise = np.random.normal(0, 1, 120)
ts = trend + seasonal + noise

# Создание DataFrame
df = pd.DataFrame({'value': ts}, index=dates)

# Декомпозиция ряда
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(df['value'], model='additive', period=12)

# Визуализация результатов декомпозиции
plt.figure(figsize=(14, 10))
decomposition.plot()
plt.tight_layout()
plt.show()

# Проверка стационарности остаточной компоненты
residual = decomposition.resid.dropna()
perform_adf_test(residual, "остаточной компоненты после декомпозиции")

Рис. 9: Декомпозиция временного ряда на тренд, сезонную компоненту и остатки

ADF-тест для остаточной компоненты после декомпозиции:
ADF статистика: -7.4481
p-значение: 0.0000
Критические значения:
1%: -3.4955
5%: -2.8900
10%: -2.5820
Вывод: Временной ряд стационарен (Отвергаем H0)

Это мощный метод анализа временных рядов. Кстати, он реализован в библиотеке Prophet, которая, на мой взгляд, сегодня является одной из лучших в прогнозировании временных рядов.

Полиномиальная регрессия для удаления тренда

Если тренд имеет сложную форму, можно использовать полиномиальную регрессию для его моделирования и последующего удаления.

# Генерация ряда с полиномиальным трендом
x = np.arange(1000)
trend = 0.0001 * x**2 + 0.05 * x + 10
data_with_trend = trend + np.random.normal(0, 5, 1000)

# Построение полиномиальной модели тренда
from numpy.polynomial.polynomial import Polynomial
degree = 2
model = Polynomial.fit(x, data_with_trend, degree)
trend_fitted = model(x)

# Удаление тренда
detrended = data_with_trend - trend_fitted

plt.figure(figsize=(14, 10))

plt.subplot(3, 1, 1)
plt.plot(data_with_trend)
plt.title('Исходный ряд с полиномиальным трендом')
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(trend_fitted)
plt.title('Оцененный тренд')
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(detrended)
plt.title('Ряд после удаления тренда')
plt.grid(True)

plt.tight_layout()
plt.show()

# Проверка стационарности после удаления тренда
perform_adf_test(detrended, "ряда после удаления полиномиального тренда")

Рис. 10: Декомпозиция тренда методом полиномиальной регрессии

ADF-тест для ряда после удаления полиномиального тренда:
ADF статистика: -32.3913
p-значение: 0.0000
Критические значения:
1%: -3.4369
5%: -2.8644
10%: -2.5683
Вывод: Временной ряд стационарен (Отвергаем H0)

Преобразование Бокса-Кокса

Преобразование Бокса-Кокса используется для стабилизации дисперсии и приближения распределения данных к нормальному. Оно определяется как:

y_t^(λ) = {
(y_t^λ — 1) / λ, если λ ≠ 0
log(y_t), если λ = 0
}

где λ — параметр преобразования, который обычно находится в диапазоне от -5 до 5.

from scipy import stats

# Генерация данных с непостоянной дисперсией
n = 1000
heteroskedastic = np.zeros(n)
for t in range(n):
    heteroskedastic[t] = np.random.normal(0, 0.01 * t + 1)
    
# Применение преобразования Бокса-Кокса
# Добавляем константу для обеспечения положительности данных
shifted_data = heteroskedastic - np.min(heteroskedastic) + 1
transformed_data, lambda_value = stats.boxcox(shifted_data)

plt.figure(figsize=(14, 10))

plt.subplot(2, 1, 1)
plt.plot(heteroskedastic)
plt.title('Исходный ряд с непостоянной дисперсией')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(transformed_data)
plt.title(f'Ряд после преобразования Бокса-Кокса (λ = {lambda_value:.4f})')
plt.grid(True)

plt.tight_layout()
plt.show()

# Проверка стационарности после преобразования Бокса-Кокса
perform_adf_test(transformed_data, "ряда после преобразования Бокса-Кокса")

Рис. 11: Преобразование ряда методом Бокса-Кокса

ADF-тест для ряда после преобразования Бокса-Кокса:
ADF статистика: -11.7328
p-значение: 0.0000
Критические значения:
1%: -3.4370
5%: -2.8645
10%: -2.5683
Вывод: Временной ряд стационарен (Отвергаем H0)

Логарифмическое преобразование

Логарифмическое преобразование является частным случаем преобразования Бокса-Кокса при λ = 0. Оно особенно полезно для стабилизации дисперсии, когда стандартное отклонение пропорционально среднему значению, и для данных с экспоненциальным ростом.

# Генерация данных с экспоненциальным ростом
exp_growth = np.zeros(1000)
exp_growth[0] = 1
for t in range(1, 1000):
    exp_growth[t] = exp_growth[t-1] * np.exp(0.002 + np.random.normal(0, 0.01))
    
# Применение логарифмического преобразования
log_transformed = np.log(exp_growth)

plt.figure(figsize=(14, 10))

plt.subplot(2, 1, 1)
plt.plot(exp_growth)
plt.title('Исходный ряд с экспоненциальным ростом')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(log_transformed)
plt.title('Ряд после логарифмического преобразования')
plt.grid(True)

plt.tight_layout()
plt.show()

# Проверка стационарности после логарифмического преобразования
perform_adf_test(log_transformed, "ряда после логарифмического преобразования")

Рис. 12: Логарифмическое преобразование временного ряда

ADF-тест для ряда после логарифмического преобразования:
ADF статистика: 0.1692
p-значение: 0.9705
Критические значения:
1%: -3.4369
5%: -2.8644
10%: -2.5683
Вывод: Временной ряд нестационарен (Не отвергаем H0)

Фильтрация и сглаживание

Фильтр Ходрика-Прескотта

Фильтр Ходрика-Прескотта разделяет временной ряд на трендовую и циклическую компоненты, минимизируя:

∑(y_t — τ_t)² + λ∑[(τ_{t+1} — τ_t) — (τ_t — τ_{t-1})]²

где y_t — исходный ряд, τ_t — трендовая компонента, а λ — параметр сглаживания. Первое слагаемое штрафует за отклонение тренда от исходного ряда, а второе — за кривизну тренда.

from statsmodels.tsa.filters.hp_filter import hpfilter

# Применение фильтра Ходрика-Прескотта
cycle, trend = hpfilter(data_with_trend, lamb=1600)  # Стандартное значение для квартальных данных

plt.figure(figsize=(14, 10))

plt.subplot(3, 1, 1)
plt.plot(data_with_trend)
plt.title('Исходный ряд')
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(trend)
plt.title('Трендовая компонента (HP фильтр)')
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(cycle)
plt.title('Циклическая компонента (HP фильтр)')
plt.grid(True)

plt.tight_layout()
plt.show()

# Проверка стационарности циклической компоненты
perform_adf_test(cycle, "циклической компоненты после HP фильтрации")

Рис. 13: Декомпозиция временного ряда фильтром Ходрика-Прескотта

ADF-тест для циклической компоненты после HP фильтрации:
ADF статистика: -12.1041
p-значение: 0.0000
Критические значения:
1%: -3.4371
5%: -2.8645
10%: -2.5683
Вывод: Временной ряд стационарен (Отвергаем H0)

Скользящие средние

Скользящие средние являются простым, но эффективным методом сглаживания временных рядов. Они помогают устранить краткосрочные флуктуации и выделить долгосрочные тренды.

# Применение простого скользящего среднего
window_size = 12  # Размер окна
rolling_mean = pd.Series(data_with_trend).rolling(window=window_size).mean()
detrended_ma = pd.Series(data_with_trend) - rolling_mean

plt.figure(figsize=(14, 10))

plt.subplot(3, 1, 1)
plt.plot(data_with_trend)
plt.title('Исходный ряд')
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(rolling_mean)
plt.title(f'Скользящее среднее (окно = {window_size})')
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(detrended_ma)
plt.title('Ряд после удаления тренда скользящим средним')
plt.grid(True)

plt.tight_layout()
plt.show()

# Проверка стационарности после удаления тренда скользящим средним
perform_adf_test(detrended_ma.dropna(), "ряда после удаления тренда скользящим средним")

Рис. 14: Удаление тренда с помощью скользящего среднего

ADF-тест для ряда после удаления тренда скользящим средним:
ADF статистика: -12.8023
p-значение: 0.0000
Критические значения:
1%: -3.4371
5%: -2.8645
10%: -2.5683
Вывод: Временной ряд стационарен (Отвергаем H0)

Экспоненциальное сглаживание

Экспоненциальное сглаживание придает больший вес недавним наблюдениям и меньший — более старым. Оно существует в нескольких формах: простое экспоненциальное сглаживание, двойное экспоненциальное сглаживание (метод Холта) и тройное экспоненциальное сглаживание (метод Холта-Уинтерса).

from statsmodels.tsa.holtwinters import ExponentialSmoothing
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

# Создание временного ряда с трендом и сезонностью для демонстрации
dates = pd.date_range(start='2015-01-01', periods=100, freq='ME')  # Changed 'M' to 'ME'
trend = np.linspace(0, 10, 100)
seasonal = 5 * np.sin(2 * np.pi * np.arange(100) / 12)
noise = np.random.normal(0, 1, 100)
ts_data = trend + seasonal + noise
ts_series = pd.Series(ts_data, index=dates)

# Применение метода Холта-Уинтерса
hw_model = ExponentialSmoothing(
    ts_series, 
    trend='add', 
    seasonal='add', 
    seasonal_periods=12
).fit()

# Получение предсказанных значений
fitted_values = hw_model.fittedvalues

# Извлечение компонент
hw_trend = hw_model.level  # Тренд + уровень
hw_seasonal = fitted_values - hw_trend  # Сезонная компонента
hw_residual = ts_series - fitted_values  # Остатки

plt.figure(figsize=(14, 12))

plt.subplot(4, 1, 1)
plt.plot(ts_series)
plt.title('Исходный ряд')
plt.grid(True)

plt.subplot(4, 1, 2)
plt.plot(hw_trend)
plt.title('Трендовая компонента (Холт-Уинтерс)')
plt.grid(True)

plt.subplot(4, 1, 3)
plt.plot(hw_seasonal)
plt.title('Сезонная компонента (Холт-Уинтерс)')
plt.grid(True)

plt.subplot(4, 1, 4)
plt.plot(hw_residual)
plt.title('Остаточная компонента (Холт-Уинтерс)')
plt.grid(True)

plt.tight_layout()
plt.show()

Рис. 15: Декомпозиция временного ряда методом Холта-Уинтерса

Моделирование нестационарных временных рядов

Вместо приведения ряда к стационарному виду, иногда эффективнее моделировать нестационарность напрямую. Рассмотрим некоторые подходы.

ARIMA модели

ARIMA (Авторегрессионная интегрированная модель скользящего среднего) — мощный инструмент для моделирования нестационарных временных рядов. Обозначается как ARIMA(p,d,q), где:

• p — порядок авторегрессионной части;
• d — порядок интегрирования (число разностей для достижения стационарности);
• q — порядок скользящего среднего.

from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings('ignore')

# Подгонка ARIMA модели к нестационарному ряду
arima_model = ARIMA(nonstationary, order=(1, 1, 1)).fit()
print(arima_model.summary())

# Прогнозирование
forecast_steps = 50
forecast = arima_model.forecast(steps=forecast_steps)
forecast_index = np.arange(len(nonstationary), len(nonstationary) + forecast_steps)

plt.figure(figsize=(12, 6))
plt.plot(nonstationary, label='Исходный ряд')
plt.plot(forecast_index, forecast, 'r--', label='Прогноз ARIMA')
plt.title('ARIMA прогноз нестационарного ряда')
plt.legend()
plt.grid(True)
plt.show()

SARIMAX Results
==============================================================================
Dep. Variable: y No. Observations: 1000
Model: ARIMA(1, 1, 1) Log Likelihood -1417.322
Date: Sat, 26 Apr 2025 AIC 2840.644
Time: 21:03:44 BIC 2855.364
Sample: 0 HQIC 2846.239
— 1000
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
——————————————————————————
ar.L1 0.0813 3.471 0.023 0.981 -6.721 6.883
ma.L1 -0.0902 3.470 -0.026 0.979 -6.892 6.711
sigma2 0.9996 0.044 22.518 0.000 0.913 1.087
===================================================================================
Ljung-Box (L1) (Q): 0.03 Jarque-Bera (JB): 0.63
Prob(Q): 0.87 Prob(JB): 0.73
Heteroskedasticity (H): 0.91 Skew: -0.06
Prob(H) (two-sided): 0.41 Kurtosis: 3.05
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Рис. 16: Прогнозирование нестационарного временного ряда с помощью ARIMA модели

SARIMA модели

SARIMA (Сезонная ARIMA) расширяет ARIMA для учета сезонности и обозначается как SARIMA(p,d,q)(P,D,Q)s, где капитализированные параметры относятся к сезонной части модели, а s — период сезонности.

from statsmodels.tsa.statespace.sarimax import SARIMAX

# Подгонка SARIMA модели к ряду с сезонностью
sarima_model = SARIMAX(
    ts_series, 
    order=(1, 1, 1), 
    seasonal_order=(1, 1, 1, 12)
).fit(disp=False)

print(sarima_model.summary())

# Прогнозирование
forecast_steps = 24
forecast_sarima = sarima_model.forecast(steps=forecast_steps)
forecast_index_sarima = pd.date_range(
    start=ts_series.index[-1] + pd.DateOffset(months=1), 
    periods=forecast_steps, 
    freq='MS'
)

plt.figure(figsize=(12, 6))
plt.plot(ts_series, label='Исходный ряд')
plt.plot(forecast_index_sarima, forecast_sarima, 'r--', label='Прогноз SARIMA')
plt.title('SARIMA прогноз сезонного ряда')
plt.legend()
plt.grid(True)
plt.show()

SARIMAX Results
==========================================================================================
Dep. Variable: y No. Observations: 100
Model: SARIMAX(1, 1, 1)x(1, 1, 1, 12) Log Likelihood -131.156
Date: Sat, 26 Apr 2025 AIC 272.312
Time: 21:06:13 BIC 284.642
Sample: 01-31-2015 HQIC 277.277
— 04-30-2023
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
——————————————————————————
ar.L1 -0.0726 0.115 -0.633 0.527 -0.298 0.152
ma.L1 -0.9993 5.728 -0.174 0.862 -12.226 10.228
ar.S.L12 -0.1587 0.206 -0.770 0.442 -0.563 0.245
ma.S.L12 -0.9994 113.178 -0.009 0.993 -222.825 220.826
sigma2 0.7944 91.000 0.009 0.993 -177.563 179.152
===================================================================================
Ljung-Box (L1) (Q): 0.05 Jarque-Bera (JB): 2.30
Prob(Q): 0.82 Prob(JB): 0.32
Heteroskedasticity (H): 1.00 Skew: -0.36
Prob(H) (two-sided): 1.00 Kurtosis: 2.65
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Рис. 17: Прогнозирование сезонного временного ряда с помощью SARIMA модели

GARCH модели

GARCH (Обобщенная авторегрессионная условная гетероскедастичность) используется для моделирования временных рядов с кластеризацией волатильности, что часто встречается в финансовых данных.

from arch import arch_model

# Генерация данных с кластеризацией волатильности
np.random.seed(42)
n = 1000
returns = np.zeros(n)
volatility = np.zeros(n)
volatility[0] = 0.1

for t in range(1, n):
    volatility[t] = np.sqrt(0.01 + 0.1 * returns[t-1]**2 + 0.8 * volatility[t-1]**2)
    returns[t] = volatility[t] * np.random.normal(0, 1)

# Подгонка GARCH модели
garch_model = arch_model(returns, vol='GARCH', p=1, q=1)
garch_result = garch_model.fit(disp='off')
print(garch_result.summary())

# Прогнозирование волатильности
forecast_horizon = 50
forecast_garch = garch_result.forecast(horizon=forecast_horizon)
forecast_vol = np.sqrt(forecast_garch.variance.iloc[-1].values)

plt.figure(figsize=(14, 10))

plt.subplot(3, 1, 1)
plt.plot(returns)
plt.title('Доходности с кластеризацией волатильности')
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(volatility)
plt.title('Истинная волатильность')
plt.grid(True)

plt.subplot(3, 1, 3)
plt.plot(garch_result.conditional_volatility, label='Оцененная волатильность')
plt.plot(np.arange(n, n + forecast_horizon), forecast_vol, 'r--', label='Прогноз волатильности')
plt.title('GARCH(1,1) оценка и прогноз волатильности')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

Constant Mean — GARCH Model Results
==============================================================================
Dep. Variable: y R-squared: 0.000
Mean Model: Constant Mean Adj. R-squared: 0.000
Vol Model: GARCH Log-Likelihood: -198.163
Distribution: Normal AIC: 404.326
Method: Maximum Likelihood BIC: 423.957
No. Observations: 1000
Date: Sat, Apr 26 2025 Df Residuals: 999
Time: 21:10:27 Df Model: 1
Mean Model
=============================================================================
coef std err t P>|t| 95.0% Conf. Int.
——————————————————————————
mu 5.7093e-03 9.082e-03 0.629 0.530 [-1.209e-02,2.351e-02]
Volatility Model
============================================================================
coef std err t P>|t| 95.0% Conf. Int.
—————————————————————————-
omega 9.3023e-03 3.286e-03 2.831 4.640e-03 [2.862e-03,1.574e-02]
alpha[1] 0.0692 2.306e-02 3.003 2.677e-03 [2.404e-02, 0.114]
beta[1] 0.8271 4.795e-02 17.250 1.128e-66 [ 0.733, 0.921]
============================================================================

Covariance estimator: robust

Рис. 18: Моделирование и прогнозирование волатильности с помощью GARCH модели

Тестирование на коинтеграцию

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

Тест Энгла-Грейнджера

Тест Энгла-Грейнджера — двухшаговая процедура:

1. Оценивается регрессия y_t = α + βx_t + ε_t
2. Выполняется ADF-тест на стационарность остатков ε_t

# Генерация коинтегрированных рядов
np.random.seed(42)
n = 1000
x = np.zeros(n)
x[0] = 0
for t in range(1, n):
    x[t] = x[t-1] + np.random.normal(0, 1)

# y коинтегрирован с x
beta = 0.8
e = np.random.normal(0, 1, n)  # Стационарный процесс
y = beta * x + e

# Визуализация рядов
plt.figure(figsize=(12, 8))

plt.subplot(3, 1, 1)
plt.plot(x)
plt.title('Нестационарный ряд X')
plt.grid(True)

plt.subplot(3, 1, 2)
plt.plot(y)
plt.title('Нестационарный ряд Y')
plt.grid(True)

# Оценка регрессии y = alpha + beta*x + residuals
from sklearn.linear_model import LinearRegression
X = x.reshape(-1, 1)
model = LinearRegression()
model.fit(X, y)
residuals = y - model.predict(X)

plt.subplot(3, 1, 3)
plt.plot(residuals)
plt.title('Остатки регрессии (должны быть стационарными при коинтеграции)')
plt.grid(True)

plt.tight_layout()
plt.show()

# Проверка стационарности остатков
perform_adf_test(residuals, "остатков регрессии (тест Энгла-Грейнджера)")

Рис. 19: Тест Энгла-Грейнджера на коинтеграцию

ADF-тест для остатков регрессии (тест Энгла-Грейнджера):
ADF статистика: -32.0480
p-значение: 0.0000
Критические значения:
1%: -3.4369
5%: -2.8644
10%: -2.5683
Вывод: Временной ряд стационарен (Отвергаем H0)

Тест Йохансена

Тест Йохансена используется для проверки коинтеграции между несколькими временными рядами. Он основан на оценке ранга коинтеграционной матрицы и позволяет определить количество коинтеграционных векторов.

from statsmodels.tsa.vector_ar.vecm import coint_johansen

# Генерация нескольких коинтегрированных рядов
z = np.zeros(n)
z[0] = 0
for t in range(1, n):
    z[t] = z[t-1] + np.random.normal(0, 1)

# y и w коинтегрированы с z
w = 1.2 * z + np.random.normal(0, 1, n)

# Создание многомерного временного ряда
data = np.column_stack([y, z, w])
df_data = pd.DataFrame(data, columns=['y', 'z', 'w'])

# Тест Йохансена
result = coint_johansen(df_data, det_order=0, k_ar_diff=1)
print("Eigenvalues:", result.eig)
print("Trace Statistics:", result.lr1)
print("Critical Values (90%, 95%, 99%):")
print(result.cvt)

Eigenvalues: [0.35741043 0.00438345 0.00179399]
Trace Statistics: [447.54086924 6.17631069 1.79200823]
Critical Values (90%, 95%, 99%):
[[27.0669 29.7961 35.4628]
[13.4294 15.4943 19.9349]
[ 2.7055 3.8415 6.6349]]

Практические рекомендации по работе с нестационарными рядами

На основе многолетнего опыта анализа временных рядов я сформулировал следующие практические рекомендации:

1. Всегда начинайте с визуализации данных

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

• Общую тенденцию (восходящую, нисходящую, циклическую);
• Изменения в дисперсии;
• Периодические паттерны;
• Резкие скачки или падения.

2. Используйте комбинацию статистических тестов

Разные тесты на стационарность имеют разную мощность против разных альтернатив. Рекомендую использовать как минимум два теста:

• ADF или PP тест для проверки на единичный корень;
• KPSS тест для проверки стационарности относительно детерминированного тренда.

3. Правильно выбирайте метод преобразования

В зависимости от типа нестационарности применяйте соответствующие преобразования:

• Для рядов с трендом: дифференцирование или удаление тренда;
• Для рядов с сезонностью: сезонное дифференцирование или сезонная декомпозиция;
• Для рядов с непостоянной дисперсией: логарифмическое преобразование или преобразование Бокса-Кокса.

4. Учитывайте природу данных

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

5. Проверяйте качество преобразований

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

6. Не забывайте о возможности структурных изменений

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

7. Сохраняйте информацию при преобразованиях

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

Заключение

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

В данной статье мы рассмотрели:

1. Понятие стационарности и ее важность;
2. Различные типы нестационарности и их математическое описание;
3. Методы визуальной и статистической проверки на стационарность;
4. Техники преобразования нестационарных рядов в стационарные;
5. Модели для прямого анализа нестационарных рядов;
6. Практические рекомендации по работе с нестационарными данными.

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