Forecasting Yahoo’s Stock Using Time Series

Introduction

Dalam studi kasus ini, saya menggunakan data yahoo stock yang akan digunakan untuk melakukan forecast harga saham yahoo pada 3 bulan ke depan, dengan mengidentifikasi data Close yang telah disediakan pada yahoo_stock.csv

Data Preparation

Library

library(readr)
library(dplyr) 
library(lubridate) 
library(ggplot2) 
library(TSstudio) 
library(padr) 
library(forecast) 
library(ggsci)
library(tidyr) 
library(zoo)

Import Data

stock <- read_csv("yahoo_stock.csv")
stock

cek tipe data:

glimpse(stock)
#> Rows: 1,825
#> Columns: 7
#> $ Date        <date> 2015-11-23, 2015-11-24, 2015-11-25, 2015-11-26, 2015-11-2…
#> $ High        <dbl> 2095.61, 2094.12, 2093.00, 2093.00, 2093.29, 2093.29, 2093…
#> $ Low         <dbl> 2081.39, 2070.29, 2086.30, 2086.30, 2084.13, 2084.13, 2084…
#> $ Open        <dbl> 2089.41, 2084.42, 2089.30, 2089.30, 2088.82, 2088.82, 2088…
#> $ Close       <dbl> 2086.59, 2089.14, 2088.87, 2088.87, 2090.11, 2090.11, 2090…
#> $ Volume      <dbl> 3587980000, 3884930000, 2852940000, 2852940000, 1466840000…
#> $ `Adj Close` <dbl> 2086.59, 2089.14, 2088.87, 2088.87, 2090.11, 2090.11, 2090…

cek missing value:

anyNA(stock)
#> [1] FALSE

Data Preprocess

Pilih data yang diperlukan:

closee <- stock %>% 
  select(c(Date,Close))
closee

Cek range data:

range(closee$Date)
#> [1] "2015-11-23" "2020-11-20"

Maka dapat dilakukan inspect periode terlewat:

all(seq.Date(from = as.Date("2015-11-23"), to = as.Date("2020-11-20"), by = "day") == closee$Date)
#> [1] TRUE

Mengurutkan data berdasarkan waktu menggunakan arrange()

closee <- closee %>% 
  arrange(closee$Date)

closee

Melakukan padding untuk memastikan interval data sama: pad() dari package padr dan replace nilai Missing Value menggunakan na.locf() agar di isi dengan nilai Close sebelum Missing:

library(padr)
library(zoo)

closee <- closee %>% 
  pad(start_val = min(closee$Date), end_val = max(closee$Date)) %>% 
  mutate(Close = na.locf(Close))

closee

Seasonality Analysis

Decomposition

Buat objek Time Series menggunakan fungsi ts() dengan besarnya frekuensi adalah 7412:

closee_ts <- ts(data = closee$Close, frequency = 7*4*12)

# melihat plot time series
closee_dc <- decompose(x= closee_ts)
closee_dc %>% autoplot()

Pola trend diatas masih belum cukup smooth, sehingga ada kemungkinan Multi-Seasonality, oleh karena itu kita perlu diubah ke dalam msts:

closee_msts <- msts(data = closee$Close,seasonal.periods = c(7*4, 7*4*12))

closee_msts_dc <- mstl(closee_msts)

closee_msts_dc %>% autoplot()

Insight:

Berdasarkan Plot di atas, kini Estimasi Trend dengan “Multiple Seasonality Time Series” lebih smooth dan jelas.

Model Fitting and Evaluation

Cross Validation

Splitting Data:

Data Test: Data yang dijadikan test untuk validasi dengan mengambil pola 3 bulan ( 7 x 4 x 3 )
Data Train: Total data kecuali Data Test (data yang akan dilatih)

test_close <- tail(closee_msts, 7*4*3)
train_close <- head(closee_msts, length(closee_msts)-length(test_close))

ARIMA Model

model_arima_msts <- stlm(train_close, method = "arima") 
arima_forecast_msts <- forecast(model_arima_msts, h=7*4*3)
plot(arima_forecast_msts)

Model Visualization

test_close %>% autoplot(series = "Actual") +
  autolayer(arima_forecast_msts$mean, series = "Predict Arima")

### Evaluation Performance Model

accuracy(object = arima_forecast_msts$mean, x = test_close)
#>                ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set 9.386997 94.43999 77.97485 0.1982769 2.272603 0.9005316  2.413182

Holt-Winters Model

Membandingkan dengan model lain dengan menggunakan Holt-Winters Model yang merupakan metode forecasting yang tepat digunakan untuk data yang memiliki efek trend dan seasonal.

model_hw_msts <- HoltWinters(train_close)

Forecasting

hw_forecast_msts <- forecast(model_hw_msts, h = 7*4*3)
plot(hw_forecast_msts)

Model Visualization

test_close %>% autoplot(series = "Actual") +
  autolayer(hw_forecast_msts$mean, series = "Predict Holt-Winters")

Evaluation Performance Model

accuracy(object = hw_forecast_msts$mean, x = test_close)
#>                 ME    RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -134.3427 162.428 136.6211 -3.991954 4.055708 0.8893351  4.296925

ETS Model

model_ets_msts <- stlm(train_close, method = "ets") 
ets_forecast_msts <- forecast(model_ets_msts, h=7*4*3)
plot(ets_forecast_msts)

### Model Visualization

test_close %>% autoplot(series = "Actual") +
  autolayer(ets_forecast_msts$mean, series = "Predict ETS")

### Evaluation Performance Model

accuracy(object = ets_forecast_msts$mean, x = test_close)
#>                ME     RMSE      MAE          MPE     MAPE     ACF1 Theil's U
#> Test set 2.404383 94.47807 78.22125 -0.005993502 2.284768 0.900688  2.423465

Visualization Actual vs Estimated

Persiapan data untuk melakukan Plotting :

accuracyData <- data.frame(Date= closee$Date %>% tail(7*4*3),
  actual = as.vector(test_close) ,
  HWForecastmsts = as.vector(hw_forecast_msts$mean) ,
  ArimaForecastmsts = as.vector(arima_forecast_msts$mean) ,
  ETSForecastmsts = as.vector(ets_forecast_msts$mean))

Visualization of Actual vs Estimated number of Close by All Models

accuracyData %>% 
 ggplot() +
  geom_line(aes(x = Date, y = actual, colour = "Actual"),size=0.5)+
  geom_line(aes(x = Date, y = HWForecastmsts, colour = "Holt Winter Model"),size=0.5)+
  geom_line(aes(x = Date, y = ArimaForecastmsts, colour = "Arima Model"),size=0.5)+
  geom_line(aes(x = Date, y = ETSForecastmsts, colour = "ETS Model"),size=0.5)+
  labs(title = "Closing Price by Actual vs All Models",x = "Date",y = "Visitor",colour = "")

# Prediction Performance

ARIMA, Holt-Winters & ETS Model

accuracy(object = arima_forecast_msts$mean, x = test_close)
#>                ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set 9.386997 94.43999 77.97485 0.1982769 2.272603 0.9005316  2.413182
accuracy(object = hw_forecast_msts$mean, x = test_close) 
#>                 ME    RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -134.3427 162.428 136.6211 -3.991954 4.055708 0.8893351  4.296925
accuracy(object = ets_forecast_msts$mean, x = test_close)
#>                ME     RMSE      MAE          MPE     MAPE     ACF1 Theil's U
#> Test set 2.404383 94.47807 78.22125 -0.005993502 2.284768 0.900688  2.423465

Dari model di atas yang memiliki nilai MAE terendah yaitu pada Model ARIMA sebesar 77.97, dengan nilai MAPE 2.27%.

Assumption Test

No-autocorrelation residual

# menggunakan Ljung-Box test
Box.test(model_arima_msts$residuals, type = "Ljung-Box")
#> 
#>  Box-Ljung test
#> 
#> data:  model_arima_msts$residuals
#> X-squared = 0.021171, df = 1, p-value = 0.8843

Pada Ljung-Box test diatas, p-value > alpha, dimana 0.88 > 0.05 yaitu gagal tolak H0, yang artinya residual/eror pada data tidak terdapat autocorrelation.

Normality of residual

shapiro.test(model_arima_msts$residuals)
#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  model_arima_msts$residuals
#> W = 0.83127, p-value < 2.2e-16

Pada Uji Normalitas ini nilai p-value < alpha, dimana 2.2e-16 < 0.05, maka residual belum berdistribusi normal.

hist(model_arima_msts$residual, breaks = 30)

Conclusion

  • Model time series dengan performa terbaik yang dilihat dari MAE (Mean Absolute Error) untuk dapat diinterpretasikan sebagai seberapa besar penyimpangan hasil prediksi terhadap nilai aktualnya adalah model ARIMA dengan MAE sebesar 77.97485 dan dilihat dari nilai MAPE yang dapat menunjukkan seberapa besar penyimpangannya dalam bentuk persentase sebesar 2.27%
  • Pada uji asumsi autokorelasi, dari hasil Box-Ljung test model ARIMA ini mendapatkan nilai p-value sebesar 0.8843 yang mengindikasikan bahwa residual tidak terdapat autokorelasi.
  • Pada uji normalitas, dari hasil saphiro test model ARIMA ini mendapatkan nilai p-value sebesar < 2.2e-16 yang mengindikasikan bahwa residual tidak terdistribusi secara normal.