Наука о данных: временные ряды

Прогнозирование озера Эри по временным рядам с использованием модели SARIMA

Линия выбора, сравнения, диагностики и прогнозирования моделей

Краткое содержание:

  • Построил модель сезонного интегрированного скользящего среднего с авторегрессией (SARIMA) для отслеживания месячных уровней воды в озере Эри в течение пяти десятилетий.
  • Предоставляется конвейер выбора моделей временных рядов с использованием следующих методов: спектральный анализ, дифференцирование, AIC, PACF и различные проверки устойчивости.
  • Модель точно предсказывает результаты после перекрестной проверки с реальными данными.

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

  1. Загрузить, очистить и преобразовать набор данных
#load the raw dataset
mydata<-read.csv(“monthly-lake-erie-levels-1921–19.csv”)
mydata_new<-mydata[-601,]#remove the last blank row
Lake_Erie_Levels_ts<-ts(mydata_new$Monthly.Lake.Erie.Levels.1921…1970,frequency = 12,start=c(1921,1))#turn into a TS model
# Plot the Dependent Variable
plot(Lake_Erie_Levels_ts,xlab=”Date”,ylab=”Levels”,
 main=”Figure 1: Monthly Lake Erie Levels from 1921 to 1970")

Как правило, мы должны внимательно следить за исходными данными на предмет выбросов, тренда, сезонности и резких скачков.

На рисунке 1 точки данных перемещаются вверх и вниз в определенном диапазоне, что свидетельствует о постоянной дисперсии.

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

Без резких разрывов.

2. Предварительная подготовка: разделение данных и проверка на сезонность и стационарность

Разобьем данные на два набора: обучающий и тестовый.

#We train a model using 90% of the data, aka. 1921.01–1965.12
training = window(Lake_Erie_Levels_ts,start=c(1921,1), end=c(1965,12))
#10% for the test set, aka. 1966.1-1970.12
test = window(Lake_Erie_Levels_ts,start=c(1966,1), end=c(1970,12))#test: 1966.1–1970.12; 10% as the test

Давайте применим спектральный анализ, чтобы определить сезонность.

install.packages(“TSA”)
require(TSA)
periodogram(Lake_Erie_Levels_ts);abline(h=0);axis(1,at=0.084)

Спектральный анализ предполагает период 1 / 0,084 ≈11,9, что довольно близко к 12.

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

library(tseries, quietly = T)
adf.test(training)

	Augmented Dickey-Fuller Test

data:  training
Dickey-Fuller = -2.0622, Lag order = 8, p-value = 0.552
alternative hypothesis: stationary

Чтобы проверить стационарность, мы проводим расширенный тест Дики-Фуллера.

Как оказалось, данные не являются стационарными.

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

training_diff1 <- diff(training,1)
ts.plot(training_diff1,ylab= expression(paste(nable,y)))
var(training_diff1)
adf.test(training_diff1)

Большой! Набор данных теперь стационарный. Среднее значение стабилизировалось с течением времени лишь с несколькими значительными всплесками.

Просто играйте в безопасную карту, давайте проверим, требуется ли второе сравнение.

training_diff2 <- diff(training,2)
var(training_diff2)
adf.test(training_diff2)

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

3. Построение и выбор модели

library(PerformanceAnalytics)
chart.ACFplus(training_diff1, maxlag = 84, elementcolor = “black”)

Чтобы определить потенциальные отношения между лагами, мы строим графики ACF и PACF после первого сравнения.

Эти два графика предполагают сильную сезонную составляющую, потому что АКФ регулярно отмирают. Другими словами, нам нужно включить сезонную составляющую.

training_diff1_seasonal12<-diff(training_diff1,12)
plot.ts(training_diff1_seasonal12)
chart.ACFplus(training_diff1_seasonal12, maxlag = 84, elementcolor = “gray”)

Он становится стационарным только с несколькими всплесками примерно в 1940 году.

На основе моделей ACF и PACF мы предварительно разработали следующие сезонные модели:

1. ACF отключается после задержки 1 и PACF уменьшается → SMA (1), или P = 0, Q = 1.
2. Всплеск ACF на задержке 12 и PACF экспоненциально затухает с сезонными задержками на 12, 24 , 36 и другие → SMA (1), или P = 0 и Q = 1.
3. ACF и PACF оба теряют сезонные задержки → SARMA (3,1), или P = 3, Q = 1

Несезонные модели:

1. PACF отключается после задержки 1 → AR (1) или p = 1, q = 0
2. PACF отключается после задержки 1 и ACF отключается после задержки 1 → p = q = 1
3. PACF отключается после задержки 3 и ACF отключается после задержки 1 → p = 3, q ​​= 1

Всего у нас есть следующие 6 потенциальных моделей с перечисленными сезонными и несезонными компонентами:

  1. Модель 1: SARIMA (1,1,0) * (0,1,1) 12
  2. Модель 2: SARIMA (1,1,1) * (0,1,1) 12
  3. Модель 3: SARIMA (3,1,1) * (0,1,1) 12
  4. Модель 4: SARIMA (1,1,0) * (3,1,1) 12
  5. Модель 5: SARIMA (1,1,1) * (3,1,1) 12
  6. Модель 6: SARIMA (3,1,1) * (3,1,1) 12

4. Сравнение моделей

########### use auto.arima to find the best model automatically
library(forecast)
model_auto= auto.arima(training, stepwise=FALSE, approximation=FALSE)#this function takes a lot of time to run!!!!
summary(model_auto)#ARIMA(1,1,2)(2,0,0)[12] with AIC=670.27 AICc=670.45 BIC=695.31
###### auto selection 
model_auto_1= auto.arima(training)
summary(model_auto_1)# ARIMA(1,1,0)(2,0,0)[12] with AIC=672.23 AICc=672.31 BIC=688.91

Встроенная функция (auto.arima) создает модель в качестве ориентира, сравнивая с выбранной вручную моделью.

model_1 = arima(training, order = c(1,1,0),seasonal = list(order=c(0,1,1),period =12));model_1#AIC=553.74
model_2 = arima(training, order = c(1,1,1),seasonal = list(order=c(0,1,1),period =12));model_2#AIC=555.48
model_3 = arima(training, order = c(3,1,1),seasonal = list(order=c(0,1,1),period =12));model_3#AIC=548.71
model_4 = arima(training, order = c(1,1,0),seasonal = list(order=c(3,1,1),period =12));model_4#AIC=556.94
model_5 = arima(training, order = c(2,1,1),seasonal = list(order=c(3,1,1),period =12));model_6#AIC=558.24
model_6 = arima(training, order = c(3,1,0),seasonal = list(order=c(3,1,1),period =12));model_7#AIC=553.14

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

5. Диагностика модели

ts.plot(residuals(model_3))# a few spikes observed.
# normality check
shapiro.test(residuals(model_3))
qqnorm(residuals(model_3));qqline(residuals(model_3),col =”blue”)#skewed to the right as can be seen from the distribution that a few data points stand above the abline.

Согласно тесту нормальности Шапиро-Уилка и графику Q-Q, остатки не имеют нормального распределения.

chart.ACFplus(residuals(model_3), maxlag = 48, elementcolor = “gray”)

На графиках ACF и PACF большинство остаточных точек попадают в 95% доверительный интервал с несколькими значительными всплесками. Остатки следуют образцу белого шума.

Box.test(residuals(model_3),lag=22,type = c(“Ljung-Box”),fitdf=4)
	Box-Ljung test

data:  residuals(model_3)
X-squared = 23.657, df = 18, p-value = 0.1666

Остатки проходят тест на линейную зависимость после прохождения теста Ljung-Box.

#install.packages(“GeneCycle”)
library(“GeneCycle”)
fisher.g.test(residuals(model_3))
[1] 0.436583

Остатки следуют гауссовскому распределению белого шума после теста Фишера.

6. Прогноз

mypred <- predict(model_3, n.ahead=120)
Upper_tr = mypred$pred + 1.96*mypred$se
Lower_tr = mypred$pred — 1.96*mypred$se
ts.plot(training, xlab=”Year”,ylab=”Level”,main = “Forecast”)
lines(Upper_tr,lty=2,col=”gray”)
lines(Lower_tr,lty=2,col=”gray”)
points(mypred$pred,col=”blue”)
lines(test,col=”red”)#plot test data

Серые пунктирные линии обозначают 95% доверительный интервал, а красные линии - реальные входные данные. Все прогнозируемые точки (синие) попадают в 95% -ный доверительный интервал. и точно согласовывать с реальными данными.

Отличный прогноз!

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

Проверить другие мои сообщения.

  1. Машинное обучение:

2. Наука о данных:

3. О жизни: