Double Moving Average (DMA)

Metode DMA merupakan metode yang menghitung rata-rata bergerak yang kedua atau menghitung rata-rata bergerak dari rata-rata bergerak SMA. DMA diperlukan untuk dapat mengakumulasi dan mengakomodir trend yang ada (Jana et al. 2015). \[ \text{Proses Pemulusan:} \\S'_t=\frac{X_t+X_{t-1}+...+X_{t-N-1}}{N} \\S''_t=\frac{S'_t+S'_{t-1}+...+S'_{t-N-1}}{N}\\ \text{Forecasting:} \\F_{t-m}=a_t+b_tm \\\text{dimana:} \\a_t=2S'_t-S''_t \\b_t=\frac{2}{N-1}(S'_t-S''_t) \]

Double Exponential Smoothing (DES)

DES merupakan metode peramalan data yang mengalami trend kenaikan dan apabila data yang digunakan semakin banyak dalam perhitungan peramalannya maka error peramalannya akan semakin kecil, begitu juga sebaliknya. \[ \text{Proses Pemulusan:} \\\tilde{y}_T=\lambda y_{T-1} \\\tilde{y}_T^{(2)}=\lambda \tilde{y}_T^{(1)}+(1-\lambda)\tilde{y}_{T-1}^{(2)} \\\text{Dengan initial value dari:} \\\tilde{y}_0^{(1)}=\hat{\beta}_{0,0}-\frac{1-\lambda}{\lambda} \hat{\beta}_{1,0} \\\tilde{y}_0^{(2)}=\hat{\beta}_{0,0}-2(\frac{1-\lambda}{\lambda})\hat{\beta}_{1,0} \\ \\\text{Forecasting:} \\e_{\tau}(1)=y_{T+1}-\hat{y}_{T+1}(T) \\\hat{y}_{T+\tau}(\tau)=(2\tilde{y}_T^{(1)}-\tilde{y}_T^{(2)})+\tau\frac{\lambda}{1-\lambda}(\tilde{y}_T^{(1)}-\tilde{y}_T^{(2)}) \\\hat{y}_{T+\tau}(\tau)=(2+\frac{\lambda}{1-\lambda}\tau)\tilde{y}_T^{(1)}-(1+\frac{\lambda}{1-\lambda}\tau)\tilde{y}_T^{(2)} \]

Data Stasioner

Data Stasioner merupakan data deret waktu yang memiliki ciri-ciri statistik dasar seperti nilai tengah dan ragam yang tetap (konstan) sepanjang waktu. Menurut Pramujo et. al (2014), keStasioneran data deret waktu memiliki arti bahwa dalam data tersebut tidak terdapat pertumbuhan atau penurunan secara signifikan dan tidak memiliki kecenderungan terhadap trend tertentu. Dengan kata lain, fluktuasi data berada di sekitar suatu nilai rata-rata yang konstan, tidak bergantung pada waktu maupun ragam dari fluktuasi tersebut. Pada dasarnya, data Stasioner dibagi menjadi dua, yaitu Stasioner dalam mean (rata-rata) dan Stasioner dalam variance (ragam). Apabila data Stasioner dalam rata-rata, maka terjadi fluktuasi data di sekitar suatu nilai rata-rata yang konstan. Sedangkan data Stasioner dalam ragam, yaitu data yang dari waktu ke waktu mempunyai fluktuasi yang tetap atau konstan dan tidak berubah-ubah (Pramujo et. al 2014).

Model Autoregressive Integrated Moving Average (ARIMA)

Metode Autoregressive Integrated Moving Average (ARIMA) adalah model peramalan yang sering digunakan dalam pemodelan time series dan menghasilkan ramalan-ramalan berdasarkan sintesis dari pola data secara historis. Metode ARIMA sering digunakan untuk peramalan di berbagai bidang, salah satunya bidang transportasi (Nurjanah, et al. 2018). Model peramalan univariate time series yang memadukan antara model Autoregressive (AR), model Moving Average (MA) dan model campuran autoregressive moving average (ARMA) untuk data yang tidak stasioner. Bentuk umum ARIMA (p, d, q) sebagai berikut:

\[ \phi_p(B)Y_t^*=\theta_0+\theta_q(B)e_t\\ \] \[ \begin{aligned} \text{Dengan:}\\ Y_t^*&= \text{ data deret waktu yang telah dilakukan differencing}\\ \phi_p(B)&=1-\phi_1B-\phi_2B^2-...-\phi_pB^p\\ \phi_q(B)&=1-\phi_1B-\phi_2B^2-...-\phi_qB^q\\ (1-B)^d&=\text{ differencing non-musiman pada orde ke-d}\\ e_t&=\text{ nilai kesalahan pada waktu ke-t} \end{aligned} \] Di mana nilai \(d>0\), untuk \(d=0\) model stasioner sama dengan \(ARMA(p,q)\).

Pembedaan (differencing) dengan lag-d dilakukan untuk menstasionerkan data. Konsep differencing adalah mengurangkan antara pengamatan Yt dengan pengamatan sebelumnya yaitu Yt-1. Metode estimasi maximum likelihood masih dianggap metode estimasi parameter ARIMA yang efisien (As’ad, et al. 2017).

Augmented Dickey-Fuller (ADF)

Uji Augmented Dickey Fuller (ADF) merupakan pengujian yang dilakukan terhadap data deret waktu (time series) untuk mengetahui apakah data deret waktu tersebut stationer atau tidak. Tidak semua proses runtun waktu dapat direpresentasikan dengan baik menggunakan model AR(1). Jika model merupakan suatu autoregresif tingkat p, dengan p ≥ 1, Dickey dan Fuller menambahkan tiga statistik-F untuk uji hipotesis gabungan (joint hypotesis) pada koefisien-koefisien model autoregresif yang terbentuk. Uji akar-akar unit metode Dickey-Fuller untuk model autoregresif tingkat p dengan p ≥ 1 dikenal sebagai uji Dickey-Fuller diperluas (Augmented Dickey-Fuller Test) (Rusdi, 2011).

Teknik Interpolasi

Teknik interpolasi memperkirakan nilai yang hilang dengan mengasumsikan hubungan dalam rentang titik data. Interpolasi adalah perkiraan suatu nilai tengah dari suatu set nilai yang diketahui. Dalam teknik statistik bergulir di mana hanya nilai-nilai sebelumnya yang dianggap memperhitungkan nilai yang hilang, teknik interpolasi memperkirakan dengan titik data masa lalu dan masa akan datang yang diketahui. Interpolasi memiliki berbagai metode:

• Linier: Mengasumsikan hubungan linier b/w rentang titik data

• Spline: Memperkirakan nilai yang meminimalkan kelengkungan keseluruhan sehingga memperoleh permukaan halus yang melewati titik input.

• Waktu: Memperkirakan nilai yang hilang dengan lebih fokus pada titik terdekat daripada titik jauh. Ada berbagai teknik interpolasi lain yang dapat digunakan untuk menghitung nilai yang hilang.

Semua ini mengasumsikan bahwa titik data yang berdekatan serupa, yang tidak selalu demikian. Menurut Fatimah (2015), berdasarkan derajat polinomial interpolasi dibagi menjadi empat macam yaitu:

1. Interpolasi Linear merupakan interpolasi untuk polinomial berderajat satu.

2. Interpolasi Kuadrat merupakan interpolasi untuk polinomial berderajat dua.

3. Interpolasi Lagrange merupakan interpolasi untuk polinomial berderajat m-1 jika diberikan n titik, yaitu 𝑥1, 𝑥2, …, 𝑥𝑛.

4. Interpolasi Newton merupakan interpolasi untuk polinomial berderajat m jika diberikan n+1 titik, yaitu 𝑥0, 𝑥1, 𝑥2,… , 𝑥𝑛.

Differencing/Pembedaan

Teknik Differencing (Pembedaan) dapat membantu menstabilkan nilai rata-rata deret waktu dengan menghilangkan perubahan tingkat deret waktu, dan karenanya menghilangkan (atau mengurangi) tren dan musiman.

Import Data

Data yang digunakan hanyalah data Close, didapatkan dengan syntax berikut:

gs4_deauth()
IHSG_2Yr <- read_sheet("https://docs.google.com/spreadsheets/d/1OIC65ByfmmT9SX7I0URCxf3cH24vpj5xVO-YwM6a3SQ/edit?usp=sharing")
IHSG_2Yr <- IHSG_2Yr[c("Date","Close")]
IHSG_2Yr

# Cek apakah ada missing value:
IHSG_2Yr[which(is.na(IHSG_2Yr$Close)),]

Dari output di atas, dapat terlihat bahwa terdapat missing value pada baris ke-87 atau pada tanggal 2 mei 2022. Hal ini melanggar salah satu syarat data time series, yaitu rentang waktu yang tidak sama. Oleh karena itu, akan dilakukan penanganan missing value.

Menduga Missing Value

Pendugaan dilakukan dengan teknik interpolasi “spline”. Nilai dugaan dari teknik ini didapatkan dengan memperhatikan pola nilai sebelum dan setelah nilai missing value. Hal ini lebih baik daripada metode lain yang hanya mengambil nilai rataan 2 nilai sebelum dan setelah, atau bahkan hanya menduplikasi nilai sebelum/setelahnya. Oleh karena itu, teknik ini dianggap paling akurat untuk menduga nilai missing value. Fungsi yang digunakan untuk mengatasi missing value adalah fungsi na_interpolation yang tersedia pada package imputeTS, yang syntax-nya adalah sebagai berikut:

IHSG_2Yr <- na_interpolation(IHSG_2Yr, option = "spline")
IHSG_2Yr[81:90,]

Setelah dilakukan penanganan, dapat terlihat bahwa nilai IHSG pada tanggal 2 mei 2022 telah terisi dengan nilai \(6840.217\).

Eksplorasi Data

Eksplorasi data diperlukan untuk melihat pola dari data time series. Selain itu, eksplorasi juga berguna dalam pembagian data training dan data testing. Data testing yang baik adalah data testing yang memiliki pola seperti data trainingnya. Maka dari itu, diperlukan eksplorasi data untuk menentukan proporsi terbaik data training vs data testing.

IHSG_2Yr$Date <- as.Date(IHSG_2Yr$Date)
ggplot(IHSG_2Yr, aes(x=Date, y=Close)) +
  geom_line() +
  scale_x_date(date_breaks = "3 month", date_labels = "%Y %b") +
  labs(title = "Plot Time Series Data IHSG mingguan",
       subtitle = "(September 2020 sd September 2022)",
       y="Harga IHSG") +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5))

Dari plot di atas, terlihat bahwa data cenderung memiliki pola trend yang naik. Terlihat pula bahwa pembagian data training-testing dapat dilakukan antara bulan april 2022 hingga juli 2022. Hal ini karena pola setelah bulan-bulan tersebut cenderung memiliki pola trend naik pula, sama seperti data trainingnya. Untuk lebih jelasnya, perhatikan penjelasan di bawah ini:

p = 0.81 #setelah dilakukan pencarian nilai proporsi data training terbaik
freq_train=as.integer(p*nrow(IHSG_2Yr))
IHSG_2Yr$Date <- as.Date(IHSG_2Yr$Date)
ggplot(IHSG_2Yr, aes(x=Date, y=Close)) +
  geom_line() + 
  scale_x_date(date_breaks = "3 month", date_labels = "%Y %b") +
  labs(title = "Plot Time Series Data IHSG mingguan",
       subtitle = "(September 2020 sd September 2022)") +
  geom_vline(aes(xintercept = Date[freq_train], 
                 col="Frequency_Train_Data"), lty=2, lwd=.7) +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5),
        legend.position = "bottom") +
  scale_color_manual(name = "", values = c(Frequency_Train_Data="red"))

Setelah dilakukan eksplorasi kembali, didapatkan titik potong yang sekiranya baik dalam membagi data training testing terdapat pada nilai proporsi data training 0.81 atau \(p=81\%\). Titik tersebut bertepatan pada tanggal 2 Mei 2022 (didapatkan dengan syntax di bawah).

IHSG_2Yr[freq_train,1]

Selanjutnya, akan dilakukan pembagian data training testing berdasarkan eksplorasi di atas.

Splitting Data

train <- 1:freq_train
y <- IHSG_2Yr$Close
y_train <- ts(y[train])
y_test <- ts(y[-train], start = 88)
ggplot(IHSG_2Yr[train,], aes(x=Date, y=Close)) +
  geom_line() +
  scale_x_date(date_breaks = "3 month", date_labels = "%Y %b") +
  labs(title = "Plot Time Series Data Training",
       subtitle = "(September 2020 sd Mei 2022)",
       y="Harga IHSG") +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5))

ggplot(IHSG_2Yr[-train,], aes(x=Date, y=Close)) +
  geom_line() +
  scale_x_date(date_breaks = "1 month", date_labels = "%Y %b") +
  labs(title = "Plot Time Series Data Testing",
       subtitle = "(Mei 2022 sd September 2022)",
       y="Harga IHSG") +
  theme(plot.title =  element_text(face = "bold", hjust=.5),
        plot.subtitle = element_text(hjust=.5))

Terlihat bahwa pola data testing tidak berbeda jauh dengan pola data training. Hal ini menandakan bahwa pembagian kedua data ini telah tepat.

Double Moving Average (DMA)

Metode pemulusan Double Moving Average (DMA) merupakan metode yang tepat untuk data IHSG Close karena data tersebut memiliki pola data tren. Fungsi di bawah ini akan me-looping nilai n dari 3 sampai dengan 11 sehingga didapatkan n optimum.

akurasi = data.frame()
for(n in 3:11){
  data.sma<-SMA(y_train, n=n)
  dma <- SMA(data.sma, n=n)
  At <- 2*data.sma - dma
  Bt <- 2/(n-1)*(data.sma - dma)
  data.ramal<- c(NA, At+Bt)
  
  error.dma = y_train-data.ramal[1:length(y_train)]
  SSE.dma = sum(error.dma[(2*n):length(y_train)]^2, na.rm = T)
  RMSE.dma = sqrt(mean(error.dma[(2*n):length(y_train)]^2, na.rm = T))
  MAPE.dma = mean(abs((error.dma[(2*n):length(y_train)]/y_train[(2*n):length(y_train)])*100), na.rm = T)
  
  akurasi.dma <- matrix(c(SSE.dma, RMSE.dma, MAPE.dma))
  row.names(akurasi.dma)<- c("SSE", "RMSE", "MAPE")
  colnames(akurasi.dma) <- c(paste("Akurasi m =", n))
  akurasi = rbind(akurasi, c(n=n, RMSE=akurasi.dma[2], MAPE=akurasi.dma[3]))
}
colnames(akurasi) = c("n", "RMSE", "MAPE")
akurasi[order(akurasi[,2]),]

Berdasarkan tabel di atas, nilai RMSE dan MAPE terkecil dimiliki oleh metode pemulusan DMA pada n = 4 yang dapat dijadikan n optimum dalam rentang nilai n dari 3 sampai 11. Dengan nilai MSE terkecil sebesar \(143.7413\) dan MAPE terkecil sebesar \(1.664693\). Syntax di bawah ini akan memodelkan pemulusan DMA dengan n optimum.

n=4
data.sma<-SMA(y_train, n=n)
dma <- SMA(data.sma, n=n)
At <- 2*data.sma - dma
Bt <- 2/(n-1)*(data.sma - dma)
data.ramal<- c(NA, At+Bt)

fcast = c()
for (i in 1:21) {
  fcast[i] = At[length(At)] + Bt[length(Bt)]*(i)
} # Ramalan DMA

data.gab <- cbind(aktual = c(y_train,rep(NA,21)), pemulusan = c(dma, rep(NA,21)), At = c(At, rep(NA,21)), Bt = c(Bt,rep(NA,21)), ramalan = c(data.ramal, fcast[-1]))
tail(data.gab)

##        aktual pemulusan At Bt  ramalan
## [103,]     NA        NA NA NA 6703.473
## [104,]     NA        NA NA NA 6678.953
## [105,]     NA        NA NA NA 6654.433
## [106,]     NA        NA NA NA 6629.913
## [107,]     NA        NA NA NA 6605.394
## [108,]     NA        NA NA NA 6580.874

# Plot
plot.ts(IHSG_2Yr$Close, xlab="Weekly", ylab="Close", main="Plot Pemulusan dan Peramalan dengan Metode DMA N Optimum (4)", ylim=c(5000,7500), type="l")
points(data.gab[,1])
lines(data.gab[,2],col="green",lwd=2)
lines(data.gab[,5],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan DMA","data peramalan"), lty=1, col=c("black","green","red"), cex=0.8)

RMSE_DMA_train <- sqrt(mean((y_train-data.ramal[-88])^2, na.rm = T))
RMSE_DMA_test <- sqrt(mean((y_test-fcast)^2))
MAPE_DMA_train <- mean(abs((y_train-data.ramal[-88])/y_train)*100, na.rm = T)
MAPE_DMA_test <- mean(abs((y_test-fcast)/y_test)*100)
Akurasi_DMA <- data.frame("DMA_Train"=c(RMSE_DMA_train,MAPE_DMA_train), "DMA_Test"=c(RMSE_DMA_test,MAPE_DMA_test))
rownames(Akurasi_DMA) <- c("RMSE","MAPE")
Akurasi_DMA

Double Exponential Smoothing (DES)

Selain DMA, Metode pemulusan lain yang cocok untuk data yang memiliki pola data tren adalah DES. Metode ini memandang bahwa semakin besar lag, maka semakin tidak berpengaruh data tersebut terhadap data aktual saat ini. Fungsi di bawah ini akan me-looping parameter DES, yaitu alpha dan beta (dari 0.1 hingga 0.9) sehingga didapatkan parameter optimum.

#Pengulangan Melihat Nilai Keakuratan
a = seq(.1,.9,.1)
b = seq(.1, .9,.1)
output = data.frame()
for (i in a) {
  for (j in b) {
    des <- HoltWinters(y_train, alpha = i, beta=j, gamma=F)
    sse <- des$SSE
    rmse <- sqrt(sse/length(y_train))
    mape <- mean(abs((des$fitted[,1]-y_train)/y_train)*100)
    akurasi <- cbind("SSE"=sse, "RMSE"=rmse, "MAPE"=mape)
    output <- rbind(output, akurasi)
  }
}
output <- cbind("Alpha" = rep(a,each=9), "Beta"=b, output)
dplyr::arrange(.data=output,RMSE) #Mengurutkan data dari nilai SSE/MSE terkecil

Setelah diurutkan berdasarkan nilai RMSE, dapat terlihat bahwa nilai alpha optimum 0.9 dan beta optimum 0.1, dengan nilai RMSE terkecilnya yaitu sebesar \(1410.87\). Kedua parameter ini yang kemudian akan digunakan untuk pemulusan.

df_des <- HoltWinters(y_train, alpha = .9, beta=.1, gamma=F)
datades <- data.frame(y_train, c(NA, NA, df_des$fitted[,1]))
colnames(datades) = c("y_train","y_train_hat")
datades

ramal_des <- forecast::forecast(df_des,h=21)
data.frame(ramal_des)

plot(ramal_des, xlab="periode waktu", ylab="Yt",
     main="Forecasting Using Double Exponential Smoothing (0.9, 0.1)")
lines(y_test, col = "black", lwd=1.8)
lines (datades[,2], col="green",lwd=1)
legend("topleft", legend = c("Data aktual", "Fitted DES", "Ramalan"), lty=1, col= c("black", "green", "blue"))

RMSE_DES_train <- sqrt(df_des$SSE/length(y_train))
RMSE_DES_test <- sqrt(mean((y_test-ramal_des$mean)^2))
MAPE_DES_train <- mean(abs((df_des$fitted[,1]-y_train)/y_train)*100)
MAPE_DES_test <- mean(abs((y_test-ramal_des$mean)/y_test)*100)
Akurasi_DES <- data.frame("DES_Train"=c(RMSE_DES_train,MAPE_DES_train), "DES_Test"=c(RMSE_DES_test,MAPE_DES_test))
rownames(Akurasi_DES) <- c("RMSE","MAPE")
Akurasi_DES

Perbandingan DMA vs DES

plot(ramal_des, xlab="periode waktu", ylab="Yt",
     main="Forecasting DMA vs DES")
lines(y_test, col = "black", lwd=1.8)
lines (datades[,2], col="green",lwd=1.5)
lines(data.gab[,5],col="red",lwd=1.5)
legend("topleft", legend = c("Aktual", "Fitted DES", "Ramalan DES", "Ramalan DMA"), lty=1, col= c("black", "green", "blue", "red"))

Akurasi <- data.frame("DMA"=Akurasi_DMA[,2], "DES"=Akurasi_DES[,2])
rownames(Akurasi) <- c("RMSE","MAPE")
Akurasi

Prosedur Eksploratif

Secara eksploratif, kestasioneran data dapat terlihat dari plot deret waktu atau melihat plot acf dari peubah yang akan diuji. Jika pada plot deret waktu menunjukan pola tren naik/turun pada kurun waktu tertentu (tidak stabil), atau, jika plot acf menunjukkan pola tails off slowly, maka data tersebut tidak stasioner.

plot.ts(y_train, ylab=expression(Y[t]), main = "Plot Time Series Data IHSG mingguan")

acf(y_train)

Dapat terlihat dari plot deret waktu dan plot acf peubah Y_train di atas, bahwa secara eksploratif dapat disimpulkan data tersebut tidak stasioner. Uji ADF akan dilakukan untuk meyakinkan kesimpulan tersebut.

Uji Formal

Secara formal, metode Augmented Dickey-Fuller (ADF) dapat memberikan hasil uji secara akurat untuk menentukan apakah sebuah data stasioner atau tidak. Dengan hipotesis sebagai berikut:

\[ \begin{aligned} H_0&: \text{ Data tidak stasioner}\\ H_1&: \text{ Data stasioner} \end{aligned} \]

adf.test(y_train)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  y_train
## Dickey-Fuller = -2.8143, Lag order = 4, p-value = 0.2414
## alternative hypothesis: stationary

Didapat nilai \(P-value>0.05\) maka tak tolak \(H_0\) pada taraf nyata 5%. Artinya, cukup bukti untuk mengatakan bahwa data tidak stasioner. Untuk mengatasi ketidakstasioneran ini, perlu dilakukan differencing atau pembedaan.

Cek ulang kestasioneran data

adf.test(y_diff)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  y_diff
## Dickey-Fuller = -3.1144, Lag order = 4, p-value = 0.1183
## alternative hypothesis: stationary

Didapat nilai \(P-value=0.1183>0.05\) maka tak tolak \(H_0\) pada taraf nyata 5% atau cukup bukti untuk mengatakan bahwa data masih tidak stasioner. Artinya, differencing pertama gagal untuk membuat data stasioner, sehingga perlu dilakukan differencing kedua untuk mengatasi ketidakstasioneran ini.

Differencing #2

y_diff2 <- diff(y_train, differences = 2)
y_diff2

## Time Series:
## Start = 3 
## End = 87 
## Frequency = 1 
##  [1] -155.943360   94.375000  145.986328  -77.178222  -40.977052    7.263185
##  [7]  191.266601  -81.774415  -14.931639  101.081542 -184.531250  100.698243
## [13]   38.148925 -261.610352   65.979493  308.397461 -163.184571 -181.862304
## [19] -378.489746  734.151855 -218.584961  -61.380859    0.452636    7.089356
## [25]   82.506836 -101.508789 -158.549316  -23.507813  242.858886  -42.704100
## [31]  -85.442872   48.146484  -46.059570   77.349121 -175.272949  240.727051
## [37]  141.053711 -186.218751 -118.708007  103.655762  -14.669922   16.227538
## [43]   15.829103   -3.485353  -60.831053  165.042967 -197.331053  -44.781251
## [49]  119.314453   74.960449 -117.602538   70.420898  -26.804200   72.461427
## [55]  168.893554 -101.354980 -141.168946  -62.791504   42.830567   78.830078
## [61]   -0.060058 -227.918946  135.662597  137.463380 -165.405762   11.957519
## [67]   57.614259  101.251952 -127.749023   40.887207 -113.833985  166.741700
## [73]   -1.664062   -7.004884  -81.857910   44.804200  -45.883301   38.088867
## [79]   15.204590   28.660156   55.847656 -107.377929  -34.623536   13.234376
## [85] -392.005536

Cek ulang kestasioneran data #2

adf.test(y_diff2)

## Warning in adf.test(y_diff2): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  y_diff2
## Dickey-Fuller = -5.6831, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Dengan implementasi yang sama seperti bagian sebelumnya, hasil di atas menunjukkan bahwa data setelah differencing telah stasioner.

ACF

acf(y_diff2)

Dari plot ACF di atas, dapat terlihat bahwa terdapat cuts off setelah lag 2. Oleh karena itu, model yang terbentuk adalah ARIMA(0,2,2).

PACF

pacf(y_diff2)

Dari plot PACF di atas, dapat terlihat bahwa terdapat cuts off setelah lag 3. Oleh karena itu, model yang terbentuk adalah ARIMA(3,2,0).

EACF

eacf(y_diff2)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 x o o o o o o o o o o  o  o  o 
## 4 o x o x o o o o o o o  o  o  o 
## 5 x o o o o o o o o o o  o  o  o 
## 6 x x o x o o o o o o o  o  o  o 
## 7 x x o x o o o o o o o  o  o  o

Dari tabel EACF di atas, dapat terlihat bahwa terdapat triangle of zero dengan titik kiri atasnya yaitu pada 0,2. Artinya, model yang terbentuk Oleh karena itu, model yang terbentuk menurut EACF adalah ARIMA(0,1,2). Namun, model ini telah didapat pada plot acf, maka triangle of zero selanjutnya (1,1) dan (2,2) akan dibentuk model. Jadi, model tambahan yang terbentuk adalah ARIMA(1,2,1) dan ARIMA(2,2,2).

Kesimpulan

Berdasarkan hasil eksploratif di atas, dapat disimpulkan bahwa terdapat tiga model yang menjadi dugaan model terbaik untuk data IHSG, yaitu ARIMA(0,2,2), ARIMA(3,2,0). ARIMA(1,2,1)., ARIMA(2,2,2). Keempat model akan dianalisis lebih lanjut pada bagian selanjutnya.

Signifikansi Parameter

arima022 <- Arima(y_train, order=c(0,2,2), method="ML")
coeftest(arima022)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.85258    0.15539 -5.4867 4.095e-08 ***
## ma2 -0.14742    0.14697 -1.0031    0.3158    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

arima320 <- Arima(y_train, order=c(3,2,0), method="ML")
coeftest(arima320)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.66173    0.10517 -6.2918 3.137e-10 ***
## ar2 -0.63274    0.10575 -5.9835 2.184e-09 ***
## ar3 -0.41186    0.10404 -3.9587 7.537e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

arima121 <- Arima(y_train, order=c(1,2,1), method="ML")
coeftest(arima121)

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ar1  0.084953   0.117483   0.7231   0.4696    
## ma1 -0.999999   0.053868 -18.5638   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

arima222 <- Arima(y_train, order=c(2,2,2), method="ML")
summary(arima222)

## Series: y_train 
## ARIMA(2,2,2) 
## 
## Coefficients:
##          ar1      ar2      ma1      ma2
##       0.0625  -0.2276  -0.9629  -0.0365
## s.e.  0.3198   0.1181   0.2504   0.3180
## 
## sigma^2 = 13306:  log likelihood = -524.51
## AIC=1059.02   AICc=1059.78   BIC=1071.23
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -12.2596 111.3019 75.13089 -0.1961632 1.221884 0.9212663
##                     ACF1
## Training set -0.01354769

coeftest(arima222)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.062542   0.319819  0.1956 0.8449603    
## ar2 -0.227617   0.118073 -1.9278 0.0538835 .  
## ma1 -0.962874   0.250448 -3.8446 0.0001207 ***
## ma2 -0.036513   0.318020 -0.1148 0.9085925    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Keakuratan Model

akurasi_model = data.frame("Model" = c("ARIMA(0,2,2)", "ARIMA(3,2,0)", "ARIMA(1,2,1)", "ARIMA(2,2,2)"),
                        "AIC" = c(arima022$aic, arima320$aic, arima121$aic,arima222$aic),
                       "BIC" = c(arima022$bic, arima320$bic, arima121$bic,arima222$bic))
akurasi_model[order(akurasi_model[,3]),]

Dari ukuran AIC dan BIC, terlihat bahwa model yang memiliki nilai AIC dan BIC terkecil adalah ARIMA(0,2,2), yaitu sebesar 1058.542 dan 1065.87. Kemudian, model yang memiliki nilai AIC dan BIC terbesar adalah ARIMA(3,2,0).

Dugaan Model terbaik

Berdasarkan hasil di atas, didapatkan tabel sebagai berikut:

Jika melihat tabel di atas, maka terdapat 2 kesimpulan:

1. Jika dilihat dari ukuran parameter model, maka model terbaiknya adalah ARIMA(3,2,0).
2. Jika dilihat dari nilai keakuratan (AIC), maka model terbaiknya adalah ARIMA(0,2,2)

Oleh karena itu, maka akan dilakukan tes lain terhadap kedua model ini yaitu perbandingan keakuratan ramalan terhadap data testing.

Perbandingan Ramalan Dua Model Terbaik

ramalan022 <- forecast(arima022, h = 21)
ramalan320 <- forecast(arima320, h = 21)
plot(ramalan022); lines(y_test)

plot(ramalan320); lines(y_test)

akurasi022 <- accuracy(ramalan022, y_test)
akurasi320 <- accuracy(ramalan320, y_test)

akurasi_model = data.frame("Model" = c("ARIMA(0,2,2)", "ARIMA(3,2,0)"),
                        "RMSE" = c(akurasi022[2,2], akurasi320[2,2]),
                       "MAPE" = c(akurasi022[2,5], akurasi320[2,5]))
akurasi_model[order(akurasi_model[,3]),]

Kesimpulan Model terbaik

Berdasarkan hasil yang telah didapat di atas, dapat disimpulkan bahwa model terbaik untuk melakukan ramalan data IHSG mingguan adalah ARIMA(0,2,2).

ARIMA(1,2,2)

arima122 <- Arima(y_train, order=c(1,2,2), method="ML")
summary(arima122)

## Series: y_train 
## ARIMA(1,2,2) 
## 
## Coefficients:
##           ar1      ma1      ma2
##       -0.3882  -0.4615  -0.5384
## s.e.   0.3733   0.3359   0.3327
## 
## sigma^2 = 13614:  log likelihood = -525.77
## AIC=1059.54   AICc=1060.04   BIC=1069.31
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -10.12641 113.2752 76.53896 -0.1609594 1.242867 0.9385323
##                     ACF1
## Training set -0.04327006

coeftest(arima122)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)
## ar1 -0.38818    0.37335 -1.0397   0.2985
## ma1 -0.46150    0.33590 -1.3739   0.1695
## ma2 -0.53841    0.33272 -1.6182   0.1056

Dari model overfitting ARIMA(1,2,2), dapat terlihat bahwa tidak ada parameter yang signifikan. Dapat terlihat jelas bahwa model ini tidak lebih baik daripada model awal.

ARIMA(0,2,3)

arima023 <- Arima(y_train, order=c(0,2,3))
summary(arima023)

## Series: y_train 
## ARIMA(0,2,3) 
## 
## Coefficients:
##           ma1      ma2     ma3
##       -0.9039  -0.2574  0.1613
## s.e.   0.1559   0.1389  0.1026
## 
## sigma^2 = 13329:  log likelihood = -525.06
## AIC=1058.11   AICc=1058.61   BIC=1067.88
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -11.58985 112.0856 75.74979 -0.1852685 1.231119 0.9288554
##                      ACF1
## Training set -0.009719192

coeftest(arima023)

## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.90389    0.15590 -5.7978 6.718e-09 ***
## ma2 -0.25738    0.13892 -1.8527   0.06392 .  
## ma3  0.16131    0.10255  1.5730   0.11572    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dari model overfitting ARIMA(0,2,3), dapat terlihat bahwa parameter ma3 yang ditambahkan tidak signifikan. Artinya, model ini juga tidak lebih baik daripada model awal.

Nilai AIC dan BIC Model

akurasi_ovrfit = data.frame("Model" = c("ARIMA(0,2,2)", "ARIMA(1,2,2)", "ARIMA(0,2,3)"),
                        "AIC" = c(arima022$aic, arima122$aic, arima023$aic),
                       "BIC" = c(arima022$bic, arima122$bic, arima023$bic))
akurasi_ovrfit[order(akurasi_ovrfit[,3]),]

Dari ukuran AIC, terlihat bahwa model yang memiliki nilai AIC terkecil adalah ARIMA(0,2,3), yaitu sebesar 1058.112. Sedangkan, model yang memiliki nilai BIC terkecil adalah ARIMA(0,2,2).

Kesimpulan

Dari signifikansi model terbaik dan ukuran BIC, dapat disimpulkan bahwa model terbaiknya adalah model awal yaitu ARIMA(0,2,2).

Prosedur Eskploratif

#Eksplorasi
best_model <- arima022
sisaan <- best_model$residuals
par(mfrow=c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "blue", lwd = 2)
plot(as.numeric(best_model$fitted), as.numeric(sisaan), xlab="Fitted", ylab="Sisaan", 
     ylim=c(-500,500), main="Residuals vs Fitted"); abline(h=0, col="red")
acf(sisaan)
plot(1:length(sisaan),sisaan,type='o', xlab="Order", ylim=c(-500,500),
     main="Residuals vs Order"); abline(h=0, col="red")

1. Normal Q-Q Plot
   Berdasarkan hasil eksplorasi di atas, terlihat bahwa banyak amatan sisaan cenderung menjauhi garis qq-plot distribusi normal. Oleh karena itu, secara eksploratif dapat disimpulkan bahwa sisaan belum cukup menyebar normal.

2. Plot Residual vs Fitted
   Berdasarkan hasil eksplorasi di atas, mayoritas titik amatan berada disekitar nol meskipun ada beberapa titik amatan yang terletak jauh dari titik nol. Selain itu, dapat terlihat bahwa lebar pita dari sisaan tidak terlalu berbeda jauh dari kiri hingga ke kanan, hal ini mengindikasikan bahwa ragam dari sisaan cukup homogen.

3. Plot ACF
   Berdasarkan hasil eksplorasi di atas, pada plot ACF tidak terdapat garis vertikal di lag tertentu yang melebihi tinggi garis biru horizontal.Oleh karena itu, secara eksploratif plot ACF dapat disimpulkan tidak terdapat autokorelasi pada model.

4. Residual vs Order
   Berdasarkan hasil eksplorasi di atas, titik amatan pada plot kebebasan sisaan mayoritas bergerak di sekitar titik nol. Namun, terdapat beberapa titik amatan yang terletak cukup jauh dari titik nol. Oleh karena itu, secara eksploratif plot residual vs order belum dapat disimpulkan apakah sisaan tersebut terdapat autokorelasi atau tidak.

Uji Formal

ks.test(sisaan,"pnorm")

## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.51689, p-value < 2.2e-16
## alternative hypothesis: two-sided

Box.test(sisaan,type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 0.11131, df = 1, p-value = 0.7387

t.test(sisaan, mu = 0, conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  sisaan
## t = -0.80804, df = 86, p-value = 0.4213
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -34.23921  14.44890
## sample estimates:
## mean of x 
## -9.895154

Plot

ramalan_best <- forecast(best_model, h = 21)
plot(ramalan_best); lines(y_test)

Dari plot di atas, dapat terlihat bahwa model ARIMA(0,2,2) sudah sangat baik dalam meramalkan nilai IHSG 21 pekan ke depan. Hal ini terlihat dari miripnya ramalan dari model dengan data aktual nilai IHSG-nya.

Nilai Akurasi

accuracy(ramalan_best, y_test)

##                     ME     RMSE       MAE         MPE     MAPE      MASE
## Training set -9.895154 113.9941  77.20645 -0.15716217 1.253738 0.9467171
## Test set      7.032481 150.1718 117.19661  0.05742668 1.691531 1.4370825
##                    ACF1 Theil's U
## Training set -0.0351606        NA
## Test set      0.5995693  1.101872

Dari nilai keakuratan di atas, terlihat bahwa tidak terjadi overfitting atau underfitting pada model. Hal ini karena nilai akurasi pada data testing tidak terlalu berbeda jauh jika dibandingkan data training. Misalnya, nilai MAPE pada data training sebesar 1.25 sedangkan pada data testing sebesar 1.69. Hanya terdapat selisih 0.44 nilai MAPE antara training testing atau MAPE data training lebih tinggi sekitar 35% dari MAPE data training. Perbedaan kedua nilai ini tidaklah terlalu signifikan.

Plot

plot.ts(IHSG_2Yr$Close,xlab="Weekly", ylab="Close", main="Plot Perbandingan Pemulusan dan Pemodelan", ylim=c(5000,7500), type="l")
points(y_train)
lines(ramal_des$mean, col="blue",lwd=2)
lines(ramalan_best$mean, col="red",lwd=2)
legend("topleft",c("data aktual","ramalan DES (0.9, 0.1)","ramalan ARIMA(0,2,2)"), lty=1, col=c("black","blue","red"), cex=0.8)

Dari plot di atas, dapat terlihat bahwa peramalan dengan menggunakan model ARIMA(0,2,2) jauh lebih mendekati plot data aktualnya bila dibandingkan dengan hasil peramalan metode DES.

Nilai Akurasi

akurasi_DES <- accuracy(ramal_des$mean, y_test)
akurasi_ARIMA <- accuracy(ramalan_best, y_test)
akurasi_Akhir <- rbind(akurasi_DES, akurasi_ARIMA[2,-6])
rownames(akurasi_Akhir) <- c("DES", "ARIMA(0,2,2)")
akurasi_Akhir

##                      ME     RMSE      MAE        MPE     MAPE      ACF1
## DES          189.695895 279.2619 246.5130 2.62875790 3.483038 0.6469238
## ARIMA(0,2,2)   7.032481 150.1718 117.1966 0.05742668 1.691531 0.5995693
##              Theil's U
## DES           2.030384
## ARIMA(0,2,2)  1.101872

Dari perbandingan nilai akurasi di atas juga terlihat bahwa ramalan pemodelan ARIMA(0,2,2) menghasilkan galat yang jauh lebih kecil jika dibandingkan galat yang dihasilkan oleh pemulusan DES. Sehingga dapat disimpulkan bahwa model ARIMA(0,2,2) merupakan metode terbaik untuk menduga harga penutupan IHSG mingguan berdasarkan harga 2 tahun terakhir.