Наука о данных: временные ряды
Прогнозирование озера Эри по временным рядам с использованием модели SARIMA
Линия выбора, сравнения, диагностики и прогнозирования моделей
Краткое содержание:
- Построил модель сезонного интегрированного скользящего среднего с авторегрессией (SARIMA) для отслеживания месячных уровней воды в озере Эри в течение пяти десятилетий.
- Предоставляется конвейер выбора моделей временных рядов с использованием следующих методов: спектральный анализ, дифференцирование, AIC, PACF и различные проверки устойчивости.
- Модель точно предсказывает результаты после перекрестной проверки с реальными данными.
В отличие от других типов статистических моделей, модели временных рядов предполагают зависимые наблюдения и убедительно аргументируют вчерашние влияния сегодня. В этом посте я обучаю модель SARIMA предсказанию месячного уровня воды в озере Эри. Я пройдусь по конвейеру, объясняя, что делать и почему мы это делаем.
- Загрузить, очистить и преобразовать набор данных
#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: SARIMA (1,1,0) * (0,1,1) 12
- Модель 2: SARIMA (1,1,1) * (0,1,1) 12
- Модель 3: SARIMA (3,1,1) * (0,1,1) 12
- Модель 4: SARIMA (1,1,0) * (3,1,1) 12
- Модель 5: SARIMA (1,1,1) * (3,1,1) 12
- Модель 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% -ный доверительный интервал. и точно согласовывать с реальными данными.
Отличный прогноз!
Я надеюсь, что вы, ребята, узнали что-то новое о временных рядах. Этот тип модели чрезвычайно полезен для исследовательских вопросов, связанных со временем и пространством.
Проверить другие мои сообщения.
- Машинное обучение:
- Классификация редких событий
- Как экспрессия генов связана с типом лейкемии с использованием PCA и иерархической кластеризации
2. Наука о данных:
- Оценка программы
- Экспериментирование и причинный вывод: Двусторонняя причинность, Почему экспериментирование?, Ловушки, Причинность vs корреляция
3. О жизни: