На динамику цен на сырую нефть влияют различные факторы. Некоторые из этих факторов являются политическими, а другие связаны с инфляцией, экономикой и т. Д. Можно ли предсказать будущие цены на нефть на основе исторических данных? Эта статья даст ответ на этот вопрос.
Во время обзора литературы по этой теме мне удалось найти предыдущую работу Chen et al. использование глубокого обучения для прогнозирования сырой нефти WTI на основе долгосрочной краткосрочной памяти LSTM.
Временные ряды
Временной ряд - это последовательность записей, хранящихся через регулярные промежутки времени.
В зависимости от частоты интервала временные ряды могут записывать годовой, квартальный, ежемесячный, ежедневный, ежечасный и т. Д. Например, еженедельные продажи в магазинах, ежедневный прогноз погоды, ежемесячный доход, почасовые цены на акции, трафик в минутах или секундах и т. Д. .
Мы используем временные ряды для прогнозирования будущего, например, мы можем прогнозировать продажи в магазине на следующий месяц на основе истории предыдущих лет. Мы можем прогнозировать инфляцию, основываясь на росте цен за последние несколько лет. Прогноз временных рядов может помочь в прогнозировании будущих продаж автомобилей, чтобы компании могли улучшить свои продуктовые линейки, цепочки поставок и т. Д.
Набор данных
Набор данных доступен с полным примером здесь. Источником набора данных является Управление энергетической информации США. CSV-файл содержит два столбца: Дата и Цена. Данные содержат ежедневные исторические цены на нефть марки Brent с 17 мая 1987 года по 30 сентября 2019 года. Давайте исследуем и настроим набор данных:
#import the csv file
oilPrices = pd.read_csv('/kaggle/input/brent-oil-prices/BrentOilPrices.csv')
#change column names to more comfortable names
oilPrices.columns=['date', 'price']
#Cast Date Column to type date
oilPrices['date'] = pd.to_datetime(oilPrices['date'])
print("Data Set:"% oilPrices.columns, oilPrices.shape)
print("Data Types:", oilPrices.dtypes)
#Check the top five records
oilPrices.head()
Как вы могли заметить, данные временных рядов не содержат значений для субботы и воскресенья, поскольку рынок закрыт по выходным. Следовательно, необходимо заполнить недостающие значения. Чтобы заполнить выходные дни, сначала сделайте дату как индекс (для метода повторной выборки), затем используйте форвардную заливку ffill (), которая присвоит значениям выходных значение пятницы. Метод повторной выборки используется для преобразования частоты и повторной выборки временных рядов. Объект должен иметь индекс, подобный datetime (DatetimeIndex, PeriodIndex или TimedeltaIndex), или передавать значения, подобные datetime, в ключевое слово on или level.
oilPrices.set_index('date', inplace=True)
oilPrices = oilPrices.resample('D').ffill().reset_index()
Убедитесь, что у нас нет нулевых значений:
oilPrices.isnull().values.any()
Давайте разделим дату на год, месяц и неделю, чтобы изучить тенденции цен на нефть:
oilPrices['year']=oilPrices['date'].dt.year oilPrices['month']=oilPrices['date'].dt.month oilPrices['week']=oilPrices['date'].dt.week
Разделение на обучение и тестирование для прогнозирования цен на нефть на 2019 год. Обратите внимание, что я обучал временные ряды с 2000 года, поскольку они продемонстрировали более высокую точность:
train = oilPrices[(oilPrices['date' ] > '2000-01-01') & (oilPrices['date' ] <= '2018-12-31')] test = oilPrices[oilPrices['date' ] >= '2019-01-01']
Годовая визуализация цен:
yearlyPrice=train.groupby(["year"])['price'].mean()
plt.figure(figsize=(16,4))
plt.title('Oil Prices')
plt.xlabel('Year')
plt.ylabel('Price')
yearlyPrice.plot()
plt.show();
Модель ARIMA
ARIMA используется для моделирования несезонных временных рядов, это означает Авторегрессивная интегрированная скользящая средняя. Я буду использовать ARIMA в качестве базовой модели, чтобы протестировать ее с моделью Prophet, чтобы увидеть, какая модель дает лучший прогноз. Полная ссылка доступна здесь.
Сначала нам нужно подготовить данные для оценщика ARIMA:
#Convert to Time Series For ARIMA Estimator series=pd.Series(data=train['price'].to_numpy(), index=train['date'])
Для прогнозирования временных рядов с помощью ARIMA необходимо установить значения трех параметров (p, d, q):
- p: порядок авторегрессивной (AR) модели (т. е. количество наблюдений с запаздыванием).
- d: степень различия.
- q: порядок модели скользящей средней (MA). По сути, это размер функции «окна» над данными временного ряда.
Найдите порядок разности (d) в модели ARIMA
Цель дифференцирования - сделать временной ряд стационарным. Для этого нам нужно измерить разницу между сегодня и вчера, а затем (завтра-вчера и т. Д.), Пока мы не дойдем до стационарных данных.
Стационарный временной ряд - это тот, в котором статистические свойства, такие как среднее значение и дисперсия, являются постоянными во времени.
Нам нужны стационарные данные для прогнозирования временных рядов. Правильный порядок разности - это минимальная разность, необходимая для получения почти стационарного ряда, который перемещается вокруг определенного среднего значения, а график функции автокорреляции (ACF) довольно быстро достигает нуля.
Чтобы проверить, есть ли у нас стационарные временные ряды, можно использовать расширенный метод Дики-Фуллера, чтобы проверить единичный корень в одномерном процессе при наличии последовательной корреляции. Если p-значение i> 0,05, которое и было истинным: 0,286247, мы продолжим поиск порядка разности.
from statsmodels.tsa.stattools import adfuller
from numpy import log
result = adfuller(series.dropna())
print('ADF Statistic: %f' % result[0])
print('p-value: %f' % result[1])
Но перед этим давайте нарисуем первые 100 записей во временном ряду, чтобы увидеть, можно ли виртуально продемонстрировать p-значение:
plt.plot(series[0:100]) plt.show()
Поскольку временной ряд не является стационарным, нам нужно провести разность, если бы временной ряд был стационарным, значение d было бы 0. Давайте сделаем разницу за один день, чтобы увидеть, сделает ли это данные стационарными:
daily_series_diff1 = series.diff(periods=1).dropna()
daily_series_diff2 = daily_series_diff1.diff(periods=1).dropna()
fig, ax = plt.subplots()
ax.plot(daily_series_diff1[0:100], label='1st Order Differencing')
ax.plot(daily_series_diff2[0:100], label='2nd Order Differencing')
plt.ylim([-3,3])
legend = ax.legend(loc='upper center', shadow=True)
plt.title('Time Series')
plt.xlabel('Date')
plt.ylabel('Diff')
plt.show()
Образец, который мы видим выше, предназначен только для первых 100 записей временного ряда, если мы посмотрим на исходный временной ряд за 18 лет, мы найдем следующее:
В целом временной ряд не кажется стационарным, но во время тестирования лучшее значение для d было равно нулю.
Временной ряд имеет стационарность, если сдвиг во времени не приводит к изменению формы распределения; единичные корни являются одной из причин нестационарности, здесь мы замечаем, что после первого дифференцирования форма распределения не меняется.
Результат показывает, что наилучшее значение для d равно 1, но давайте сделаем заключительный тест, выполнив:
- Расширенный тест Дики-Фуллера (ADF): это тест единичного корня на стационарность. Корни единиц могут привести к непредсказуемым результатам в вашем анализе временных рядов.
- Квятковски – Филлипс – Шмидт – Шин (KPSS): тест выясняет, является ли временной ряд стационарным вокруг среднего или линейного тренда или является нестационарным из-за в единичный корень.
- Тест Филлипса – Перрона (PP) : это тест единичного корня, это модификация теста Дики Фуллера и исправляет автокорреляцию и гетероскедастичность в ошибках
#Number of differences required for a stationary series
from pmdarima.arima.utils import ndiffs
y=series
# augmented Dickey–Fuller test (adf test)
print("ADF Test: ",ndiffs(y, test='adf'))
# Kwiatkowski–Phillips–Schmidt–Shin (KPSS) test
print("KPSS Test: ",ndiffs(y, test='kpss'))
# Phillips–Perron (PP) test:
print("PP Test: ",ndiffs(y, test='pp'))
По результатам трех тестов был выставлен 1, это было сделано только для того, чтобы убедиться, что решение о выборе порядка разности значений d = 1 было наиболее правильным.
Порядок авторегрессии Срок P
Следующий шаг - определить, требует ли модель условий дополненной реальности. Чтобы узнать необходимое количество членов AR, нам нужно проверить график частичной автокорреляции (PACF):
plt.rcParams.update({'figure.figsize':(12,3), 'figure.dpi':120})
from statsmodels.graphics.tsaplots import plot_pacf
fig, axes = plt.subplots(1, 2, sharex=True)
plot_pacf(daily_series_diff1, lags=10, ax=axes[0], title="Partial Autocorrelation 1st Order Differencing")
plot_pacf(daily_series_diff2, lags=10, ax=axes[1], title="Partial Autocorrelation 2nd Order Differencing")
plt.xlabel('Lag')
plt.ylabel('PACF')
plt.show()
PACF - это корреляция между рядом и его запаздыванием после исключения вкладов от промежуточных запаздываний. Итак, PACF как бы передает чистую корреляцию между запаздыванием и серией. Таким образом, мы узнаем, нужна ли эта задержка в термине AR или нет. Если мы рассмотрим PACF после разности 1-го порядка p = 0 или разности 2-го порядка, то p = 1.
Определение порядка члена скользящей средней q
Точно так же, как мы смотрели на график PACF для количества членов AR, мы можем посмотреть на график ACF для количества членов MA. Срок скользящей средней - это технически ошибка запаздывающего прогноза.
ACF сообщает, сколько членов MA необходимо для устранения любой автокорреляции в стационарном ряду.
Если автокорреляция положительна для большого количества лагов (10 или более), то ряд нуждается в дальнейшем дифференцировании. С другой стороны, если автокорреляция с запаздыванием 1 слишком отрицательна, то ряды, вероятно, слишком сильно различаются. Я сделал два различия, чтобы увидеть автокорреляцию ниже:
plt.rcParams.update({'figure.figsize':(12,3), 'figure.dpi':120})
from statsmodels.graphics.tsaplots import plot_acf
fig, axes = plt.subplots(1, 2, sharex=True)
plot_acf(daily_series_diff1, lags=20, ax=axes[0], title="Autocorrelation 1st Order Differencing")
plot_acf(daily_series_diff2, lags=20, ax=axes[1], title="Autocorrelation 1st Order Differencing")
plt.xlabel('Lag')
plt.ylabel('ACF')
plt.show()
Для 1-го разложения q = 0. Только одно запаздывание значительно выше линии значимости при втором различении. Итак, давайте установим q на 1. Это было решено после тестирования. Итак, три параметра (p, d, q) равны (1,0,1)
Мы также можем случайным образом подобрать два значения параметров для модели, чтобы найти наиболее подходящее значение, путем измерения информационного критерия Акаике (AIC), который оценивает качество каждой модели:
import pmdarima as pm
model = pm.auto_arima(series, start_p=1, start_q=1,
test='adf', # use adftest to find optimal 'd'
max_p=3, max_q=3, # maximum p and q
m=1, # frequency of series
d=None, # let model determine 'd'
seasonal=False, # No Seasonality
start_P=0,
D=0,
trace=True,
error_action='ignore',
suppress_warnings=True,
stepwise=True)
print(model.summary())
Итак, окончательная форма модели такова:
from statsmodels.tsa.arima_model import ARIMA model = ARIMA(series, order=(1, 0, 1)).fit(transparams=False) print(model.summary())
Заметим, что P ›| z | разумно меньше 0,05. Теперь давайте спрогнозируем цены на нефть марки Brent на период начало = '1/1/2019', конец = '30.09.2019'.
ARIMA_Predict = model.predict(start='1/1/2019', end='9/30/2019')
Модель пророка
Prophet - это инструмент, разработанный Facebook для создания высококачественных прогнозов для данных временных рядов, которые имеют множественную сезонность с линейным или нелинейным ростом.
Краткое описание того, как работает Prophet
По сути, Prophet представляет собой аддитивную модель со следующими компонентами:
y(t) = g(t) + s(t) + h(t) + ϵₜ
- g (t) моделирует тенденцию, которая описывает долгосрочное увеличение или уменьшение данных. Prophet включает в себя две модели тенденций, модель насыщающего роста и кусочно-линейную модель, в зависимости от типа задачи прогнозирования.
- s (t) моделирует сезонность с помощью ряда Фурье, который описывает, как на данные влияют сезонные факторы, такие как время года (например, больше запросов на гоголь-моголь во время зимних праздников. )
- h (t) моделирует влияние праздников или крупных событий, которые влияют на временные ряды бизнеса (например, запуск нового продукта, Черная пятница, Суперкубок и т. д.)
- ϵₜ представляет собой термин неприводимой ошибки
Полная информация о Пророке доступна в статье Масштабный прогноз. Теперь давайте импортируем и обучим модель пророка:
from fbprophet import Prophet
d={'ds':train['date'],'y':train['price']}
df_pred=pd.DataFrame(data=d)
model = Prophet(daily_seasonality=False)
model.fit(df_pred)
Теперь спрогнозируем будущие цены на нефть с 1 января 2019 года по 30 сентября 2019 года, 273 дня.
future = model.make_future_dataframe(periods=273)
forecast = model.predict(future)
plt.figure(figsize=(18, 6))
model.plot(forecast, xlabel = 'Date', ylabel = 'Price')
plt.title('Brent Oil Price Prediction');
Возьмем только образец 2019 года:
forecast2019 = forecast[(forecast['ds' ] >= '2019-01-01') & (forecast['ds' ] <= '2019-09-30')]
Оценка
Поскольку мы использовали модели ARIMA и Prophet для прогнозирования будущих цен на нефть, позвольте нам изучить как прогнозируемые, так и фактические цены на обеих моделях, а затем оценить результаты:
fig, ax = plt.subplots()
ax.plot(forecast2019['ds'], ARIMA_Predict, label='Predicted Prices')
ax.plot(test['date'], test['price'], label='Original Prices')
plt.ylim([0,100])
legend = ax.legend(loc='upper center', shadow=True)
plt.title('ARIMA Model Brent Oil Prices Forecast 2019')
plt.xlabel('Month')
plt.ylabel('Price')
plt.show()
fig, ax = plt.subplots()
ax.plot(forecast2019['ds'], forecast2019['yhat'], label='Predicted Prices')
ax.plot(test['date'], test['price'], label='Original Prices')
plt.ylim([0,100])
legend = ax.legend(loc='upper center', shadow=True)
plt.title('Prophet Model Brent Oil Prices Forecast 2019')
plt.xlabel('Month')
plt.ylabel('Price')
plt.show()
Первоначальное наблюдение показывает, что прогнозируемые цены с помощью Пророка немного выше, чем фактические цены, в то время как прогнозируемые цены в ARIMA ниже, но прогноз обеих моделей насколько близок к фактическим ценам, давайте углубимся в меры ошибок.
Средняя абсолютная ошибка (MAE): MAE измеряет среднюю величину ошибок в наборе прогнозов без учета их направления. MAE - это среднее расстояние по вертикали и горизонтали между каждой точкой и идентификационной линией. Абсолютные различия между предсказанными и наблюдаемыми значениями суммируются, а затем вычисляется среднее значение, при котором все индивидуальные различия имеют равный вес.
from sklearn.metrics import mean_absolute_error
maeARIMA=mean_absolute_error(test['price'],ARIMA_Predict)
maeProphet=mean_absolute_error(test['price'],forecast2019['yhat'])
print('Mean Absolute Error ARIMA = {}'.format(round(maeARIMA, 2)))
print('Mean Absolute Error Prophet = {}'.format(round(maeProphet, 2)))
Средняя абсолютная ошибка для ARIMA составляет [12.23]; в то время как Средняя абсолютная ошибка для Пророка составляет [7.24]
Среднеквадратичная ошибка (MSE): MSE - это правило квадратичной оценки, которое также измеряет среднюю величину ошибки. Это квадратный корень из среднего квадрата разницы между прогнозируемыми и наблюдаемыми значениями.
mseProphet = mean_squared_error(test['price'],forecast2019['yhat'])
mseARIMA = mean_squared_error(test['price'],ARIMA_Predict)
print('The Mean Squared Error of ARIMA forecasts is {}'.format(round(mseARIMA, 2)))
print('The Mean Squared Error of Prophet forecasts is {}'.format(round(mseProphet, 2)))
Среднеквадратичная ошибка прогнозов ARIMA составляет [172,32], тогда как для модели Пророка она была [72,04]. Мое наблюдение за разницей между MAE и MSE заключается в том, что разница между предсказанными и наблюдаемыми значениями возводится в квадрат перед усреднением, поэтому, если абсолютное значение разницы между предсказанными и наблюдаемыми больше единицы, MSE покажет более высокую ошибку, чем MAE, в то время как, если абсолютное значение разницы меньше единицы, MSE будет меньше MAE. MSE - это функция риска, соответствующая ожидаемому значению квадрата потери ошибок. Я использовал оба метода в своей оценке для демонстрации.
Среднеквадратичная ошибка (RMSE): мера ошибки прогнозирования, которая позволяет получить более реалистичный корень MSE, поскольку измеряет стандартное отклонение остатков (остаток - это расстояние по вертикали. между точкой данных и линией регрессии), более подробную информацию можно найти в этой статье.
rmseProphet = sqrt(mseProphet)
rmseARIMA = sqrt(mseARIMA)
print('The Root Mean Squared Error of ARIMA forecasts is {}'.format(round(rmseARIMA, 2)))
print('The Root Mean Squared Error of Prophet forecasts is {}'.format(round(rmseProphet, 2)))
Среднеквадратичная ошибка прогнозов ARIMA составляет [13,13], а для Пророка - [8,49]
Выводы
Статья направлена на исследование использования прогнозов временных рядов для прогнозирования будущих цен на нефть. Прогнозируемые цены для обеих моделей ARIMA и Prophet не далеки от фактических цен, несмотря на колебания цен на нефть из-за многих факторов. Модель Prophet в этом эксперименте продемонстрировала немного более высокую точность, чем модель ARIMA.
