Прогноз индекса Московской биржи - ARIMA

Загрузка необходимых пакетов:

library(tidyverse)
library(readxl)
library(ggplot2)
library(forecast)
library(tseries)

Загрузка набора данных:

stock_data <- read_excel("stockk1.xlsx")
names(stock_data) <- c("IMEX", "USD", "SP", "Oil", "P", "IQ")

Преобразование фрейма данных в вид фрейма временных рядов:

stock_data_ts<-ts(stock_data,start = c(2000, 3), frequency = 12)
class(stock_data_ts)

[1] "mts"    "ts"     "matrix"

head(stock_data_ts)

           IMEX   USD      SP   Oil      P    IQ
Mar 2000 198.99 28.65 1366.42 29.04 101.04 110.8
Apr 2000 257.30 28.60 1498.58 24.15 100.64 109.0
May 2000 244.88 28.40 1452.43 23.33 100.89 112.0
Jun 2000 221.10 28.23 1420.60 29.03 101.75 110.8
Jul 2000 199.13 28.05 1454.60 31.30 102.55 110.9
Aug 2000 209.47 27.82 1430.83 25.43 101.79 113.5

plot(stock_data_ts)

Извлечение временного ряда Московской биржи (IMEX):

IMEX <- stock_data_ts[,1]
#This is the start of the time series
start(IMEX)

[1] 2000    3

#This is the end of the time series
end(IMEX)

[1] 2018    3

summary(IMEX)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  135.7   570.8  1386.9  1179.8  1663.1  2292.2 

plot(IMEX)
#This will fit in a line
abline(reg=lm(IMEX~time(IMEX)))

#This will print the cycle across years.
cycle(IMEX)

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2000           3   4   5   6   7   8   9  10  11  12
2001   1   2   3   4   5   6   7   8   9  10  11  12
2002   1   2   3   4   5   6   7   8   9  10  11  12
2003   1   2   3   4   5   6   7   8   9  10  11  12
2004   1   2   3   4   5   6   7   8   9  10  11  12
2005   1   2   3   4   5   6   7   8   9  10  11  12
2006   1   2   3   4   5   6   7   8   9  10  11  12
2007   1   2   3   4   5   6   7   8   9  10  11  12
2008   1   2   3   4   5   6   7   8   9  10  11  12
2009   1   2   3   4   5   6   7   8   9  10  11  12
2010   1   2   3   4   5   6   7   8   9  10  11  12
2011   1   2   3   4   5   6   7   8   9  10  11  12
2012   1   2   3   4   5   6   7   8   9  10  11  12
2013   1   2   3   4   5   6   7   8   9  10  11  12
2014   1   2   3   4   5   6   7   8   9  10  11  12
2015   1   2   3   4   5   6   7   8   9  10  11  12
2016   1   2   3   4   5   6   7   8   9  10  11  12
2017   1   2   3   4   5   6   7   8   9  10  11  12
2018   1   2   3                                    

#Box plot across months will give us a sense on seasonal effect
boxplot(IMEX~cycle(IMEX))

Вывод:

• Тренд показывает, что IMEXs увеличиваются год за годом => временой ряд IMEX нестацонарный
• среднее значение в каждом месяце не сильно отличается => сезонный эффект слабый

стационарность ряда.

IMEX(t) = (mean + trend * t) + error

=> Differencing method: Для исключения нестацинарности ряда

Применеия ARIMA для прогнозирования индекса Московской биржи:

ARIMA имеет следуюшие параметры:

• AR: p

• I: d (количество “difference” для приведения исходного нестационарного ряда в стационарный)

• MA: q

=> IMEX(t) – IMEX(t-1) = ARMA (p , q)

Деление на train and test set:

train_imex <- window(IMEX, start = c(2000, 3), end = c(2014, 12))
test_imex <- window(IMEX, start = 2015)
length(test_imex)

[1] 39

Применение модели по autoarima:

fit1 <- auto.arima(train_imex)
pred1 <- forecast(fit1, h = length(test_imex))
plot(pred1)
lines(test_imex, col = "red")

legend('bottomright',c("Actual","Forecasted"),lty=c(1,1),lwd=c(1.5,1.5),col=c('red','blue'))

Через экспоненту:

fit2 <- auto.arima(log(train_imex))
pred2 <- forecast(fit2, h = length(test_imex))
plot(pred2)
lines(log(test_imex), col = "red")

legend('bottomright',c("Actual","Forecasted"),lty=c(1,1),lwd=c(1.5,1.5),col=c('red','blue'))

Добавление BoxCox преобразования:

L <- BoxCox.lambda(IMEX)
fit3 <- auto.arima(train_imex, lambda = L)
pred3 <- forecast(fit3, h = length(test_imex), lambda = L)
plot(pred3)
lines(test_imex, col = "red")

legend('bottomright',c("Actual","Forecasted"),lty=c(1,1),lwd=c(1.5,1.5),col=c('red','blue'))

Ручной подбор модели:

Проверка стационарности ряда:

adf.test(diff(log(IMEX)), alternative="stationary", k = 0)

p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  diff(log(IMEX))
Dickey-Fuller = -12.097, Lag order = 0, p-value = 0.01
alternative hypothesis: stationary

Ряд стационарный после 1 difference => d = 1. Когда мы указываем d = 1 в модели ARIMA то ARIMA сама получит стационарный ряд в разности логарифмов

ACF plots:

acf(log(IMEX))

Распад диаграммы ACF происходит очень медленно => ряд нестационарный (выше уже обсудили)

View ACF and PACF after regressing on the difference:

acf(diff(log(IMEX)))

pacf(diff(log(IMEX)))

fit4 <- arima(log(train_imex), c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12))
pred4 <- forecast(fit4, h = length(test_imex))
accuracy(log(test_imex), pred4$mean)

                ME      RMSE        MAE      MPE     MAPE      ACF1 Theil's U
Test set 0.0947091 0.1177479 0.09982258 1.230632 1.298106 0.7568122  5.179196

plot(pred4)
lines(log(test_imex), col = "red")

legend('bottomright',c("Actual","Forecasted"),lty=c(1,1),lwd=c(1.5,1.5),col=c('red','blue'))

Возврат прогоза через экспоненту к прогнозу исходных значений индекса МБ

imex_pred <- exp(pred4$mean)
accuracy(imex_pred, test_imex)

                ME     RMSE      MAE       MPE    MAPE      ACF1 Theil's U
Test set -195.5966 249.9426 205.8375 -10.20452 10.7077 0.8229259  2.626428

plot(imex_pred)
lines(test_imex, col = "red")

df <- data.frame(actual = test_imex, forecast = imex_pred)
df