Moving Average & Exponential Smoothing Forecasts

Data

Data yang digunakan adalah harga tutup harian komoditas emas dalam USD. Data terbagi menjadi training set (1 Januari 2001 - 31 Desember 2019 [4748 hari]) dan testing set (1 Januari 2020 - 21 Januari 2022 [752 hari]).

Metodologi

Data akan dieksplorasi dan didekomposisi untuk dilihat pergerakan tren dan musimannya. Kemudian, data-data yang hilang nilainya akan diimputasi menggunakan metode interpolasi. Metode peramalan deret waktu menggunakan SES, DES, SMA, dan DMA pada dataset untuk dicari parameter-parameter terbaiknya. Pemilihan model terbaik akan didasarkan pada metrik MSE dan MAE-nya.

Pembersihan Data

Penyesuaian format data dilakukan dengan cara mengubah tipe data kolom X pada df menjadi datetime. Kemudian, dibuatlah objek time-series dari data tanggal dan harga tutup pada df. Lalu, data hilang diimputasi menggunakan interpolasi.

df <- read.csv("GOLD.csv")
df$X <- as.Date(df$X)
str(df)

## 'data.frame':    5449 obs. of  7 variables:
##  $ X            : Date, format: "2000-08-30" "2000-08-31" ...
##  $ GC.F.Open    : num  274 275 277 NA 276 ...
##  $ GC.F.High    : num  274 278 277 NA 276 ...
##  $ GC.F.Low     : num  274 275 277 NA 276 ...
##  $ GC.F.Close   : num  274 278 277 NA 276 ...
##  $ GC.F.Volume  : int  0 0 0 NA 2 0 125 0 0 0 ...
##  $ GC.F.Adjusted: num  274 278 277 NA 276 ...

library(xts)

## Warning: package 'xts' was built under R version 4.1.2

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(tsbox)

## Warning: package 'tsbox' was built under R version 4.1.2

ts.data <- ts_ts(xts(df$GC.F.Close, df$X))
class(ts.data)

## [1] "ts"

sum(is.na(ts.data))

## [1] 2447

library(imputeTS)

## Warning: package 'imputeTS' was built under R version 4.1.2

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## 
## Attaching package: 'imputeTS'

## The following object is masked from 'package:zoo':
## 
##     na.locf

ts.data.imp <- na_interpolation(ts.data)
plot(ts.data.imp, main = "Gold Closing Price", col = "black")

Tidak seluruh data akan relevan untuk memprediksi waktu di masa yang akan datang, sehingga tidak seluruh data pula akan digunakan. Data yang digunakan untuk pemodelan akan dimulai dari minimal 15 tahun terakhir (1 Januari 2007). Terakhir, data dibagi menjadi data latih dan data uji. Pembagian data adalah latih 86% (Jan 2007 - Des 2019) dan uji 14% (Jan 2020 - 21 Jan 2022).

ts.data.sel <- window(ts.data.imp, start = 2007)
plot(ts.data.sel, main = "Selected Gold Closing Price", ylab = "Price (US$)")

ts.data.train <- window(ts.data.sel, end = 2020)
ts.data.test <- window(ts.data.sel, start = 2020)

start(ts.data.train); end(ts.data.train)

## [1] 2007.001

## [1] 2019.998

start(ts.data.test); end(ts.data.test)

## [1] 2020

## [1] 2022.057

length(ts.data.train)  # 4748

## [1] 4748

h <- length(ts.data.test); h  # 752

## [1] 752

Eksplorasi Data

Data uji cenderung merupakan kelanjutan dari pola data latih, sehingga pemilihan rentang waktu sudah cukup baik.

plot.train <- ts.data.sel
plot.train[time(ts.data.sel) >= start(ts.data.test)] <- NA
plot.test <- ts.data.test

plot.zoo(cbind(plot.train, plot.test), main = "Train-Test Split",
         plot.type = "single", col = c("black", "blue"), ylab = "Price (US$)")
legend("bottomright", legend = c("train", "test"), col = c("black", "blue"),
       lty = 1)

Hasil dekomposisi menunjukkan bahwa data mengandung unsur musiman yang berulang setiap tahunnya. Hal tersebut menunjukkan bahwa data memiliki siklus tahunan.

plot(decompose(ts.data.train))

Prasyarat Analisis

library(forecast)

## Warning: package 'forecast' was built under R version 4.1.2

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(Metrics)

## Warning: package 'Metrics' was built under R version 4.1.2

## 
## Attaching package: 'Metrics'

## The following object is masked from 'package:forecast':
## 
##     accuracy

library(TTR)

## Warning: package 'TTR' was built under R version 4.1.2

Fungsi myTS digunakan untuk mengubah data bertipe vektor menjadi objek ts dengan satuan penanggalan yang menyesuaikan dengan data uji. Sementara itu, fungsi myMetrics digunakan untuk secara otomatis mengeluarkan metrik MSE dan MAE untuk hasil-hasil prediksi.

myTS <- function (x) {
    return(ts(x, start = start(ts.data.test), end = end(ts.data.test),
              frequency = frequency(ts.data.test)))
}

myMetrics <- function (yTrue, yPred) {
    colName <- deparse(substitute(yPred))
    
    yTrue <- as.vector(yTrue)
    yPred <- as.vector(yPred)
    
    result <- as.data.frame(c(mse(yTrue, yPred), mae(yTrue, yPred)))
    colnames(result) <- colName
    rownames(result) <- c("MSE", "MAE")
    
    return(result)
}

Simple Moving Average

Parameter dari SMA adalah n. Parameter yang akan divalidasi adalah 4, 8, 12, 16, dan 20. Hasil uji menunjukkan bahwa metrik terbaik dimiliki oleh SMA(12).

sma4 <- forecast(SMA(ts.data.train, n = 4), h = h)$mean

## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series

sma8 <- forecast(SMA(ts.data.train, n = 8), h = h)$mean

## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series

sma12 <- forecast(SMA(ts.data.train, n = 12), h = h)$mean

## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series

sma16 <- forecast(SMA(ts.data.train, n = 16), h = h)$mean

## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series

sma20 <- forecast(SMA(ts.data.train, n = 20), h = h)$mean

## Warning in ets(object, lambda = lambda, biasadj = biasadj,
## allow.multiplicative.trend = allow.multiplicative.trend, : Missing values
## encountered. Using longest contiguous portion of time series

myMetrics(ts.data.test, sma4)

##          sma4
## MSE 80819.461
## MAE   265.055

myMetrics(ts.data.test, sma8)

##          sma8
## MSE 68083.104
## MAE   240.697

myMetrics(ts.data.test, sma12)  # best

##          sma12
## MSE 65130.2898
## MAE   235.2418

myMetrics(ts.data.test, sma16)

##          sma16
## MSE 72150.4135
## MAE   249.3069

myMetrics(ts.data.test, sma20)

##           sma20
## MSE 682481.9213
## MAE    650.3323

Double Moving Average

myDMA <- function (data, m) {
    result <- list()
    smooth.sma <- SMA(data, n = m)  # SMA
    smooth.dma <- SMA(smooth.sma, n = m)  # DMA = SMA of SMA
    
    # DMA forecast formula
    At <- 2*smooth.sma - smooth.dma
    Bt <- 2/(m - 1)*(smooth.sma - smooth.dma)
    
    data.dma <- At + Bt  
    forecast.dma <- c(NA, smooth.dma)
    
    bind <- cbind(data = c(data, NA),
                  smooth.dma = c(smooth.dma, NA),
                  forecast = forecast.dma)
    
    f = At[length(At)] + Bt[length(Bt)]*(1:h)
    
    future <- data.frame(data = rep(NA, h),
                         smooth.dma = rep(NA, h),
                         forecast = f)
    
    return(ts(future$forecast, start = start(ts.data.test),
              end = end(ts.data.test), frequency = frequency(ts.data.test)))
}

Fungsi myDMA dibuat untuk menggabungkan dua fungsi SMA yang berbeda dan menghasilkan pemodelan DMA yang sesuai. Parameter dari DMA adalah m. Parameter yang akan divalidasi adalah 7, 30, 90, 120, dan 365. Hasil uji menunjukkan bahwa metrik terbaik dimiliki oleh DMA(365).

dma7 <- myDMA(ts.data.train, m = 7)
dma30 <- myDMA(ts.data.train, m = 30)
dma90 <- myDMA(ts.data.train, m = 90)
dma120 <- myDMA(ts.data.train, m = 120)
dma365 <- myDMA(ts.data.train, m = 365)

myMetrics(ts.data.test, dma7)

##            dma7
## MSE 2725720.065
## MAE    1377.857

myMetrics(ts.data.test, dma30)

##          dma30
## MSE 39254.3904
## MAE   166.1042

myMetrics(ts.data.test, dma90)

##           dma90
## MSE 170149.3060
## MAE    391.9683

myMetrics(ts.data.test, dma120)

##         dma120
## MSE 36209.5287
## MAE   162.2277

myMetrics(ts.data.test, dma365)  # best

##         dma365
## MSE 23989.6959
## MAE   121.0686

plot(plot.train, ylim = c(min(plot.train, na.rm = T), max(plot.test)),
     ylab = "Price (US$)", main = "Best SMA vs Best DMA")
lines(plot.test, col = "blue")
lines(myTS(sma12), col = "purple", lwd = 3, lty = "dotted")
lines(dma365, col = "red", lwd = 3, lty = "dotted")
legend("bottomright", legend = c("SMA (12)", "DMA (365)"),
       col = c("purple", "red"), lty = "dotted", lwd = 3)

Simple Exponential Smoothing

Parameter dari SES adalah alpha. Parameter yang akan divalidasi adalah 0.2, 0.4, 0.6, 0.8, dan optimum. Alpha optimum pada kasus ini adalah 0.99. Hasil uji menunjukkan bahwa metrik terbaik dimiliki oleh SES(alpha = 0.99).

ses20 <- ses(ts.data.train, h = h, alpha = .20)$mean
ses40 <- ses(ts.data.train, h = h, alpha = .40)$mean
ses60 <- ses(ts.data.train, h = h, alpha = .60)$mean
ses80 <- ses(ts.data.train, h = h, alpha = .80)$mean

sesOpt <- ses(ts.data.train, h = h, alpha = NULL)
paste("alpha optimum: ", sesOpt$model$par[['alpha']])

## [1] "alpha optimum:  0.999899988982947"

sesOpt <- sesOpt$mean

myMetrics(ts.data.test, ses20)  # alpha .2

##          ses20
## MSE 88286.5546
## MAE   278.5275

myMetrics(ts.data.test, ses40)  # alpha .4

##          ses40
## MSE 84026.5034
## MAE   270.8966

myMetrics(ts.data.test, ses60)  # alpha .6

##          ses60
## MSE 82906.6152
## MAE   268.8613

myMetrics(ts.data.test, ses80)  # alpha .8

##          ses80
## MSE 82300.2877
## MAE   267.7533

myMetrics(ts.data.test, sesOpt)  # alpha .99  # best

##         sesOpt
## MSE 81759.9137
## MAE   266.7619

Double Exponential Smoothing

Parameter dari DES adalah alpha dan beta. Parameter yang akan divalidasi adalah kombinasi dari: alpha (0.25, 0.5, dan 0.75) dan beta (0.25, 0.5, dan 0.75). Alpha dan beta optimum pada kasus ini masing-masing adalah 1 dan 0. Hasil uji menunjukkan bahwa metrik terbaik dimiliki oleh DES(alpha = 0.5, beta = 0.75).

des2525 <- forecast(HoltWinters(ts.data.train, alpha = .25,
                                beta = .25, gamma = FALSE), h = h)$mean
des2550 <- forecast(HoltWinters(ts.data.train, alpha = .25,
                                beta = .50, gamma = FALSE), h = h)$mean
des2575 <- forecast(HoltWinters(ts.data.train, alpha = .25,
                                beta = .75, gamma = FALSE), h = h)$mean

des5025 <- forecast(HoltWinters(ts.data.train, alpha = .50,
                                beta = .25, gamma = FALSE), h = h)$mean
des5050 <- forecast(HoltWinters(ts.data.train, alpha = .50,
                                beta = .50, gamma = FALSE), h = h)$mean
des5075 <- forecast(HoltWinters(ts.data.train, alpha = .50,
                                beta = .75, gamma = FALSE), h = h)$mean

des7525 <- forecast(HoltWinters(ts.data.train, alpha = .75,
                                beta = .25, gamma = FALSE), h = h)$mean
des7550 <- forecast(HoltWinters(ts.data.train, alpha = .75,
                                beta = .50, gamma = FALSE), h = h)$mean
des7575 <- forecast(HoltWinters(ts.data.train, alpha = .75,
                                beta = .50, gamma = FALSE), h = h)$mean

desOpt <- HoltWinters(ts.data.train, alpha = NULL, beta = NULL, gamma = FALSE)
paste("alpha optimum:", desOpt$alpha)

## [1] "alpha optimum: 1"

paste("beta optimum:", desOpt$beta)

## [1] "beta optimum: 0"

desOpt <- forecast(desOpt, h = h)$mean

myMetrics(ts.data.test, des2525)  # alpha = .25; beta = .25

##         des2525
## MSE 2483712.141
## MAE    1311.799

myMetrics(ts.data.test, des2550)  # alpha = .25; beta = .50

##         des2550
## MSE 1588844.254
## MAE    1038.109

myMetrics(ts.data.test, des2575)  # alpha = .25; beta = .75

##        des2575
## MSE 60121.5424
## MAE   191.0497

myMetrics(ts.data.test, des5025)  # alpha = .50; beta = .25

##          des5025
## MSE 1113194.8199
## MAE     857.5382

myMetrics(ts.data.test, des5050)  # alpha = .50; beta = .50

##        des5050
## MSE 76980.8417
## MAE   213.1952

myMetrics(ts.data.test, des5075)  # alpha = .50; beta = .75  # best

##        des5075
## MSE 23858.5390
## MAE   117.7018

myMetrics(ts.data.test, des7525)  # alpha = .75; beta = .25

##         des7525
## MSE 910520.0957
## MAE    769.0104

myMetrics(ts.data.test, des7550)  # alpha = .75; beta = .50

##         des7550
## MSE 423558.1539
## MAE    504.0681

myMetrics(ts.data.test, des7575)  # alpha = .75; beta = .75

##         des7575
## MSE 423558.1539
## MAE    504.0681

myMetrics(ts.data.test, desOpt)  # alpha = 1; beta = 0

##         desOpt
## MSE 81759.6475
## MAE   266.7615

plot(plot.train, ylim = c(min(plot.train, na.rm = T), max(plot.test)),
     ylab = "Price (US$)", main = "Best SES vs Best DES")
lines(plot.test, col = "blue")
lines(sesOpt, col = "purple", lwd = 3, lty = "dotted")
lines(des5075, col = "red", lwd = 3, lty = "dotted")
legend("bottomright", legend = c("SES (.99)", "DES (.50; .75)"), col = c("purple", "red"),
       lty = "dotted", lwd = 3)

Seasonal Model

Walaupun pada eksplorasi data telah ditunjukkan bahwa data memiliki siklus musiman, fungsi HoltWinters tidak dapat digunakan.

try(addOpt <- HoltWinters(ts.data.train, seasonal = "additive"))

## Error in decompose(ts(x[1L:wind], start = start(x), frequency = f), seasonal) : 
##   time series has no or less than 2 periods

try(mulOpt <- HoltWinters(ts.data.train, seasonal = "multiplicative"))

## Error in decompose(ts(x[1L:wind], start = start(x), frequency = f), seasonal) : 
##   time series has no or less than 2 periods

Conclusion

Hasil uji metrik menunjukkan bahwa DES(.5, .75) adalah model terbaik. Akan tetapi, model DES(365) juga dapat dipilih sebagai alternatif karena metriknya tidak berbeda dan pendekatannya jauh lebih sederhana. Plot hasil prediksi ditunjukkan pada grafik di bawah.

myMetrics(ts.data.test, sma12)

##          sma12
## MSE 65130.2898
## MAE   235.2418

myMetrics(ts.data.test, dma365)

##         dma365
## MSE 23989.6959
## MAE   121.0686

myMetrics(ts.data.test, sesOpt)

##         sesOpt
## MSE 81759.9137
## MAE   266.7619

myMetrics(ts.data.test, des5075)  # best

##        des5075
## MSE 23858.5390
## MAE   117.7018

plot(ts.data.test, ylim = c(min(plot.test), max(plot.test)),
     main = "Gold Price", ylab = "Price (US$)")
lines(myTS(sma12), col = "purple", lty = "dashed", lwd = 2)
lines(dma365, col = "red", lty = "dashed", lwd = 2)
lines(sesOpt, col = "orange", lty = "dashed", lwd = 2)
lines(des5075, col = "blue", lty = "dashed", lwd = 2)
legend("topleft", legend = c("Test Data", "SMA", "DMA", "SES", "DES"),
       col = c("black", "purple", "red", "orange", "blue"),
       lty = c(1, 2, 2, 2, 2), cex = .75)