27 KiB

Raw Permalink Blame History Unescape Escape

Прогнозиране на времеви редове с ARIMA

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

🎥 Кликнете върху изображението по-горе за видео: Кратко въведение в моделите ARIMA. Примерът е направен в R, но концепциите са универсални.

Тест преди лекцията

Въведение

В този урок ще откриете специфичен начин за изграждане на модели с ARIMA: AutoRegressive Integrated Moving Average. Моделите ARIMA са особено подходящи за данни, които показват нестационарност.

Общи концепции

За да работите с ARIMA, трябва да знаете някои основни концепции:

🎓 Стационарност. В статистически контекст стационарността се отнася до данни, чието разпределение не се променя при изместване във времето. Нестационарните данни показват колебания, дължащи се на тенденции, които трябва да бъдат трансформирани, за да бъдат анализирани. Например сезонността може да въведе колебания в данните и може да бъде премахната чрез процес на "сезонно диференциране".
🎓 Диференциране. Диференцирането на данни, отново в статистически контекст, се отнася до процеса на трансформиране на нестационарни данни, за да станат стационарни чрез премахване на тяхната неконстантна тенденция. "Диференцирането премахва промените в нивото на времевия ред, елиминирайки тенденцията и сезонността и съответно стабилизирайки средната стойност на времевия ред." Статия от Shixiong et al

ARIMA в контекста на времеви редове

Нека разгледаме частите на ARIMA, за да разберем по-добре как помага за моделиране на времеви редове и за правене на прогнози.

AR - Авторегресивен. Авторегресивните модели, както подсказва името, гледат "назад" във времето, за да анализират предишни стойности в данните и да правят предположения за тях. Тези предишни стойности се наричат "лагове". Пример би бил данни, показващи месечни продажби на моливи. Общата стойност на продажбите за всеки месец би се считала за "еволюираща променлива" в набора от данни. Този модел се изгражда като "еволюиращата променлива от интерес се регресира върху собствените си лагови (т.е. предишни) стойности." wikipedia
I - Интегриран. За разлика от подобните модели 'ARMA', 'I' в ARIMA се отнася до неговия интегриран аспект. Данните се "интегрират", когато се прилагат стъпки на диференциране, за да се елиминира нестационарността.
MA - Подвижна средна. Аспектът на подвижната средна в този модел се отнася до изходната променлива, която се определя чрез наблюдение на текущите и миналите стойности на лаговете.

В крайна сметка: ARIMA се използва за създаване на модел, който да пасне възможно най-точно на специалната форма на данни от времеви редове.

Упражнение - изграждане на ARIMA модел

Отворете папката /working в този урок и намерете файла notebook.ipynb.

Стартирайте notebook-а, за да заредите библиотеката statsmodels за Python; ще ви е необходима за ARIMA модели.
Заредете необходимите библиотеки.

Сега заредете още няколко библиотеки, полезни за визуализация на данни:

import os
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import datetime as dt
import math

from pandas.plotting import autocorrelation_plot
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.preprocessing import MinMaxScaler
from common.utils import load_data, mape
from IPython.display import Image

%matplotlib inline
pd.options.display.float_format = '{:,.2f}'.format
np.set_printoptions(precision=2)
warnings.filterwarnings("ignore") # specify to ignore warning messages

Заредете данните от файла /data/energy.csv в Pandas dataframe и ги разгледайте:
```
energy = load_data('./data')[['load']]
energy.head(10)
```
Начертайте всички налични данни за енергия от януари 2012 до декември 2014. Не би трябвало да има изненади, тъй като видяхме тези данни в предишния урок:
```
energy.plot(y='load', subplots=True, figsize=(15, 8), fontsize=12)
plt.xlabel('timestamp', fontsize=12)
plt.ylabel('load', fontsize=12)
plt.show()
```
Сега, нека изградим модел!

Създаване на тренировъчни и тестови набори от данни

Сега данните ви са заредени, така че можете да ги разделите на тренировъчен и тестов набор. Ще обучите модела си върху тренировъчния набор. Както обикновено, след като моделът приключи обучението, ще оцените неговата точност, използвайки тестовия набор. Трябва да се уверите, че тестовият набор обхваща по-късен период от време спрямо тренировъчния набор, за да гарантирате, че моделът не получава информация от бъдещи времеви периоди.

Отделете двумесечен период от 1 септември до 31 октомври 2014 за тренировъчния набор. Тестовият набор ще включва двумесечния период от 1 ноември до 31 декември 2014:
```
train_start_dt = '2014-11-01 00:00:00'
test_start_dt = '2014-12-30 00:00:00'
```
Тъй като тези данни отразяват дневната консумация на енергия, има силен сезонен модел, но консумацията е най-сходна с консумацията през по-скорошни дни.
Визуализирайте разликите:
```
energy[(energy.index < test_start_dt) & (energy.index >= train_start_dt)][['load']].rename(columns={'load':'train'}) \
    .join(energy[test_start_dt:][['load']].rename(columns={'load':'test'}), how='outer') \
    .plot(y=['train', 'test'], figsize=(15, 8), fontsize=12)
plt.xlabel('timestamp', fontsize=12)
plt.ylabel('load', fontsize=12)
plt.show()
```
Следователно, използването на сравнително малък времеви прозорец за обучение на данните би трябвало да е достатъчно.

Забележка: Тъй като функцията, която използваме за настройка на ARIMA модела, използва вътрешна валидация по време на настройката, ще пропуснем валидационните данни.

Подготовка на данните за обучение

Сега трябва да подготвите данните за обучение, като извършите филтриране и скалиране на данните. Филтрирайте набора от данни, за да включва само необходимите времеви периоди и колони, и скалирайте, за да гарантирате, че данните са проектирани в интервала 0,1.

Филтрирайте оригиналния набор от данни, за да включва само споменатите времеви периоди за всеки набор и само необходимата колона 'load' плюс датата:
```
train = energy.copy()[(energy.index >= train_start_dt) & (energy.index < test_start_dt)][['load']]
test = energy.copy()[energy.index >= test_start_dt][['load']]

print('Training data shape: ', train.shape)
print('Test data shape: ', test.shape)
```
Можете да видите формата на данните:
```
Training data shape:  (1416, 1)
Test data shape:  (48, 1)
```

Скалирайте данните, за да бъдат в диапазона (0, 1).

scaler = MinMaxScaler()
train['load'] = scaler.fit_transform(train)
train.head(10)

Визуализирайте оригиналните спрямо скалираните данни:

energy[(energy.index >= train_start_dt) & (energy.index < test_start_dt)][['load']].rename(columns={'load':'original load'}).plot.hist(bins=100, fontsize=12)
train.rename(columns={'load':'scaled load'}).plot.hist(bins=100, fontsize=12)
plt.show()

Оригиналните данни

Скалираните данни

Сега, след като сте калибрирали скалираните данни, можете да скалирате тестовите данни:
```
test['load'] = scaler.transform(test)
test.head()
```

Имплементиране на ARIMA

Време е да имплементирате ARIMA! Сега ще използвате библиотеката statsmodels, която инсталирахте по-рано.

Сега трябва да следвате няколко стъпки:

Дефинирайте модела, като извикате SARIMAX() и подадете параметрите на модела: параметрите p, d и q, както и параметрите P, D и Q.
Подгответе модела за тренировъчните данни, като извикате функцията fit().
Направете прогнози, като извикате функцията forecast() и зададете броя на стъпките (хоризонта), които да прогнозирате.

🎓 Какво означават всички тези параметри? В ARIMA модел има 3 параметъра, които се използват за моделиране на основните аспекти на времеви редове: сезонност, тенденция и шум. Тези параметри са:

p: параметър, свързан с авторегресивния аспект на модела, който включва минали стойности. d: параметър, свързан с интегрирания аспект на модела, който влияе върху количеството диференциране (🎓 помните ли диференцирането 👆?) за прилагане към времевия ред. q: параметър, свързан с аспекта на подвижната средна в модела.

Забележка: Ако вашите данни имат сезонен аспект - какъвто е случаят тук - използваме сезонен ARIMA модел (SARIMA). В този случай трябва да използвате друг набор от параметри: P, D и Q, които описват същите асоциации като p, d и q, но съответстват на сезонните компоненти на модела.

Започнете, като зададете предпочитаната стойност за хоризонта. Нека опитаме 3 часа:
```
# Specify the number of steps to forecast ahead
HORIZON = 3
print('Forecasting horizon:', HORIZON, 'hours')
```
Изборът на най-добрите стойности за параметрите на ARIMA модела може да бъде предизвикателство, тъй като е донякъде субективен и отнема време. Може да обмислите използването на функцията auto_arima() от библиотеката pyramid.

Засега опитайте някои ръчни селекции, за да намерите добър модел.

order = (4, 1, 0)
seasonal_order = (1, 1, 0, 24)

model = SARIMAX(endog=train, order=order, seasonal_order=seasonal_order)
results = model.fit()

print(results.summary())

Таблица с резултати се отпечатва.

Създадохте първия си модел! Сега трябва да намерим начин да го оценим.

Оценка на модела

За да оцените модела си, можете да извършите така наречената walk forward валидация. На практика, моделите за времеви редове се преобучават всеки път, когато се появят нови данни. Това позволява на модела да направи най-добрата прогноза на всяка стъпка.

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

Забележка: Трябва да запазите прозореца на тренировъчния набор фиксиран за по-ефективно обучение, така че всеки път, когато добавите ново наблюдение към тренировъчния набор, да премахнете наблюдението от началото на набора.

Този процес предоставя по-устойчива оценка за това как моделът ще се представи на практика. Въпреки това, той идва с изчислителната цена на създаването на толкова много модели. Това е приемливо, ако данните са малки или ако моделът е прост, но може да бъде проблем при мащаб.

Walk-forward валидацията е златният стандарт за оценка на модели за времеви редове и се препоръчва за вашите собствени проекти.

Първо, създайте тестова точка за данни за всяка стъпка на HORIZON.

test_shifted = test.copy()

for t in range(1, HORIZON+1):
    test_shifted['load+'+str(t)] = test_shifted['load'].shift(-t, freq='H')

test_shifted = test_shifted.dropna(how='any')
test_shifted.head(5)

		load	load+1	load+2
2014-12-30	00:00:00	0.33	0.29	0.27
2014-12-30	01:00:00	0.29	0.27	0.27
2014-12-30	02:00:00	0.27	0.27	0.30
2014-12-30	03:00:00	0.27	0.30	0.41
2014-12-30	04:00:00	0.30	0.41	0.57

Данните се изместват хоризонтално според точката на хоризонта.

Направете прогнози за тестовите данни, използвайки този подход с плъзгащ прозорец в цикъл с дължина на тестовите данни:

%%time
training_window = 720 # dedicate 30 days (720 hours) for training

train_ts = train['load']
test_ts = test_shifted

history = [x for x in train_ts]
history = history[(-training_window):]

predictions = list()

order = (2, 1, 0)
seasonal_order = (1, 1, 0, 24)

for t in range(test_ts.shape[0]):
    model = SARIMAX(endog=history, order=order, seasonal_order=seasonal_order)
    model_fit = model.fit()
    yhat = model_fit.forecast(steps = HORIZON)
    predictions.append(yhat)
    obs = list(test_ts.iloc[t])
    # move the training window
    history.append(obs[0])
    history.pop(0)
    print(test_ts.index[t])
    print(t+1, ': predicted =', yhat, 'expected =', obs)

Можете да наблюдавате как се извършва обучението:

2014-12-30 00:00:00
1 : predicted = [0.32 0.29 0.28] expected = [0.32945389435989236, 0.2900626678603402, 0.2739480752014323]

2014-12-30 01:00:00
2 : predicted = [0.3  0.29 0.3 ] expected = [0.2900626678603402, 0.2739480752014323, 0.26812891674127126]

2014-12-30 02:00:00
3 : predicted = [0.27 0.28 0.32] expected = [0.2739480752014323, 0.26812891674127126, 0.3025962399283795]

Сравнете прогнозите с действителното натоварване:

eval_df = pd.DataFrame(predictions, columns=['t+'+str(t) for t in range(1, HORIZON+1)])
eval_df['timestamp'] = test.index[0:len(test.index)-HORIZON+1]
eval_df = pd.melt(eval_df, id_vars='timestamp', value_name='prediction', var_name='h')
eval_df['actual'] = np.array(np.transpose(test_ts)).ravel()
eval_df[['prediction', 'actual']] = scaler.inverse_transform(eval_df[['prediction', 'actual']])
eval_df.head()

Резултат

		timestamp	h	prediction	actual
0	2014-12-30	00:00:00	t+1	3,008.74	3,023.00
1	2014-12-30	01:00:00	t+1	2,955.53	2,935.00
2	2014-12-30	02:00:00	t+1	2,900.17	2,899.00
3	2014-12-30	03:00:00	t+1	2,917.69	2,886.00
4	2014-12-30	04:00:00	t+1	2,946.99	2,963.00

Наблюдавайте прогнозата за почасовите данни, сравнена с действителното натоварване. Колко точна е тя?

Проверка на точността на модела

Проверете точността на модела си, като тествате неговата средна абсолютна процентна грешка (MAPE) за всички прогнози.

🧮 Покажи ми математиката

MAPE се използва за показване на точността на прогнозата като съотношение, дефинирано от горната формула. Разликата между действителната и прогнозната стойност се дели на действителната.

"Абсолютната стойност в това изчисление се сумира за всяка прогнозирана точка във времето и се дели на броя на точките n." wikipedia

Изразете уравнението в код:

if(HORIZON > 1):
    eval_df['APE'] = (eval_df['prediction'] - eval_df['actual']).abs() / eval_df['actual']
    print(eval_df.groupby('h')['APE'].mean())

Изчислете MAPE за една стъпка:

print('One step forecast MAPE: ', (mape(eval_df[eval_df['h'] == 't+1']['prediction'], eval_df[eval_df['h'] == 't+1']['actual']))*100, '%')

MAPE за прогнозата на една стъпка: 0.5570581332313952 %

Отпечатайте MAPE за многократна прогноза:
```
print('Multi-step forecast MAPE: ', mape(eval_df['prediction'], eval_df['actual'])*100, '%')
```
```
Multi-step forecast MAPE:  1.1460048657704118 %
```
Най-добре е да имате ниска стойност: имайте предвид, че прогноза с MAPE от 10 е с отклонение от 10%.

Но както винаги, е по-лесно да видите този вид измерване на точността визуално, затова нека го начертаем:

 if(HORIZON == 1):
    ## Plotting single step forecast
    eval_df.plot(x='timestamp', y=['actual', 'prediction'], style=['r', 'b'], figsize=(15, 8))

else:
    ## Plotting multi step forecast
    plot_df = eval_df[(eval_df.h=='t+1')][['timestamp', 'actual']]
    for t in range(1, HORIZON+1):
        plot_df['t+'+str(t)] = eval_df[(eval_df.h=='t+'+str(t))]['prediction'].values

    fig = plt.figure(figsize=(15, 8))
    ax = plt.plot(plot_df['timestamp'], plot_df['actual'], color='red', linewidth=4.0)
    ax = fig.add_subplot(111)
    for t in range(1, HORIZON+1):
        x = plot_df['timestamp'][(t-1):]
        y = plot_df['t+'+str(t)][0:len(x)]
        ax.plot(x, y, color='blue', linewidth=4*math.pow(.9,t), alpha=math.pow(0.8,t))

    ax.legend(loc='best')

plt.xlabel('timestamp', fontsize=12)
plt.ylabel('load', fontsize=12)
plt.show()

🏆 Много хубав график, показващ модел с добра точност. Браво!

🚀Предизвикателство

Разгледайте начините за тестване на точността на модел за времеви серии. В този урок разглеждаме MAPE, но има ли други методи, които можете да използвате? Проучете ги и ги опишете. Полезен документ може да бъде намерен тук

Тест след лекцията

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

Този урок обхваща само основите на прогнозиране на времеви серии с ARIMA. Отделете време да задълбочите знанията си, като разгледате този репозиторий и различните му типове модели, за да научите други начини за изграждане на модели за времеви серии.

Задача

Нов ARIMA модел

Отказ от отговорност:
Този документ е преведен с помощта на AI услуга за превод Co-op Translator. Въпреки че се стремим към точност, моля, имайте предвид, че автоматизираните преводи може да съдържат грешки или неточности. Оригиналният документ на неговия роден език трябва да се счита за авторитетен източник. За критична информация се препоръчва професионален човешки превод. Ние не носим отговорност за недоразумения или погрешни интерпретации, произтичащи от използването на този превод.

27 KiB Raw Permalink Blame History Unescape Escape