Глава 12. Анализ временных рядов
Временной ряд — это последовательность измерений некоторой системы, меняющаяся во времени. Многие инструменты, которые мы использовали в предыдущих главах, такие как регрессия, также могут быть использованы для работы с временны́ми рядами. Но существуют и отдельные методы, особенно полезные для такого рода данных.
Для примера возьмем два датасета — с данными о производстве электроэнергии из возобновляемых источников в США с 2001 по 2024 год и с погодными данными за тот же период. Мы разработаем методы для разделения временного ряда на долгосрочный тренд и повторяющуюся сезонную составляющую. Для моделирования и прогнозирования трендов будем использовать линейную регрессию. Мы также опробуем широко используемую модель для анализа данных временных рядов, которая носит название «интегрированная модель авторегрессии скользящего среднего» (autoregressive integrated moving average), или сокращенно ARIMA.
Данные об электрогенерации
В качестве примера данных временных рядов используем датасет Управления энергетической информации (Energy Information Administration) США. В нем содержатся данные о месячном объеме выработки электроэнергии с использованием возобновляемых источников за период с 2001 по 2024 год. Инструкции по загрузке данных, как обычно, приведены в Jupyter-блокноте этой главы.
После загрузки данные необходимо преобразовать, чтобы привести в удобный для работы формат:
elec = (
pd.read_csv("Net_generation_for_all_sectors.csv", skiprows=4)
.drop(columns=["units", "source key"])
.set_index("description")
.replace("--", np.nan)
.transpose()
.astype(float)
)
В отформатированном датасете каждый столбец представляет собой последовательность суммарных месячных показателей выработки в гигаватт-часах (ГВт-ч). Вот названия столбцов для разных источников электроэнергии, или «секторов»:
elec.columns
Index(['Net generation for all sectors (Чистая выработка во всех секторах)',
'United States (США)',
'United States : all fuels (utility-scale) (США: все виды ископаемого
топлива (промышленного масштаба))', 'United States : nuclear (США:
атомная энергетика)',
'United States : conventional hydroelectric (США: традиционная
гидроэнергетика)',
'United States : other renewables (США: другие возобновляемые источники
энергии)', 'United States : wind (США: ветрогенерация)',
'United States : all utility-scale solar (США: вся солнечная генерация
промышленного масштаба)', 'United States : geothermal (США:
геотермальная генерация)',
'United States : biomass (США: биомасса)',
'United States : hydro-electric pumped storage (США: гидроаккумулирующие
электростанции)',
'United States : all solar (США: вся солнечная генерация)',
'United States : small-scale solar photovoltaic (США: малая солнечная
фотоэлектрическая генерация'],
dtype='object', name='description')
Метки индекса содержат месяц и год. Вот первые двенадцать из них:
elec.index[:12]
Index(['Jan (Январь) 2001', 'Feb (Февраль) 2001', 'Mar (Март) 2001',
'Apr (Апрель) 2001', 'May (Май) 2001', 'Jun (Июнь) 2001',
'Jul (Июль) 2001', 'Aug (Август) 2001', 'Sep (Сентябрь) 2001',
'Oct (Октябрь) 2001', 'Nov (Ноябрь) 2001', 'Dec (Декабрь) 2001'],
dtype='object')
Будет проще работать с этими данными, если заменить строки объектами Timestamp библиотеки Pandas. Сгенерируем с помощью функции date_range последовательность объектов Timestamp начиная с января 2001 года. Частотный код "ME", расшифровывающийся как «конец месяца» (month end), указывает этой функции, что необходимо подставить последний день каждого месяца:
elec.index = pd.date_range(start="2001-01", periods=len(elec), freq="ME")
elec.index[:6]
DatetimeIndex(['2001-01-31', '2001-02-28', '2001-03-31', '2001-04-30',
'2001-05-31', '2001-06-30'],
dtype='datetime64[ns]', freq='ME')
Теперь индекс представляет собой объект DataTimeIndex с типом данных datetime64[ns], определенным в библиотеке NumPy. 64 означает, что каждая метка использует 64 бита, а ns — что она имеет наносекундную точность.
Декомпозиция
В качестве первого примера рассмотрим, как менялась выработка электроэнергии атомными электростанциями (АЭС) за период с января 2001 года по июнь 2024 года, и разложим временные ряды на долгосрочный тренд и периодическую составляющую. Ниже представлен график общего объема месячной выработки электроэнергии на АЭС в США:
nuclear = elec["United States : nuclear"]
nuclear.plot(label="атомная генерация", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Видно, что есть подъемы и спады, но их трудно рассмотреть, поскольку наблюдаются большие изменения от месяца к месяцу. Чтобы лучше понять долгосрочный тренд, можно воспользоваться методами rolling и mean для вычисления скользящего среднего:
trend = nuclear.rolling(window=12).mean()
Параметр размера окна window=12 задает перекрывающиеся интервалы в 12 месяцев: первый интервал охватывает 12 измерений, начиная с первого, второй интервал — начиная со второго, и т.д. Для каждого интервала мы вычисляем средний объем выработки.
Вот как выглядит кривая тренда вместе с исходными данными:
nuclear.plot(label="атомная генерация", **actual_options)
trend.plot(label="тренд", **trend_options)
decorate(ylabel="Выработка (ГВт-ч)")
Тренд все еще довольно изменчив. Мы могли бы сгладить его еще сильнее, увеличив параметр window, но пока оставим 12-месячное окно.
Если вычесть тренд из исходных данных, то в результате получится временной ряд с исключенным трендом (detrended), это подразумевает, что его долгосрочное среднее значение почти является константой. Вот как он выглядит:
detrended = (nuclear - trend).dropna()
detrended.plot(label="за вычетом тренда", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Похоже, существует повторяющаяся ежегодная закономерность. И в этом есть логика, поскольку спрос на электроэнергию зависит от сезона: зимой она используется для обогрева, а летом — для кондиционирования воздуха. Чтобы описать эту закономерность, можно взять месячную часть объектов datetime в индексе, сгруппировать данные по месяцам и вычислить средний объем производства. Вот как выглядят среднемесячные значения:
monthly_averages = detrended.groupby(detrended.index.month).mean()
monthly_averages.plot(label="среднемесячное значение", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
На оси x отложены месяцы с января (1) по декабрь (12). Генерация электроэнергии максимальна в самые холодные и самые теплые месяцы года, минимальна в апреле и октябре.
Мы можем использовать среднемесячные значения из monthly_averages для вычисления сезонной составляющей данных. Она представляет собой ряд той же длины, что и в переменной nuclear. Элементы ряда для каждого месяца равны среднему значению за этот месяц. Вот как выглядит сезонная составляющая:
seasonal = monthly_averages[nuclear.index.month]
seasonal.index = nuclear.index
seasonal.plot(label="сезонная составляющая", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Каждый промежуток в 12 месяцев идентичен остальным.
Сумма тренда и сезонной составляющей представляет собой ожидаемое (expected) значение для каждого месяца:
expected = trend + seasonal
Вот как она выглядит в сравнении с оригинальным рядом:
expected.plot(label="ожидаемое значение", **pred_options)
nuclear.plot(label="фактическое значение", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Если вычесть эту сумму из исходного ряда, результатом будет остаток (residual), представляющий собой отклонение от ожидаемого значения для каждого месяца:
resid = nuclear - expected
resid.plot(label="остаток", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Остаток можно считать суммой всех возможных факторов, которые влияют на производство энергии, но не объясняются долгосрочным трендом или сезонной составляющей. Среди прочего, эта сумма включает погоду, оборудование, находящееся на техобслуживании, и изменение спроса в связи с конкретными событиями. Поскольку остаток — это сумма многих непредсказуемых, а иногда и неизвестных факторов, его часто рассматривают как случайную величину.
Вот как выглядит распределение остатков:
from thinkstats import plot_kde
plot_kde(resid.dropna())
decorate(xlabel="Остаток (ГВт-ч)")
Оно напоминает колоколообразную кривую нормального распределения, что согласуется с допущением, что остаток является суммой множества случайных вкладов.
Чтобы количественно оценить, насколько хорошо эта модель описывает исходный ряд, вычислим коэффициент детерминации, который показывает, насколько дисперсия остатков меньше по сравнению с дисперсией исходного ряда:
rsquared = 1 - resid.var() / nuclear.var()
rsquared
0.9054559977517084
Коэффициент R2 около 0.91 означает, что на долгосрочный тренд и сезонную составляющую приходится 91 % изменчивости ряда. В данном случае R2 существенно выше, чем значения, которые нам встретились в предыдущей главе. Это характерно для данных временных рядов, особенно в сценариях, подобных этому, когда мы построили модель так, чтобы она имела сходство с данными.
Процедура, которую мы только что рассмотрели, называется сезонной декомпозицией. В пакете StatsModels есть функция seasonal_decompose, которая выполняет эти действия:
from statsmodels.tsa.seasonal import seasonal_decompose
decomposition = seasonal_decompose(nuclear, model="additive", period=12)
Аргумент model="additive" задает аддитивную модель, поэтому ряд разбивается на сумму тренда, сезонной составляющей и остатка. Вскоре вы также познакомитесь с мультипликативной моделью. Аргумент period=12 задает продолжительность сезонного компонента в 12 месяцев.
Функция возвращает содержащий эти три составляющих объект. В Jupyter-блокноте этой главы реализована функция, которая строит их графики:
plot_decomposition(nuclear, decomposition)
Результаты аналогичны тем, что мы получили самостоятельно, а небольшие отличия можно объяснить особенностями реализации.
Такая сезонная декомпозиция позволяет получить представление о структуре временного ряда. Как мы увидим в следующем разделе, она также полезна в построении прогнозов.
Предсказание данных
Результаты сезонной декомпозиции можно использовать для предсказания будущего. В демонстративных целях с помощью функции ниже разделим временной ряд на обучающий (training series) — его мы будем использовать для генерации предсказаний — и тестовый (test series), который будем использовать для проверки точности предсказаний:
def split_series(series, n=60):
training = series.iloc[:-n]
test = series.iloc[-n:]
return training, test
При n=60 длина тестового ряда составляет пять лет, начиная с июля 2019 года:
training, test = split_series(nuclear)
test.index[0]
Timestamp('2019-07-31 00:00:00')
Теперь представьте, что сейчас июнь 2019 года и вас просят составить пятилетний прогноз производства электроэнергии на АЭС. Чтобы ответить на этот вопрос, сначала с помощью обучающих данных построим модель, а затем воспользуемся ею для генерации предсказаний. Начнем с сезонной декомпозиции обучающих данных:
decomposition = seasonal_decompose(training, model="additive", period=12)
trend = decomposition.trend
Теперь применим к тренду линейную модель. Объясняющая переменная months равна количеству месяцев от начала ряда:
import statsmodels.formula.api as smf
months = np.arange(len(trend))
data = pd.DataFrame({"trend": trend, "months": months}).dropna()
results = smf.ols("trend ~ months", data=data).fit()
Вот сводная таблица результатов:
from thinkstats import display_summary
display_summary(results)
| coef | std err | t | P > | t | | [0.025 | 0.975] | |
| Intercept | 6.482e+04 | 131.524 | 492.869 | 0.000 | 6.46e+04 | 6.51e+04 |
| months | 10.9886 | 1.044 | 10.530 | 0.000 | 8.931 | 13.046 |
R-squared: 0.3477
Величина коэффициента R2 около 0.35 свидетельствует, что модель не очень хорошо соответствует данным. Мы сможем лучше в этом разобраться, если построим линию регрессии. Воспользуемся методом predict для вычисления ожидаемых значений для обучающих и тестовых данных:
months = np.arange(len(training) + len(test))
df = pd.DataFrame({"months": months})
pred_trend = results.predict(df)
pred_trend.index = nuclear.index
Вот график тренда и линия регрессии:
trend.plot(**trend_options)
pred_trend.plot(label="линейная модель", **model_options)
decorate(ylabel="Выработка (ГВт-ч)")
Линейная модель не учитывает многих деталей, но, похоже, в целом наблюдается тенденция к повышению выработки.
Далее воспользуемся сезонной составляющей из результатов декомпозиции, чтобы получить объект Series со среднемесячными значениями (monthly averages):
seasonal = decomposition.seasonal
monthly_averages = seasonal.groupby(seasonal.index.month).mean()
Сезонную составляющую можно предсказать, если по датам линии регрессии выбрать элементы из последовательности monthly_averages:
pred_seasonal = monthly_averages[pred_trend.index.month]
pred_seasonal.index = pred_trend.index
Наконец, чтобы сгенерировать предсказания, добавим к тренду сезонную составляющую:
pred = pred_trend + pred_seasonal
Вот обучающие данные и полученные предсказания:
pred.plot(label="предсказание", **pred_options)
training.plot(label="обучающие данные", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Предсказанные значения достаточно хорошо согласуются с обучающими данными, и прогноз выглядит обоснованным в предположении, что долгосрочный тренд сохранится.
Теперь из позиции в будущем посмотрим, насколько точным оказался наш прогноз. Вот предсказанные и фактические значения за пятилетний период с июля 2019 года:
forecast = pred[test.index]
forecast.plot(label="предсказанное значение", **pred_options)
test.plot(label="фактическое значение", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Прогноз для первого года довольно неплох, но производство на АЭС в 2020 году оказалось ниже, чем ожидалось — возможно, из-за пандемии COVID-19, — и оно так и не вернулось к уровню долгосрочного тренда.
Чтобы количественно оценить точность предсказаний, используем среднюю абсолютную процентную ошибку (mean absolute percentage error, MAPE), которую вычисляет следующая функция:
def MAPE(predicted, actual):
ape = np.abs(predicted - actual) / actual
return np.mean(ape) * 100
В этом примере предсказания в среднем смещены на 3.81 %:
MAPE(forecast, test)
3.811940747879257
Мы вернемся к этому примеру позже в этой главе и посмотрим, можно ли улучшить результаты с помощью другой модели.
Мультипликативная модель
В аддитивной модели, которую мы использовали в предыдущем разделе, предполагается, что временной ряд представляет собой сумму долгосрочного тренда, сезонной составляющей и остатка. Это подразумевает, в свою очередь, что величина сезонной составляющей и остатков не меняется с течением времени.
В качестве примера того, как нарушается это предположение, рассмотрим объемы малой (small-scale) солнечной генерации электроэнергии начиная с 2014 года:
solar = elec["United States : small-scale solar photovoltaic"].dropna()
solar.plot(label="солнечная генерация", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
За этот период общий объем выработки увеличился в несколько раз. Амплитуда сезонных колебаний, очевидно, также возросла.
Предположение о том, что амплитуда сезонных колебаний и случайной изменчивости пропорциональна величине тренда, наводит на мысль о модели, альтернативной аддитивной, в которой временной ряд является произведением трех составляющих.
Чтобы опробовать эту мультипликативную модель, разделим этот ряд на обучающий и тестовый наборы:
training, test = split_series(solar)
И вызовем функцию seasonal_decompose, указывая аргумент model="multiplicative":
decomposition = seasonal_decompose(training, model="multiplicative", period=12)
Вот как выглядят полученные составляющие:
plot_decomposition(training, decomposition)
Теперь сезонная составляющая и остаток являются мультипликативными множителями. И картина такова, что сезонная составляющая колеблется примерно от уровня 25 % ниже тренда до 25 % выше. Независимо от этого остаток обычно составляет менее 5 %, за исключением небольшого временного отрезка в начальный период. Составляющие модели можно извлечь следующим образом:
trend = decomposition.trend
seasonal = decomposition.seasonal
resid = decomposition.resid
Коэффициент R2 у этой модели очень высок:
rsquared = 1 - resid.var() / training.var()
rsquared
0.9999999992978134
Солнечная генерация главным образом зависит от количества падающего на панели солнечного света, поэтому вполне логично, что выработка так строго следует годовому циклу. Для предсказания долгосрочного тренда используем квадратичную модель:
months = range(len(training))
data = pd.DataFrame({"trend": trend, "months": months}).dropna()
results = smf.ols("trend ~ months + I(months**2)", data=data).fit()
В формульной записи пакета Patsy подстрока I(months**2) добавляет к модели квадратичный член, поэтому его не нужно вычислять явно. Вот результаты подгонки модели:
display_summary(results)
| coef | std err | t | P > | t | | [0.025 | 0.975] | |
| Intercept | 766.1962 | 13.494 | 56.782 | 0.000 | 739.106 | 793.286 |
| months | 22.2153 | 0.938 | 23.673 | 0.000 | 20.331 | 24.099 |
| I(months ** 2) | 0.1762 | 0.014 | 12.480 | 0.000 | 0.148 | 0.205 |
R-squared: 0.9983
У линейного и квадратичного коэффициентов p-значения очень малы, что позволяет предположить, что квадратичная модель отражает больше информации о тренде, чем линейная. Коэффициент R2 также очень велик.
Теперь с помощью этой модели можно посчитать ожидаемые значения тренда для моментов времени в прошлом и будущем:
months = range(len(solar))
df = pd.DataFrame({"months": months})
pred_trend = results.predict(df)
pred_trend.index = solar.index
Вот как выглядят полученные результаты:
pred_trend.plot(label="квадратичная модель", **model_options)
trend.plot(**trend_options)
decorate(ylabel="Выработка (ГВт-ч)")
Квадратичная модель хорошо передает тенденцию в прошлом. Теперь можно использовать сезонную составляющую для предсказания сезонной изменчивости в будущем:
monthly_averages = seasonal.groupby(seasonal.index.month).mean()
pred_seasonal = monthly_averages[pred_trend.index.month]
pred_seasonal.index = pred_trend.index
Наконец, чтобы рассчитать ретроспективные оценки (retrodiction) для значений в прошлом и предсказания для будущего, перемножим тренд и сезонную составляющую:
pred = pred_trend * pred_seasonal
Вот полученные результаты на одном графике с обучающими данными:
training.plot(label="обучающие данные", **actual_options)
pred.plot(label="предсказание", **pred_options)
decorate(ylabel="Выработка (ГВт-ч)")
Ретроспективные оценки хорошо согласуются с обучающими данными, а прогноз кажется правдоподобным. Теперь посмотрим, насколько он оказался верным. Вот предсказанные значения вместе с тестовыми данными:
future = pred[test.index]
future.plot(label="предсказанное значение", **pred_options)
test.plot(label="фактическое значение", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Предсказания для первых трех лет очень хорошо соответствуют данным. В дальнейшем, как видно, фактический рост превысил ожидания.
В этом примере сезонная декомпозиция хорошо показала себя при моделировании и прогнозировании солнечной генерации, но в предыдущем примере с АЭС она была не очень эффективна. В следующем разделе мы опробуем другой подход — авторегрессию.
Авторегрессия
Авторегрессия основана на том принципе, что будущее похоже на прошлое. Например, во временных рядах, которые мы рассматривали до сих пор, присутствует четкий годовой цикл. Поэтому при составлении прогноза на июнь следующего года можно ориентироваться на июнь прошлого.
Чтобы проверить, что этот принцип работает, вернемся к ряду nuclear, содержащему ежемесячный объем выработки электроэнергии на АЭС США, и посчитаем разность значений за один и тот же месяц в следующих один за другим годах, которая называется разностью «год к году» (year-over-year difference):
diff = (nuclear - nuclear.shift(12)).dropna()
diff.plot(label="разница год к году", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Величина этой разности существенно меньше значений исходного ряда, что приводит нас к следующему принципу авторегрессии: возможно, будет проще предсказать эту разность, чем исходные значения.
С этой целью посмотрим, существует ли корреляция между последовательными элементами в ряду разностей. Если да, можно использовать эту корреляцию для предсказания последующих значений на основе предыдущих.
Сначала создадим объект DataFrame, поместив исходные разности в первый столбец и их же со сдвигом на 1, 2 и 3 месяца — в последующие столбцы. Назовем последние lag1, lag2 и lag3 соответственно, поскольку в них содержатся ряды с лагом (lag), или временно́й задержкой:
df_ar = pd.DataFrame({"diff": diff})
for lag in [1, 2, 3]:
df_ar[f"lag{lag}"] = diff.shift(lag)
df_ar = df_ar.dropna()
Вот значения корреляции между столбцами:
df_ar.corr()[["diff"]]
| diff | |
| diff | 1.000000 |
| lag1 | 0.562212 |
| lag2 | 0.292454 |
| lag3 | 0.222228 |
Эти корреляции называются запаздывающими корреляциями (lagged correlations) или автокорреляциями. Префикс «авто» указывает на то, что мы считаем корреляцию ряда с самим собой. Как частный случай, корреляция между diff и lag1 называется последовательной корреляцией (serial correlation), поскольку это корреляция между следующими один за другим элементами в ряду.
Сила этих корреляций достаточна, чтобы предположить, что они могут оказаться полезными для предсказаний. Поэтому включим их в модель множественной регрессии. Функция ниже использует столбцы объекта DataFrame для создания формулы в формате Patsy: первый столбец используется в качестве зависимой переменной, а остальные столбцы — в качестве объясняющих переменных:
def make_formula(df):
"""Составить формулу Patsy из заголовков столбцов."""
y = df.columns[0]
xs = " + ".join(df.columns[1:])
return f"{y} ~ {xs}"
Вот результаты работы линейной модели, предсказывающей следующее значение последовательности на основе трех предыдущих значений:
formula = make_formula(df_ar)
results_ar = smf.ols(formula=formula, data=df_ar).fit()
display_summary(results_ar)
| coef | std err | t | P > | t | | [0.025 | 0.975] | |
| Intercept | 24.2674 | 114.674 | 0.212 | 0.833 | –201.528 | 250.063 |
| lag1 | 0.5847 | 0.061 | 9.528 | 0.000 | 0.464 | 0.706 |
| lag2 | –0.0908 | 0.071 | –1.277 | 0.203 | –0.231 | 0.049 |
| lag3 | 0.1026 | 0.062 | 1.666 | 0.097 | –0.019 | 0.224 |
R-squared: 0.3239
Теперь воспользуемся методом predict для генерации предсказаний значений ряда в прошлом. Вот как выглядит эта ретроспективная оценка в сравнении с фактическими данными:
pred_ar = results_ar.predict(df_ar)
pred_ar.plot(label="предсказанное значение", **pred_options)
diff.plot(label="разность", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Предсказанные значения местами близки к фактическим, но коэффициент R2 всего около 0.319, так что модель оставляет желать лучшего:
resid_ar = (diff - pred_ar).dropna()
R2 = 1 - resid_ar.var() / diff.var()
R2
0.3190252265690783
Один из способов улучшить предсказания — вычислить остатки от предсказаний этой модели и использовать другую модель для предсказания этих остатков, что является еще одним принципом авторегрессии.
Скользящее среднее
Представьте, что сейчас июнь 2019 года и вас просят сделать предсказание на июнь 2020-го. В качестве начального приближения можно принять, что значение этого года повторится в следующем году.
Теперь предположим, что сейчас май 2020 года и вас просят пересмотреть ваш прогноз на июнь 2020-го. Вы могли бы использовать результаты за последние три месяца и автокорреляционную модель из предыдущего раздела, чтобы предсказать разность «год к году».
Наконец, допустим, что вы проверяете предсказания последних нескольких месяцев и видите, что все они занижены. Это означает, что предсказание на следующий месяц также может быть слишком низким, поэтому вы, скорее всего, пересмотрите его в сторону повышения. Основная предпосылка здесь в том, что близкие по времени ошибки предсказаний предсказывают будущие ошибки предсказаний.
Чтобы убедиться в этом, создадим объект DataFrame, где остатки из модели авторегрессии занимают первый столбец, а запаздывающие версии остатков — остальные столбцы. В этом примере я буду использовать лаг в 1 и 6 месяцев:
df_ma = pd.DataFrame({"resid": resid_ar})
for lag in [1, 6]:
df_ma[f"lag{lag}"] = resid_ar.shift(lag)
df_ma = df_ma.dropna()
Воспользуемся классом ols для построения модели авторегрессии остатков. Эта часть модели называется скользящим средним (moving average), поскольку она уменьшает изменчивость предсказаний способом, аналогичным эффекту скользящего среднего. Мне этот термин не кажется информативным, тем не менее он общеупотребим.
Как бы там ни было, вот сводная таблица модели авторегрессии остатков:
formula = make_formula(df_ma)
results_ma = smf.ols(formula=formula, data=df_ma).fit()
display_summary(results_ma)
| coef | std err | t | P > | t | | [0.025 | 0.975] | |
| Intercept | –14.0016 | 114.697 | –0.122 | 0.903 | –239.863 | 211.860 |
| lag1 | 0.0014 | 0.062 | 0.023 | 0.982 | –0.120 | 0.123 |
| lag6 | –0.1592 | 0.063 | –2.547 | 0.011 | –0.282 | –0.036 |
R-squared: 0.0247
Коэффициент R2 довольно мал, так что, похоже, эта часть модели не очень подходит. Но малость p-значения для лага в шесть месяцев означает, что эта переменная содержит больше информации, чем случайный вклад.
Теперь с помощью этой модели можно сгенерировать ретроспективные оценки для остатков:
pred_ma = results_ma.predict(df_ma)
А сейчас, чтобы сгенерировать ретроспективные оценки для разностей год к году, добавим поправку из второй модели к ретроспективным оценкам первой:
pred_diff = pred_ar + pred_ma
Коэффициент R2 для суммы двух моделей, равный примерно 0.332, совсем ненамного лучше результата без поправки скользящего среднего (0.319):
resid_ma = (diff - pred_diff).dropna()
R2 = 1 - resid_ma.var() / diff.var()
R2
0.3315101001391231
Далее мы с помощью этих разностей год к году будем создавать ретроспективные оценки исходных значений.
Ретроспективная оценка с авторегрессией
Начнем генерировать ретроспективные оценки с того, что поместим значения разности год к году в объект Series, в котором индекс соответствует индексу исходных данных:
pred_diff = pd.Series(pred_diff, index=nuclear.index)
Используя метод isna для проверки значений NaN, обнаруживаем, что первый 21 элемент нового объекта Series отсутствует:
n_missing = pred_diff.isna().sum()
n_missing
21
Это объясняется тем, что мы сначала сдвинули временной ряд на 12 месяцев, чтобы рассчитать разность год к году, затем эти разности сдвинули на 3 месяца для первой модели авторегрессии и остатки первой модели — еще на 6 месяцев для второй модели. Каждый раз, когда мы сдвигаем ряд подобным образом, мы теряем несколько значений в начале, а сумма всех этих сдвигов равна 21.
Поэтому, прежде чем мы сможем генерировать ретроспективные оценки, мы должны подготовить наш ряд, скопировав первый 21 элемент из исходного в новый объект Series:
pred_series = pd.Series(index=nuclear.index, dtype=float)
pred_series.iloc[:n_missing] = nuclear.iloc[:n_missing]
Теперь можно исполнить следующий цикл, который заполняет элементы с индекса 21 (22-го элемента) до конца последовательности. Каждый элемент — это сумма значения из предыдущего года и предсказанной разности год к году:
for i in range(n_missing, len(pred_series)):
pred_series.iloc[i] = pred_series.iloc[i - 12] + pred_diff.iloc[i]
Теперь заменим скопированные элементы на NaN, чтобы не ставить себе в заслугу то, что мы идеально «предсказали» первое 21 значение:
pred_series[:n_missing] = np.nan
Вот как выглядит ретроспективная оценка в сравнении с оригинальными данными:
pred_series.plot(метка="предсказанное значение", **pred_options)
nuclear.plot(label="фактическое значение", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Она выглядит довольно неплохо, а коэффициент R2 равен примерно 0.86:
resid = (nuclear - pred_series).dropna()
R2 = 1 - resid.var() / nuclear.var()
R2
0.8586566911201015
Модель, которую мы использовали для расчета этих ретроспективных оценок, называется SARIMA и является одной из моделей семейства ARIMA. Каждая составляющая этих аббревиатур соотносится с определенным элементом модели:
• S означает сезонный (seasonal), поскольку первым шагом было вычисление разностей между значениями, разделенными одним периодом сезонности.
• AR расшифровывается как авторегрессия (autoregression), которую мы использовали для моделирования запаздывающих корреляций в разностях.
• I означает интегрированный (integrated), так как итеративный процесс, который мы использовали для вычисления элементов pred_series, аналогичен интегрированию в математическом анализе.
• MA расшифровывается как скользящее среднее (moving average) — общепринятое название для второй модели авторегрессии, которую мы применяли к остаткам первой.
Модели ARIMA — это мощные и универсальные инструменты для моделирования данных временных рядов.
ARIMA
Модуль StatsModel содержит библиотеку tsa, название которой является сокращением от «анализ временных рядов» (time series analysis). В ней реализована функция ARIMA, которая подбирает модели ARIMA и генерирует прогнозы.
Чтобы подогнать модель SARIMA, которую мы строили в предыдущих разделах, вызовем эту функцию с двумя кортежами — order и seasonal_order — в качестве аргументов. Вот значения в кортеже order, соответствующие модели из предыдущих разделов:
order = ([1, 2, 3], 0, [1, 6])
Значения параметра порядка (order) определяют:
• Какие лаги следует включить в модель авторегрессии (AR): первые три числа в данном примере.
• Сколько раз модель должна считать разность между последовательными элементами. В нашем примере это 0, потому что взамен мы вычисляли сезонную разницу, — мы к этому еще вернемся.
• Какие лаги должны быть включены в модель скользящего среднего (MA): в данном примере это один и шесть.
А вот значения параметра сезонного порядка (seasonal_order):
seasonal_order = (0, 1, 0, 12)
Первый и третий элементы равны 0, так как в этой модели отсутствуют и сезонная авторегрессия (AR), и сезонное скользящее среднее (MA). Второй элемент, равный 1, указывает модели считать сезонные разности, а последний элемент задает период сезонности.
Вот как вызвать функцию ARIMA для создания и подгонки этой модели:
import statsmodels.tsa.api as tsa
model = tsa.ARIMA(nuclear, order=order, seasonal_order=seasonal_order)
results_arima = model.fit()
display_summary(results_arima)
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
| ar.L1 | 0.0458 | 0.379 | 0.121 | 0.904 | –0.697 | 0.788 |
| ar.L2 | –0.0035 | 0.116 | –0.030 | 0.976 | –0.230 | 0.223 |
| ar.L3 | 0.0375 | 0.049 | 0.769 | 0.442 | –0.058 | 0.133 |
| ma.L1 | 0.2154 | 0.382 | 0.564 | 0.573 | –0.533 | 0.964 |
| ma.L6 | –0.0672 | 0.019 | –3.500 | 0.000 | –0.105 | -0.030 |
| sigma2 | 3.473e+06 | 1.9е-07 | 1.83е+13 | 0.000 | 3.47e+06 | 3.47e+06 |
В таблице результатов приведены оценки коэффициентов для трех лагов в модели авторегрессии (AR), двух лагов в модели скользящего среднего (MA) и дисперсия остатков sigma2.
Из объекта results_arima можно извлечь значение поля fittedvalues, в котором содержатся ретроспективные оценки. По той же причине, по которой в начале вычисленного нами ряда отсутствовали значения, в начале fittedvalues есть некорректные значения. Удалим их:
fittedvalues = results_arima.fittedvalues[n_missing:]
Подогнанные (fitted) значения похожи на те, которые мы уже вычислили, но не совпадают в точности — вероятно, из-за того что данная реализация ARIMA по-другому обрабатывает начальные условия:
fittedvalues.plot(label="модель ARIMA", **pred_options)
nuclear.plot(label="фактическое значение", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Коэффициент R2 тоже близок, но не идентичен:
resid = fittedvalues - nuclear
R2 = 1 - resid.var() / nuclear.var()
R2
0.8262717330784233
Функция ARIMA позволяет легко экспериментировать с различными вариантами данной модели.
Предсказание с помощью ARIMA
В объекте, возвращаемом функцией ARIMA, имеется метод get_forecast, который генерирует предсказания. В демонстрационных целях разделим временной ряд на обучающую и тестовую части и применим ту же модель к обучающему ряду:
training, test = split_series(nuclear)
model = tsa.ARIMA(training, order=order, seasonal_order=seasonal_order)
results_training = model.fit()
Полученный результат можно использовать, чтобы создать прогноз для тестового ряда:
forecast = results_training.get_forecast(steps=len(test))
У полученного объекта имеется атрибут forecast_mean, содержащий массив средних значений, и функция, возвращающая доверительный интервал (confidence interval, ci):
forecast_mean = forecast.predicted_mean
forecast_ci = forecast.conf_int()
forecast_ci.columns = ["lower", "upper"]
Чтобы представить результаты графически и сравнить их с фактическими данными временного ряда, выполним следующий код:
plt.fill_between(
forecast_ci.index,
forecast_ci.lower,
forecast_ci.upper,
lw=0,
color="gray",
alpha=0.2,
)
plt.plot(forecast_mean.index, forecast_mean, label="прогноз", **pred_options)
plt.plot(test.index, test, label="фактическое значение", **actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Фактические значения почти полностью укладываются в доверительный интервал прогнозных значений. Вот величина метрики MAPE для предсказаний модели:
MAPE(forecast_mean, test)
3.381754924564627
Предсказания в среднем смещены на 3.38 %, что несколько лучше результатов, полученных в случае сезонной декомпозиции (3.81 %).
Модель ARIMA универсальнее, чем сезонная декомпозиция, и часто позволяет делать более точные предсказания. В данном временном ряду автокорреляция не очень сильна, поэтому выигрыш от использования ARIMA умеренный.
Глоссарий
Временной ряд (time series)
Датасет, в котором каждое значение связано с определенным временем и часто представляет собой измерения, проводимые через регулярные промежутки времени.
Сезонная декомпозиция (seasonal decomposition)
Метод разделения временного ряда на долгосрочный тренд, повторяющуюся сезонную составляющую и остаток.
Обучающий ряд (training series)
Часть временного ряда, используемая для подгонки модели.
Тестовый ряд (test series)
Часть временного ряда, используемая для проверки точности сгенерированных моделью предсказаний.
Ретроспективная оценка (retrodiction)
Предсказание значения, наблюдавшегося в прошлом, часто используемое для тестирования или валидации модели.
Окно (window)
Набор следующих друг за другом значений во временном ряду, используемый для вычисления скользящего среднего.
Скользящее среднее (moving average)
Метод сглаживания колебаний временных рядов, состоящий в усреднении значений в перекрывающихся окнах.
Последовательная корреляция (serial correlation)
Корреляция между следующими друг за другом элементами временного ряда.
Автокорреляция (autocorrelation)
Корреляция между временным рядом и его сдвинутой, или запаздывающей (lagged), версией.
Лаг (lag)
Величина сдвига в последовательной корреляции или автокорреляции.
Упражнения
Упражнение 12.1
В качестве примера сезонной декомпозиции смоделируем среднемесячные значения температуры поверхности в США. Используем датасет проекта «Наш мир в данных» (Our World in Data), который содержит «значения температуры [в градусах Цельсия] воздуха, измеренные на высоте двух метров над уровнем земли, охватывающие поверхности суши, моря, а также внутриматериковые водные поверхности» для большинства стран мира за период с 1950 по 2024 год. Инструкции по загрузке данных приведены в Jupyter-блокноте этой главы.
Данные можно считать так:
temp = pd.read_csv("monthly-average-surface-temperatures-by-year.csv")
В следующей ячейке выбираются данные по США с 2001 года до конца измерений и сохраняются в объекте Series библиотеки Pandas:
temp_us = temp.query("Code == 'USA'")
columns = [str(year) for year in range(2000, 2025)]
temp_series = temp_us.loc[:, columns].transpose().stack()
temp_series.index = pd.date_range(start="2000-01", periods=len(temp_series),
freq="ME")
Вот как выглядят полученные результаты:
temp_series.plot(label="среднемесячное значение", **actual_options)
decorate(ylabel="Температура поверхности (°C)")
Вполне ожидаемо наблюдается четкая сезонная закономерность. Проведите аддитивную сезонную декомпозицию с периодом в 12 месяцев. Подгоните линейную модель к линии тренда. Каково среднегодовое повышение температуры поверхности за этот временной интервал? При желании повторите тот же анализ для других интервалов или данных из других стран.
Упражнение 12.2
Ранее в этой главе мы использовали мультипликативную сезонную декомпозицию для моделирования малой солнечной генерации электроэнергии с 2014 по 2019 год и прогнозирования производства на период с 2019 по 2024 год. Теперь проделаем то же самое для солнечной генерации промышленного масштаба (utility-scale solar power). Вот как выглядит ее временной ряд:
util_solar = elec["United States : all utility-scale solar"].dropna()
util_solar = util_solar[util_solar.index.year >= 2014]
util_solar.plot(**actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
С помощью функции split_series разделите эти данные на обучающий и тестовый ряды. Выполните мультипликативную декомпозицию обучающего ряда за 12-месячный период. Подгоните линейную или квадратичную модель под тренд и сгенерируйте пятилетний прогноз, включая сезонную составляющую. Постройте график прогноза вместе с тестовым рядом и вычислите среднюю абсолютную процентную ошибку (MAPE).
Упражнение 12.3
Посмотрим, насколько хорошо модель ARIMA подходит для моделирования объема электроэнергии, произведенной гидроэлектростанциями в США. Вот как выглядит временной ряд с 2001 по 2024 год:
hydro = elec["United States : conventional hydroelectric"]
hydro.plot(**actual_options)
decorate(ylabel="Выработка (ГВт-ч)")
Подгоните модель SARIMA с периодом сезонности в 12 месяцев к этим данным. Поэкспериментируйте с разными величинами лага в авторегрессии и скользящим средним модели — проверьте, удастся ли вам найти комбинацию, которая максимизирует коэффициент R2 модели. Создайте пятилетний прогноз и постройте его график вместе с доверительным интервалом.
Известна также как автокорреляция первого порядка. — Примеч. пер.
В случае последовательной корреляции лаг всегда равен 1. — Примеч. пер.
