1 Introduction

Halo, ini merupakan Learn By Building yang dilakukan di Algoritma Data Science. Kali ini saya akan melakukan analis terhadap data Time Series dan membuat model forecasting. Saya akan melakukan prediksi harga closing Bitcoin dan membuat model forecasting, tak lupa juga akan melakukan evaluasi pada model tersebut.

Bitcoin sendiri adalah bentuk mata uang digital atau kripto yang diperkenalkan pada tahun 2009 oleh seseorang atau kelompok yang menggunakan nama samaran Satoshi Nakamoto. Ini adalah mata uang kripto pertama yang diperkenalkan dan terkenal sampai hari ini.Berbeda dengan mata uang tradisional yang diatur oleh pemerintah atau bank sentral, Bitcoin beroperasi di atas teknologi blockchain. Blockchain adalah buku besar terdesentralisasi yang mencatat semua transaksi Bitcoin. Bitcoin menggunakan teknologi kriptografi untuk mengamankan transaksi dan mengontrol penciptaan unit baru.

Beberapa ciri utama Bitcoin termasuk desentralisasi (tidak dikendalikan oleh satu otoritas sentral), pseudonimitas (identitas pengguna terlindungi oleh alamat kripto), dan keterbatasan pasokan (hanya akan ada sekitar 21 juta Bitcoin yang dapat ditambang).Bitcoin dapat digunakan sebagai alat pembayaran, investasi, atau sebagai bentuk penyimpan nilai. Meskipun popularitasnya terus tumbuh, Bitcoin juga telah menjadi subjek perdebatan karena fluktuasi harganya yang signifikan dan kekhawatiran seputar regulasi dan keamanan.

Kita sudah mengetahui mengenai apa itu bitcoin, mari kita melakukan prediksi menggunakan data yang ada. Data harga Bitcoin saya dapatkan melalui Kaggle : https://www.kaggle.com/datasets/abhishek14398/bitcoin-prediction-dataset-bullrun/code

2 Data Preparation

2.1 Memanggil Library

# load library

library(dplyr) # data wrangling
## 
## 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(lubridate) # date manipulation
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(padr) # complete data frame
## Warning: package 'padr' was built under R version 4.3.2
library(zoo) # Missing value imputation
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(forecast) # time series library
## Warning: package 'forecast' was built under R version 4.3.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(TTR) # for Simple moving average function
library(MLmetrics) # calculate error
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:base':
## 
##     Recall
library(tseries) # adf.test
## Warning: package 'tseries' was built under R version 4.3.2
library(fpp) # data for forecasting: principles and practice
## Warning: package 'fpp' was built under R version 4.3.2
## Loading required package: fma
## Warning: package 'fma' was built under R version 4.3.2
## Loading required package: expsmooth
## Warning: package 'expsmooth' was built under R version 4.3.2
## Loading required package: lmtest
library(TSstudio) # mempercantik visualisasi timeseries
## Warning: package 'TSstudio' was built under R version 4.3.2
library(ggplot2)
library(tidyr)

2.2 Read data

btc <- read.csv(file = "BTC_USD_Prediction.csv")

#Mengecheck data
glimpse(btc)
## Rows: 2,787
## Columns: 7
## $ X                   <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, …
## $ Currency            <chr> "BTC", "BTC", "BTC", "BTC", "BTC", "BTC", "BTC", "…
## $ Date                <chr> "2014-03-14", "2014-03-15", "2014-03-16", "2014-03…
## $ Closing.Price..USD. <dbl> 124.6550, 126.4550, 109.5848, 119.6747, 122.3387, …
## $ X24h.Open..USD.     <dbl> 125.3047, 124.6550, 126.4550, 109.5848, 119.6747, …
## $ X24h.High..USD.     <dbl> 125.7517, 126.7585, 126.6657, 119.6750, 122.9363, …
## $ X24h.Low..USD.      <dbl> 123.56349, 124.63383, 84.32833, 108.05816, 119.005…

2.3 Data Wrangling

Agar dapat dilakukan analisis dan pemodelan dari data time series, terdapat beberapa syarat utama yang harus dipenuhi, antara lain:

  1. Data harus terurut berdasarkan periode waktunya (dari data terlampau sampai ke data terbaru).
  2. Interval waktunya harus tetap/sama.
  3. Tidak boleh ada waktu yang terlewat untuk setiap interval.
  4. Tidak boleh ada data yang missing

Berhubung kita ingn membuat data time series, saya akan membuang beberapa kolom yang tidak digunakan dalam pembuatan TS, kolom - kolom tersebut adalah X,Currency,X24h.Open..USD.,X24h.High..USD.,X24h.Low..USD.

Selain itu, untuk memudahkan membaca kolom, kita akan mengubah nama kolom Closing.Price..USD. menjadiClosingPriceUSD, dilanjut menggunakan arrange untuk mengurutkan interval waktu pada data

#Menggunakan `select` untuk membuang kolom
btc <- btc %>% select(-X,-`Currency`, -`X24h.Open..USD.`, -`X24h.High..USD.`, -`X24h.Low..USD.`) %>% 
 #menggunakan `rename` untuk mengubah nama kolom
   rename( ClosingPriceUSD = `Closing.Price..USD.`) %>% 
  #Mengurutkan data berdasarkan kolom date
  arrange(Date)

Selanjutnya, kita akan mengubah tipe kolom Date menjadi tipe data date dengan menggunakan as.Date

btc$Date <- as.Date(btc$Date)

#fungsi `tail` untuk mengecheck data
tail(btc, n = 100)

Lalu, kita mengecheck apakah ada data interval waktu yang terlewat

# Cek periode yang terlewat
all(btc$Date %in% seq(as.Date("2014-03-14"), as.Date("2021-10-29"), by = "days"))
## [1] TRUE

2.4 Chechking Missing Value

colSums(is.na(btc))
##            Date ClosingPriceUSD 
##               0               0

Okey, ternyata data kita sudah memenuhi syarat untuk dijadikan object TS alias time series

2.5 Object ts

Kita akan membuat object ts dari data btc dan disimpan dalam object btc_ts dengan data yang terkumpul adalah data harian dengan pola yang ingin dilihat adalah bulanan

  • data: harian
  • pola: bulanan
  • frequency: 30
btc_ts <- ts(btc$ClosingPriceUSD, frequency = 30, start = c(2014, 4), end = c(2021,10))

Kita memulai dari tahun 2015

btc_ts <- window(btc_ts,start = 2015)

3 EDA

Kita akan melihat grafik atau visualisasi object dari btc_ts dengan autoplot

btc_ts %>% autoplot()

📈 Insight : Dalam frame bulanan, Bitcoin mengalami kenaikan tren dari tahun 2014 sampai dengan 2016, setelah itu mengalami tren penurunan sampai 2021. Walaupun terdapat tren kenaikan kembali, selain itu terdapat seasonal pada data,terdapat multiseasonality pada data pula.

Untuk kita bisa melihat lebih lanjut mengenai komponen trend, seasonal, dan error kita bisa melakukan decomposition menggunakan decompose

btc_ts %>% decompose(type = "multiplicative") %>% autoplot()

4 Cross Validation

Untuk time series, cross validation memiliki ketentuan:

  • Data train: data yang lebih lampau
  • Data test : data yang lebih baru
  • Jumlah data train > jumlah data test (tidak harus sekian persen)

Oleh karena itu mari kita pisahkan data kita seperti berikut:

  • Data train: 2015 Januari - 2019 Oktober (58 bulan) simpan dengan nama object btc_train
  • Data test : 2019 November - 2021 oktober (24 bulan) simpan dengan nama object btc_test
# train = selain 24 data terakhir
btc_train <- head(btc_ts, n = -24)

# test = 24 data terakhir
btc_test <- tail(btc_ts, n = 24)

5 Build model

Data kita memiliki Trend dan juga Seasonal, itulah mengapa kita memutuskan untuk membuat model SARIMA

5.1 SARIMA

Seasonal ARIMA adalah metode arima dimana object time series yang ada memiliki pola seasonal. Tahapan dalam melakukan pemodelan menggunakan SARIMA mirip dengan tahap pemodelan ARIMA dengan sedikit perbedaan. Sebelum membuat model SARIMA, kita perlu membuat data kita menjadi stationery terlebih dahulu karena data btc_ts kira belum stationery tapi kita bisa check terlebih dahulu

5.1.1 Stasionarity

Untuk mengetahui secara pasti apakah data kita sudah stasioner atau belum, kita dapat melakukan uji asumsi dengan adf.test() (Augmented Dickey-Fuller test)

Menguji stasioneritas data:

  • H0: data tidak stasioner
  • H1: data stasioner
  • alpha = 0.05

Kondisi: - p-value > alpha: data belum stasioner - p-value < alpha: data sudah stasioner.

Dari hasil perhitungan di atas, apakah data kita sudah dianggap stasioner? belum, karena p-value > alpha (0.05)

adf.test(btc_train)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  btc_train
## Dickey-Fuller = -2.5311, Lag order = 5, p-value = 0.3547
## alternative hypothesis: stationary

Dari hasil perhitungan di atas, Data kita belum stationary, karena p-value (0.3547) > alpha (0.05)

5.1.1.1 Seasonal Differencing

Untuk membuat data menjadi stasioner, solusinya adalah dengan differencing. Differencing adalah proses mengurangi data saat ini dengan data sebelumnya. Data sebelumnya disebut dengan lag

Untuk differencing seasonal, kita akan gunakan parameter lag yang nilainya mengikuti frekuensi data yang kita miliki

# melakukan satu kali differencing untuk menghilangkan efek seasonal (1x)
btc_train %>% 
  diff(lag = 30) %>% # differencing seasonal
  diff() %>%         # differencing trend (differencing biasa)
  autoplot()

Saat ini data kita menurut adf.test() sudah stationary, dan kita bisa melihat visualisasi dengan autoplot()

btc_train %>% 
  diff(lag = 30) %>% # differencing seasonal
  diff() %>%  #diff trend
  adf.test()
## Warning in adf.test(.): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  .
## Dickey-Fuller = -4.2626, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

5.1.1.2 Plot PACF & ACF

Tentukan p dan q serta P dan Q dari plot PACF & ACF menggunakan tsdisplay().

  • p dan q: 0 sampai 5 lag pertama yang keluar dari garis biru
  • P dan Q: lag yang keluar dari garis biru sesuai dengan kelipatan frekuensi
# untuk menampilkan banyak lag yang muncul di plot acf dan pacf bisa tambahkan `lag.max`
btc_train %>% 
  diff(lag = 30) %>% # diff untuk seasonal
  diff() %>% # diff trend
  tsdisplay(lag.max = 90)

Untuk menentukan p dan q kita melihat 5 lag pertama

  • p: 0,3,5
  • d: 1
  • q: 0,3,5

Untuk P dan Q (seasonal) kita melihat lag berdasarkan kelipatan frequency.

  • P : 0,30
  • D : 1
  • Q : 0,30

5.1.2 Build Model

Dalam membuat model ARIMA kita dapat menggunakan dua fungsi, yaitu Arima() dan auto.arima().

Fungsi Arima() digunakan untuk membuat model secara manual, parameter:

  • y = object time series
  • order = nilai (p,d,q)
  • seasonal = nilai (P,D,Q)
# order = (p,d,q), seasonal = (P,D,Q)
btc_sarima1 <- Arima(y = btc_train, order = c(3,1,0), seasonal = c(0,1,0))
btc_sarima2 <- Arima(y = btc_train, order = c(5,1,5), seasonal = c(0,1,0))
# auto.arima dengan d = 1
btc_auto <- auto.arima(y = btc_ts, d = 1, D = 1)
btc_auto
## Series: btc_ts 
## ARIMA(0,1,0)(2,1,0)[30] 
## 
## Coefficients:
##          sar1     sar2
##       -0.8040  -0.3880
## s.e.   0.0794   0.0917
## 
## sigma^2 = 3768:  log likelihood = -890.25
## AIC=1786.5   AICc=1786.65   BIC=1795.7

6 Forecast

# forecast, h = 24
btc_forecast_1 <- forecast(btc_sarima1, h = 24)
btc_forecast_2 <- forecast(btc_sarima2, h = 24)
btc_forecast_auto <- forecast(btc_auto, h = 24)

Visualisasi

btc_train %>% 
  autoplot() +
  autolayer(btc_forecast_1$mean, series = "forecast manual 1") +
  autolayer(btc_forecast_2$mean, series = "forecast manual 2") +
  autolayer(btc_forecast_auto$mean, series = "forecast auto") +
  autolayer(btc_test, series = "data test")

btc_test %>% 
  autoplot() +
  autolayer(btc_forecast_1$mean, series = "forecast manual 1") +
  autolayer(btc_forecast_2$mean, series = "forecast manual 2") +
  autolayer(btc_forecast_auto$mean, series = "forecast auto")

Insight:

  • Dari hasil visualisasi di dapatkan beberapa pola gagal ditangkap dengan baik oleh model
  • ARIMA baik ketika forecast period cenderung pendek

7 Evaluation

Menggunakan accuracy() dari package forecast dan melihat nilai MAPEnya

accuracy(btc_forecast_1)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -7.700573 75.63095 47.02639 -1.498383 6.663547 0.2033946
##                    ACF1
## Training set 0.02030032
accuracy(btc_forecast_2)
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -4.758779 67.25362 43.82558 -0.9698258 6.199298 0.1895507
##                    ACF1
## Training set 0.01552822
accuracy(btc_forecast_auto)
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set -7.28199 55.79604 35.77795 -1.248912 5.331073 0.1704362 0.01172023

Dari hasil evaluasi, model dengan MAPE terkecil adalah btc_forecast_auto, yaitu ARIMA(0,1,0)(2,1,0)

8 Assumption

Asumsi pada time series diujikan untuk mengukur apakah residual model yang terbentuk sudah cukup baik untuk menggambarkan informasi pada data.

8.1 No-autocorrelation residual

Untuk mengecek ada/tidaknya autokorelasi pada hasil forecasting time series bisa menggunakan uji Ljung-box dengan menggunakan fungsi Box.test(residual model, type = "Ljung-Box")

  • \(H_0\): residual tidak berkorelasi
  • \(H_1\): residual berkorelasi

Yang diinginkan: p-value > 0.05 (alpha), no-autocorrelation

Box.test(btc_forecast_auto$residuals, type = "Ljung-Box")
## 
##  Box-Ljung test
## 
## data:  btc_forecast_auto$residuals
## X-squared = 0.026513, df = 1, p-value = 0.8707

Kesimpulan: model kita sudah memenuhi asumsi no-autocorrelation