Wczytanie danych z pliku. W danych zawarty jest opis, dlatego pierwsze 15 wierszy zostało pominiete.

library(readxl)
library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(zoo)
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
df = read_excel('Dow_Jones.xlsx', skip=15, col_names = c("month", "index_value"))
df

Podstawowe statystyki:

summary(df)
##      month                         index_value    
##  Min.   :1968-08-01 00:00:00.00   Min.   :     0  
##  1st Qu.:1974-08-01 00:00:00.00   1st Qu.:  3672  
##  Median :1980-08-01 00:00:00.00   Median :  3754  
##  Mean   :1980-08-26 05:33:50.44   Mean   :  4133  
##  3rd Qu.:1986-10-01 00:00:00.00   3rd Qu.:  3850  
##  Max.   :1992-10-01 00:00:00.00   Max.   :120000  
##  NA's   :2                        NA's   :2
str(df)
## tibble [291 × 2] (S3: tbl_df/tbl/data.frame)
##  $ month      : POSIXct[1:291], format: "1968-08-01" "1968-09-01" ...
##  $ index_value: num [1:291] 3645 3626 3634 3620 3607 ...
nrow(df)
## [1] 291

Po przeglądnieciu statysyk można zauważyc, że mamy tutaj do czynienia z szeregiem czasowym określającym wartość indeksu Dow-Jones na koniec każdego miesiąca. W data framie występują dwie wartości brakujące. Zobaczmy jak wygląda przebieg tego szeregu.

ggplot(data=df, aes(x=month, y=index_value)) + geom_line()
## Warning: Removed 2 row(s) containing missing values (geom_path).

ggplot(data=df,aes(y=index_value)) + geom_boxplot()
## Warning: Removed 2 rows containing non-finite values (stat_boxplot).

ggplot(data=df, aes(x=index_value)) + geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing non-finite values (stat_bin).

Na każdym z trzech powyżej przedstawionych wykresów można zauważyć, występujące wartości odstające w pobliżu 120000 oraz 0, które zabużają wartość indeksu. Wartości te są abstrakcyjne, ponieważ w historii indeks nigdy nie osiągnął 0 ani 120000. W następnych krokach szereg zostanie oczyszczony z wartości odstających metodą wskaźnika IQR

Q1 = quantile(df$index_value, .25, na.rm=TRUE)
Q3 = quantile(df$index_value, .75, na.rm=TRUE)
IQR = IQR(df$index_value, na.rm=TRUE)

df = df %>%
  mutate(index_value = ifelse(df$index_value > (Q1 - 1.5*IQR) &
                              df$index_value < (Q3 + 1.5*IQR), index_value, NA) 
         )
sum(is.na(df$index_value))
## [1] 5

Rozkład indeksu nie należy do rozkładu normalnego

ggplot(data=df, aes(x=index_value)) + geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 5 rows containing non-finite values (stat_bin).

ggplot(data=df,aes(y=index_value)) + geom_boxplot()
## Warning: Removed 5 rows containing non-finite values (stat_boxplot).

Wartości odstające zostały usunięte, aczkolwiek w zbiorze nadal występują braki danych co widać też na wykresie poniżej.

ggplot(data=df, aes(x=month, y=index_value)) + geom_line()
## Warning: Removed 2 row(s) containing missing values (geom_path).

Ponieważ jest to szereg czasowy i są to dane historyczny pozbycie sie 5 miesięcy ze zbioru może wiąząć się z niepoprawnymi wynikami w predykcji, tego nie chcemy, dlatego wartości odstające zostaną usunięte.Do zastąpienie brakujących wartości uzyjemy funkcji aproksymacyjnej z pakietu zoo.

i_ts = ts(df$index_value,start=c(1968,8),end=c(1992,10),frequency=12)
i_ts =  na.approx(i_ts);
sum(is.na(i_ts))
## [1] 0
autoplot(i_ts)

Z wykresu liniowego możemy wysunąć stwierdzenie, że występuje trend rosnący z pewną powtażającą się sezonowością addytywną. Zobaczmy wykresy ACF oraz PACF potwierdzenia hipotezy.

acf(df$index_value, na.action = na.pass)

Na wykresie ACF można zauwazyć powolny spadek wraz ze wzortem lagów, oznacza to że zbiór posiada trend. Co do sezonowatosci ciezko tutaj stwierdzić, ponieważ brak jest lagów ktore odbiegały by rozmiarem od innych, aczkolwiek można zauważyc, że peaki tworzą “delikatne fale”

pacf(df$index_value, na.action = na.pass)

Na wykrecie PACF mozemy zauwzyc znaczaca korelacje na 1 lagu, która znacząco spadana na juz na drugi logu.

Dokonajmy dekompozycji szeregu

decomposed =  decompose(i_ts)

plot(decomposed)

Usunmy trend oraz sezonowość, zeby zobrazować sobie jak te dwie cechy wplywaja na nasz szereg. Mozna zauwzyc, ze po odjeciu sezonowatosci od szeregu szereg się nie zmienil, może to oznaczyc jego brak wystepowania.

tsWithoutTrend = i_ts/decomposed$trend
plot(i_ts-decomposed$trend)

plot(i_ts-decomposed$seasonal)

Sprawdzmy jak wygladaja dane na danych intervałach 1,2,5,10 lat, w celu potwierdzenia braku sezonowatosci

startYear = 1968

while (startYear < 1992) {
  plot(window(i_ts, start=c(startYear,1), end=c(startYear+1,1)))
  startYear =startYear+1
}
## Warning in window.default(x, ...): 'start' value not changed

startYear = 1968

while (startYear < 1992) {
  plot(window(i_ts, start=c(startYear,1), end=c(startYear+2,1)))
  startYear =startYear+2
}
## Warning in window.default(x, ...): 'start' value not changed

startYear = 1968

while (startYear < 1992) {
  plot(window(i_ts, start=c(startYear,1), end=c(startYear+5,1)))
  startYear =startYear+5
}
## Warning in window.default(x, ...): 'start' value not changed

## Warning in window.default(x, ...): 'end' value not changed

startYear = 1968

while (startYear < 1998) {
  plot(window(i_ts, start=c(startYear,1), end=c(startYear+10,1)))
  startYear =startYear+10
}
## Warning in window.default(x, ...): 'start' value not changed

## Warning in window.default(x, ...): 'end' value not changed

W zadnym okresie nie wystepuje podobienstwo pomiedzy wykresami, zatem definitywnie utwierdzam się w przekonaniu, że w szeregu nie wystepuje sezonowość. Wiec szereg ma trend rosnacy bez sezonowosci.

Przejdzmy do modelowania szeregu

Podział danych na testowe i treningowe z założeniem 4 elementowego zbioru testowego

train = head(i_ts,length(i_ts)-4)
test = tail(i_ts, 4)

is_no_trend = diff(i_ts)
train_no_trend = head(is_no_trend,length(is_no_trend)-4)
test_no_trend = tail(is_no_trend, 4)

Model nr 1. Modelowanie szeregu za pomoca metody Simple Exponential Smoothing. Metoda SES wykorzystywana jest do szeregow nie wykazujących trendu oraz sezonowosci. W przypadku naszych danych trend jest widoczny, dlatego musimy się go pozbyć

autoplot(train_no_trend)

Jako pierwszy model wykonajamy sobie model z podstawowymi parametrami(pozwalamy modelowi wyestymowac najlepszy za nas) w oknie o wielkosci 4

basic_model_ses = ses(train_no_trend,  h=4)
basic_model_ses$model
## Simple exponential smoothing 
## 
## Call:
##  ses(y = train_no_trend, h = 4) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
## 
##   Initial states:
##     l = 0.9572 
## 
##   sigma:  22.4518
## 
##      AIC     AICc      BIC 
## 3401.310 3401.396 3412.278
autoplot(basic_model_ses)

autoplot(train_no_trend) +
  autolayer(basic_model_ses$mean,
            color = "red") 

 autoplot(test_no_trend) +
  autolayer(basic_model_ses, alpha=0.5) +
  ggtitle("Predicted vs actuals")

Rozkład residuów modelu sprowadza się do rozkładu normalnego

hist(basic_model_ses$residuals)

plot(basic_model_ses$fitted, basic_model_ses$residuals)
title("fitted vs residuals")

Zoptymalizujemy parametr alpha w modelu metodą minimalizacji wskaznika RMSE w celu uzuskania lepszego przedziału rozwiązań

# comparing our model
alpha = seq(.01, .99, by = .01)
RMSE = c()
for(i in seq_along(alpha)) {
  fit = ses(train_no_trend, alpha = alpha[i],
             h = 4)
  RMSE[i] = accuracy(fit,
                      test_no_trend)[2,2]
}
alpha_fit = data_frame(alpha, RMSE)
## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## ℹ Please use `tibble()` instead.
alpha_min = filter(alpha_fit,
                    RMSE == min(RMSE))
 
alpha_min

Najniższe RMSE, które i tak jest relatywnie duże zostało uzyskane zostało za pomocą alpha = 0.99

ggplot(alpha_fit, aes(alpha, RMSE)) +
  geom_line() +
  geom_point(data = alpha_min,
             aes(alpha, RMSE),
             size = 2, color = "red")

Model z parametrem alpha = 0.99

final_model_ses = ses(train_no_trend, alpha=.99, h=4)
autoplot(final_model_ses)

final_model_ses$model
## Simple exponential smoothing 
## 
## Call:
##  ses(y = train_no_trend, h = 4, alpha = 0.99) 
## 
##   Smoothing parameters:
##     alpha = 0.99 
## 
##   Initial states:
##     l = -18.755 
## 
##   sigma:  30.8717
## 
##      AIC     AICc      BIC 
## 3581.475 3581.518 3588.787

Porównanie modeli

Przedstawienie sredniej modelu na wykresie szeregu

autoplot(train_no_trend) +
  autolayer(final_model_ses$mean,
            color = "red") 

 autoplot(test_no_trend) +
  autolayer(final_model_ses, alpha=0.5) +
  ggtitle("Predicted vs actuals")

Rozkład residuów modelu sprowadza się do rozkładu normalnego

hist(final_model_ses$residuals)

plot(final_model_ses$fitted, final_model_ses$residuals)
title("fitted vs residuals")

accuracy(basic_model_ses, test_no_trend)
##                         ME     RMSE      MAE      MPE     MAPE      MASE
## Training set  1.979394e-04 22.37315 16.72553     -Inf      Inf 0.7043591
## Test set     -1.670725e+01 26.31197 25.22862 101.4891 101.4891 1.0624484
##                     ACF1 Theil's U
## Training set  0.04340117        NA
## Test set     -0.11303131 0.9850212
accuracy(final_model_ses, test_no_trend)
##                       ME     RMSE      MAE     MPE    MAPE      MASE       ACF1
## Training set  0.04575685 30.76358 23.18247     NaN     Inf 0.9762792 -0.5027572
## Test set     -9.95059484 22.63188 21.85030 90.9787 90.9787 0.9201775 -0.1130313
##              Theil's U
## Training set        NA
## Test set     0.7653049

Porównianie otrzymanych modeli na wykresie

p1 <- autoplot(basic_model_ses) +
  coord_cartesian(ylim = c(-120, 120)) +
  ggtitle("Original SES Model") 
p2 <- autoplot(final_model_ses) +
  coord_cartesian(ylim = c(-120, 120)) +
  ggtitle("Optimal SES Model") 


gridExtra::grid.arrange(p1, p2, nrow = 1)

Dzięki optymalizacji udało się zopytamalizować RMSE z 25 do 22 oraz wskaznik predykcji mape ze 100 do 90, co nadal nie jest zbyt zadawalącym wynikiem. W Obu przypadkach resiuda sprowadzaja sie do rozkladu normalnego, natomiast wariancja w modelu numer 1 jest blizsza 0 ze wzledu na mniejszy przedział ufności możliwych wyników. Ze wzgledu na to, że dany model oparty jest o szereg z usunietym trendem, wybrał bym model numer 2, przede wszystkim ze wzgledu na lepszy MAPE oraz RMSE, ale i wiekszy przedzial ufnosci, na wykresie można zuważyć można wiele odchylen, przez co moim zdaniem taki scenariusz jest bardziej pradopodobny.

Model nr 2. Model Holta.

Drugi model wykonany z zastosowaniem metody Holta. Wybór padł na tą metodę ze względu na występowanie trendu w szeregu, przez co idealnie pasuję do przedstawonych danych. Tak samo jak w przypadku pierwszym wykonajmy podstawowy model w oknie 4 elementowym w celu porownania wyników

holt_model_basic = holt(train, h=4)
holt_model_basic$model
## Holt's method 
## 
## Call:
##  holt(y = train, h = 4) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 6e-04 
## 
##   Initial states:
##     l = 3625.6483 
##     b = 0.8272 
## 
##   sigma:  22.5246
## 
##      AIC     AICc      BIC 
## 3418.029 3418.243 3436.327
autoplot(holt_model_basic)

autoplot(test) +
  autolayer(holt_model_basic, alpha= 0.5) +
  ggtitle("Predicted vs actuals")

hist(holt_model_basic$residuals)

plot(holt_model_basic$fitted, holt_model_basic$residuals)
title("fitted vs residuals")

autoplot(train) +
  autolayer(holt_model_basic$mean,
            color = "red")

W modelu Holta jako parametr, który pozwala pozbyc się bledów z danych treningowych jest parametr beta, ktory przymuje wartosci od 0.0001 do 5

beta = seq(.0001, .5, by = .001)
RMSE = c()
for(i in seq_along(beta)) {
  fit = holt(train, 
              beta = beta[i],
              h = 4)
  RMSE[i] = accuracy(fit,
                      test)[2,2]
}

beta_fit = data_frame(beta, RMSE)
beta_min = filter(beta_fit,
                   RMSE == min(RMSE))
 
beta_min
ggplot(beta_fit, aes(beta, RMSE)) +
  geom_line() +
  geom_point(data = beta_min,
             aes(beta, RMSE),
             size = 2, color = "red") +
  ggtitle("beta vs RMSE")

final_model_holt = holt(train,
                    beta = 1e-04,
                    h = 4)

final_model_holt$model
## Holt's method 
## 
## Call:
##  holt(y = train, h = 4, beta = 1e-04) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 3630.49 
##     b = 1.0142 
## 
##   sigma:  22.5059
## 
##      AIC     AICc      BIC 
## 3415.552 3415.693 3430.189
autoplot(test) +
  autolayer(final_model_holt, alpha= 0.5) +
  ggtitle("Predicted vs actuals")

hist(final_model_holt$residuals)

plot(final_model_holt$fitted, final_model_holt$residuals)
title("fitted vs residuals")

autoplot(train) +
  autolayer(final_model_holt$mean,
            color = "red")

Statystyki opisujące modele

accuracy(holt_model_basic, test)
##                       ME     RMSE      MAE          MPE      MAPE      MASE
## Training set   0.1874348 22.36712 16.73841  0.003120699 0.4454027 0.2807260
## Test set     -21.9043234 37.61722 30.47325 -0.567907550 0.7856141 0.5110779
##                    ACF1 Theil's U
## Training set 0.04082631        NA
## Test set     0.20194560  1.525088
accuracy(final_model_holt, test)
##                        ME     RMSE      MAE          MPE      MAPE      MASE
## Training set  -0.01010764 22.34851 16.71322 -0.002187622 0.4447283 0.2803034
## Test set     -22.28559922 37.97747 30.77827 -0.577721842 0.7934909 0.5161934
##                    ACF1 Theil's U
## Training set 0.04156291        NA
## Test set     0.20226164  1.541196
p1 <- autoplot(holt_model_basic) +
  coord_cartesian(ylim = c(3500, 4030)) +
  ggtitle("Original Holt's model") 
 
p2 <- autoplot(final_model_holt) +
  coord_cartesian(ylim = c(3500, 4030)) +
  ggtitle("Optimal Holt's model")

gridExtra::grid.arrange(p1, p2, nrow=1)

Porównąjąc model podstawowy do zopytamlizowanego, można zauwazyc, ze przedział ufnosci dla mozliwych wyników w przypadku modelu zoptymalizowanego jest mniejszy. RMSE w przypadku pierwszym jak i drugim jest na poziomnym poziomie = 37, natomiast wskaznik MAPE predykcji jest na poziomie ~ 8%, co swiadczy o dobrej predykcji dalszych wyników. W obu przypadkach rozkład residuów sprowadza się do rozkładu normalnego, oba modele znacząco się od siebie nie rożnią.

PODSUMOWANIE WYNIKÓW

Podsumowując, w modelu z wykorzystaniem SES otrzymaliśmy wniki z lepszym RMSE = 22, w przypadku modelu Holta wynosi ono ~38, natomiast wznacząco różnią się wskazniki predkycji MAPE ponieważ w pierwszym przypadku najlepszy wynik uzysakliśmy z MAPE=90, w modelu holta jest o wiele lepiej = 0.8. Podejżewam, ale nie jestem w stanie potwierdzić tej hipotezy, ze tak slaby MAPE w modelu nr 1 jest wynikiem usunięcia trendu z danych, co może wiazać się z wiekszą trudnością dalszych predykcji(w przypadku numer 2 znamy trend, co juz w pewnym stopniu ułatwia dalsze predykcje). Ze względu na brak potrzeby usuwania trendu z szeregu uważam, że nie mniej jednak algorytm Holta dał lepsze wyniki.