Halo!!

Pada LBB ini saya ingin menggunakan data pergerakan saham BMRI yang saya dapat dari https://www.kaggle.com/tiwill/saham-idx untuk kepentingan forecasting dengan menggunakan metode time series. Let’s get started.

Preparation

library(tidyverse)
library(lubridate) # date manipulation
library(padr) # TS padding
library(zoo) # imputasi missing value TS
library(fpp) # TS dataset
library(TSstudio) # TS interactive viz

library(forecast) # algoritma forecasting
library(TTR) # SMA function
library(tseries) # adf.test
rm(list=ls())
bmri <- read.csv("BMRI.csv")
head(bmri)

EDA & Data Preprocessing

Cek Tipe Data

glimpse(bmri)
## Rows: 460
## Columns: 25
## $ date                  <chr> "2019-07-29T00:00:00", "2019-07-30T00:00:00", "2~
## $ previous              <dbl> 7750, 7800, 7950, 7975, 7775, 7675, 7425, 7250, ~
## $ open_price            <dbl> 7825, 7800, 7875, 7875, 7700, 7600, 7300, 7275, ~
## $ first_trade           <dbl> 7825, 7800, 7900, 7875, 7700, 7600, 7300, 7275, ~
## $ high                  <dbl> 7825, 7975, 7975, 7900, 7750, 7650, 7400, 7400, ~
## $ low                   <dbl> 7675, 7725, 7800, 7725, 7600, 7325, 7100, 7275, ~
## $ close                 <dbl> 7800, 7950, 7975, 7775, 7675, 7425, 7250, 7350, ~
## $ change                <dbl> 50, 150, 25, -200, -100, -250, -175, 100, 150, -~
## $ volume                <dbl> 27603200, 43537800, 29504400, 48365500, 51863500~
## $ value                 <dbl> 214424905000, 343432807500, 233733622500, 376774~
## $ frequency             <dbl> 4210, 5010, 4017, 5692, 5284, 7669, 9479, 6849, ~
## $ index_individual      <dbl> 2350.3, 2395.5, 2403.1, 2342.8, 2312.7, 2237.3, ~
## $ offer                 <dbl> 7800, 7950, 7975, 7800, 7700, 7425, 7300, 7375, ~
## $ offer_volume          <dbl> 446400, 4295200, 6718100, 2181600, 282800, 26450~
## $ bid                   <dbl> 7775, 7925, 7925, 7775, 7675, 7400, 7250, 7350, ~
## $ bid_volume            <dbl> 300, 42200, 1000, 433500, 24000, 301400, 2680800~
## $ listed_shares         <dbl> 4.62e+10, 4.62e+10, 4.62e+10, 4.62e+10, 4.62e+10~
## $ tradeble_shares       <dbl> 4.62e+10, 4.62e+10, 4.62e+10, 4.62e+10, 4.62e+10~
## $ weight_for_index      <dbl> 4.62e+10, 4.62e+10, 4.62e+10, 4.62e+10, 4.62e+10~
## $ foreign_sell          <dbl> 18401800, 23262600, 22946100, 17684800, 20649100~
## $ foreign_buy           <dbl> 10264800, 18254500, 14466600, 10010700, 30265500~
## $ delisting_date        <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ~
## $ non_regular_volume    <dbl> 2890100, 759370, 554718, 7082271, 118606, 166182~
## $ non_regular_value     <dbl> 22439393700, 5859249700, 4447772400, 56296815128~
## $ non_regular_frequency <dbl> 25, 14, 13, 10, 5, 4, 4, 6, 3, 18, 7, 6, 9, 2, 1~

Kolom date perlu diubah menjadi tipe data datetime. Lalu untuk kepentingan forecasting time series, kita hanya akan ambil kolom date dan close. Karena pandemi COVID-19 menyerang Indonesia pertama kali pada Februari-Maret 2020 dan menyebabkan harga market drop, maka data yang digunakan pada forecasting kali ini dimulai dari April 2020.

Menyiapkan Data Frame Baru

bmri_new <- bmri %>% 
  select(date, close) %>% 
  mutate(date=ymd_hms(date)) %>% 
  arrange(date)
bmri_new <- bmri_new %>% 
  filter(date >= "2020-04-01")
head(bmri_new)

Cek missing value

bmri_new %>% 
  pad() %>% 
  anyNA()
## pad applied on the interval: day
## [1] TRUE

Ternyata terdapat missing values, mari kita lihat dimana letak missing values tersebut.

bmri_new %>% pad()
## pad applied on the interval: day

Terdapat missing values setiap 5 hari (kemungkinan hari sabtu dan minggu dimana pasar saham Indonesia tutup), mari kita isi nilainya di 2 hari tersebut dengan harga hari tepat di atasnya dengan asumsi hari sabtu dan minggu tidak terjadi transaksi (karena pasar tutup) dan harga menjadi stagnan.

bmri_new <- bmri_new %>%
  pad() %>% 
  fill(close, .direction = "down")
## pad applied on the interval: day

Kemudian kita cek lagi data bmri_new dengan pad

bmri_new %>% 
  pad() %>% 
  anyNA()
## pad applied on the interval: day
## [1] FALSE

Sampai sini data kita sudah memenuhi syarat time series yaitu :

[x] Data tidak boleh ada yang NA [x] Data harus terurut berdasarkan waktu [x] Data tidak boleh ada waktu yang missing

Mari masuk tahap berikutnya.

Membuat Object Time Series

bmri_ts <- ts(bmri_new$close, start = c(2019, 8), frequency = 30)

Membuat Data Train dan Data Test

bmri_test <- bmri_ts %>% tail(30)
bmri_train <- bmri_ts %>% head(-30)

Decompose

bmri_dc <- bmri_train %>% 
  decompose()

autoplot(bmri_dc)

Seasonal Analysis

autoplot(bmri_dc$x - bmri_dc$trend)

Trend Analysis

autoplot(bmri_dc$x - bmri_dc$seasonal)



Dari hasil analisis kedua plot di atas, tidak ditemukan adanya seasonal dan trend dalam pergerakan BMRI. Oleh karena itu selanjutnya kita dapat membuat model Holtwinters ataupun ARIMA untuk melakukan prediksi time series terhadap pergerakan harga BMRI.

Modelling

Holtwinters


Untuk model Holtwinters ini kita dapat set beta dan gamma menjadi FALSE

bmri_hw <- HoltWinters(bmri_train, beta = F, gamma = F)


## ARIMA

bmri_arima <-  auto.arima(bmri_train)

Forecasting

Model Holtwinters

hw_forecast <- forecast(object = bmri_hw, h = 30)
hw_forecast %>% 
  autoplot() +
  autolayer(bmri_test)

Model ARIMA

arima_forecast <- forecast(object = bmri_arima, h = 30)
arima_forecast %>% 
  autoplot() +
  autolayer(bmri_test)

Evaluation

Model Comparison

compare_forecast <- function(x, test){
  lapply(x, function(x) forecast::accuracy(x, test)["Test set", ]) %>%
    lapply(., function(x) x %>% t() %>% as.data.frame) %>% 
    bind_rows() %>% 
    mutate(model = names(x)) %>%
    select(model, everything())}

forecast_list <- list(
  "Holt_winters" = hw_forecast,
  "ARIMA" = arima_forecast)

compare_forecast(forecast_list, bmri_test)

Dari hasil komparasi kedua model di atas, terlihat dari nilai MAPE model Holtwinters sedikit lebih kecil dibanding model ARIMA. Oleh karena itu dapat kita katakan bahwa model Holtwinters sedikit lebih baik dalam memprediksi dalam case kali ini.

Visualisasi

test_forecast(actual = bmri_ts, 
              forecast.obj = hw_forecast, 
              train = bmri_train, 
              test = bmri_test)

Uji Asumsi


Asumsi pada time series diujikan untuk mengukur apakah residual yang peroleh dari hasil modeling sudah cukup baik untuk menggambarkan dan menangkap informasi pada data. Kita akan menggunakan tes pada residual data, mengapa menggunakan residual data? Karena kita dapat mendapatkan informasi dari data aktual maupun dari hasil prediksi menggunakan model.

Metode forecasting yang baik menghasilkan nilai residual berikut ini:
- Residual yang tidak berkorelasi. Apabila terdapat residual yang berkorelasi, artinya masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast.
- Residual berdistribusi normal di sekitaran nilai 0.

Untuk memastikan bahwa residual yang dihasilkan memiliki kriteria diatas, maka terdapat asumsi untuk mengujinya.

No Auto-correlation Residuals


Untuk melakukan tes asumsi NO-Auto Correlation, lebih pasti jika melakukan Uji Ljung-box dengan hipotesis:
- H0 : tidak terdapat autokorelasi pada residual
- H1 : terdapat autokorelasi pada residual

Harapan: p-value > alpha (0.05)

Box.test(bmri_arima$residuals, type = 'Ljung-Box')
## 
##  Box-Ljung test
## 
## data:  bmri_arima$residuals
## X-squared = 0.16514, df = 1, p-value = 0.6845

Normality Residuals


Untuk melakukan tes asumsi normality, lebih pasti jika melakukan uji shapiro test menggunakan fungsi shapiro.test(model$residuals) dengan hipotesis:
- H0 : residual menyebar normal
- H1 : residual tidak menyebar normal

Harapan: p-value > alpha (0.05)

shapiro.test(bmri_arima$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  bmri_arima$residuals
## W = 0.93153, p-value = 9.074e-13

Conclusion


1. Dari hasil analisis di atas, model Holtwinters sedikit lebih baik performanya dengan MAPE sebesar 2.814821 dibandingkan dengan model ARIMA.
2. Dari hasil uji asumsi, tidak ada auto korelasi pada residual (p-value >0.05), tetapi residual forecast masih tidak terdistribusi secara normal.