Project MPDW : Peramalan Saham Unilever Menggunakan Pemodelan ARIMA ARCH-GARCH

Data yang kelompok 12 gunakan adalah data persentase saham Unilever di Indonesia 1 Januari 2018 - 30 Oktober 2023.

Library

library(ggplot2)
library(tsibble)
## Warning: package 'tsibble' was built under R version 4.3.3
## Registered S3 method overwritten by 'tsibble':
##   method               from 
##   as_tibble.grouped_df dplyr
## 
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(tseries)
## Warning: package 'tseries' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(MASS)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(TTR)
## Warning: package 'TTR' was built under R version 4.3.3
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.3.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.3.3
## 
## Attaching package: 'zoo'
## The following object is masked from 'package:tsibble':
## 
##     index
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(orcutt)
## Warning: package 'orcutt' was built under R version 4.3.3
library(HoRM)
## Warning: package 'HoRM' was built under R version 4.3.3
library(dLagM)
## Warning: package 'dLagM' was built under R version 4.3.3
## Loading required package: nardl
## Warning: package 'nardl' was built under R version 4.3.3
## Loading required package: dynlm
## Warning: package 'dynlm' was built under R version 4.3.3
## 
## Attaching package: 'dLagM'
## The following object is masked from 'package:forecast':
## 
##     forecast
library(dynlm)
library(MLmetrics)
## Warning: package 'MLmetrics' was built under R version 4.3.3
## 
## Attaching package: 'MLmetrics'
## The following object is masked from 'package:dLagM':
## 
##     MAPE
## The following object is masked from 'package:base':
## 
##     Recall
library(car)
## Warning: package 'car' was built under R version 4.3.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.3.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(TSA)
## Warning: package 'TSA' was built under R version 4.3.3
## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(aTSA)
## 
## Attaching package: 'aTSA'
## The following object is masked from 'package:dLagM':
## 
##     forecast
## The following object is masked from 'package:forecast':
## 
##     forecast
## The following objects are masked from 'package:tseries':
## 
##     adf.test, kpss.test, pp.test
## The following object is masked from 'package:graphics':
## 
##     identify
library(rugarch)
## Warning: package 'rugarch' was built under R version 4.3.3
## Loading required package: parallel
## 
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
## 
##     sigma
library(imputeTS)
## Warning: package 'imputeTS' was built under R version 4.3.3
## 
## Attaching package: 'imputeTS'
## The following object is masked from 'package:zoo':
## 
##     na.locf
## The following object is masked from 'package:tseries':
## 
##     na.remove
library(tidyr)

Impor Data

#install.packages("rio") #install jika belum ada
library(rio)
## Warning: package 'rio' was built under R version 4.3.3
datasaham <- import("https://raw.githubusercontent.com/aidara11/mpdw/main/Pertemuan%2011/Data%20Saham%20Unilever%202018-2023.csv")
head(datasaham)
##         Date  Open  High   Low Close Adj Close   Volume
## 1 2018-01-01 11180 11620 10635 10800  9254.643 44448000
## 2 2018-01-08 10800 10940 10700 10850  9297.488 43629500
## 3 2018-01-15 10850 11000 10850 10890  9331.764 59268000
## 4 2018-01-22 10890 11155 10860 10910  9348.902 44043000
## 5 2018-01-29 11000 11140 10500 11005  9430.310 59799500
## 6 2018-02-05 10790 11080 10700 11080  9494.579 75593000

Cek data

str(datasaham)         #struktur data
## 'data.frame':    305 obs. of  7 variables:
##  $ Date     : IDate, format: "2018-01-01" "2018-01-08" ...
##  $ Open     : num  11180 10800 10850 10890 11000 ...
##  $ High     : num  11620 10940 11000 11155 11140 ...
##  $ Low      : num  10635 10700 10850 10860 10500 ...
##  $ Close    : num  10800 10850 10890 10910 11005 ...
##  $ Adj Close: num  9255 9297 9332 9349 9430 ...
##  $ Volume   : int  44448000 43629500 59268000 44043000 59799500 75593000 35915500 50283000 56477500 71074000 ...
dim(datasaham)         #dimensi data
## [1] 305   7
sum(is.na(datasaham))  #cek data kosong
## [1] 0

Data persentase saham Unilever Indonesia tahun 2018 - 2023 terdiri dari 305 amatan. Peubah yang akan digunakan adalah peubah date dan open. Data saham Unilever sudah berformat numerik dan tidak terdapat data kosong.

Ubah menjadi time series

datasaham.ts <- ts(datasaham$Close)

Plot Data Penuh

plot.ts(datasaham.ts, lty=1, xlab="Periode", ylab="saham", main="Plot Persentase Saham Unilever di Indonesia tahun 2018-2023")

Berdasarkan plot data deret waktu, terlihat bahwa data cenderung memiliki trend yang turun. Berdasarkan pola data, pembagian data latih dan data uji ditetapkan dengan proporsi 80%:20%.

Plot Data Latih

sahamtrain <- datasaham.ts[1:244]
train.ts <- ts(sahamtrain)
plot.ts(train.ts, lty=1, xlab="Periode", ylab="saham", main="Plot Train Persentase Saham Unilever di Indonesia tahun 2018-2023")

Berdasarkan plot data deret waktu pada data latih, terlihat bahwa data cenderung memiliki trend yang turun dan cenderung tidak bergerak pada nilai tengah tertentu. Hal ini mengindikasikan bahwa data tidak stasioner dalam rataan.

Plot Data Uji

sahamtest <- datasaham.ts[245:305]
test.ts <- ts(sahamtest)
plot.ts(test.ts, lty=1, xlab="Periode", ylab="saham", main="Plot Test Persentase Saham Unilever di Indonesia tahun 2018-2023")

Berdasarkan plot data deret waktu pada data uji, terlihat bahwa data cenderung memiliki trend yang turun dan cenderung tidak bergerak pada nilai tengah tertentu. Hal ini mengindikasikan bahwa data uji tidak stasioner dalam rataan.

Cek pola data musiman

seasonplot(datasaham.ts,12, 
  main="Saham Unilever Indonesia 01/01/2018 - 30/10/2023", 
  xlab = "Bulan ke-",
  ylab = "Close Price",
  year.labels = TRUE, col=rainbow(18))

Plot data menunjukkan bahwa data tidak berpola musiman.

Uji Stasioneritas Data

Plot ACF

acf(train.ts)

Berdasarkan plot ACF, terlihat bahwa plot ACF data train menurun secara perlahan (tails of slowly) dan tidak membentuk gelombang sinus. Hal ini juga menjadi indikasi bahwa data tidak stasioner dalam rataan.

Uji ADF

tseries::adf.test(train.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.ts
## Dickey-Fuller = -2.2889, Lag order = 6, p-value = 0.4543
## alternative hypothesis: stationary

\(H_0\) : Data tidak stasioner dalam rataan \(H_1\) : Data stasioner dalam rataan

Berdasarkan uji ADF tersebut, didapat p-value sebesar 0.4543 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa data train tidak stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga ketidakstasioneran model kedepannya harus ditangani.

Plot Box-Cox

index <- seq(1:244)
bc = boxcox(train.ts~index, lambda = seq(-2, 4, by=.01))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda
## [1] 1.88
#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]
##  [1] 1.60 1.61 1.62 1.63 1.64 1.65 1.66 1.67 1.68 1.69 1.70 1.71 1.72 1.73 1.74
## [16] 1.75 1.76 1.77 1.78 1.79 1.80 1.81 1.82 1.83 1.84 1.85 1.86 1.87 1.88 1.89
## [31] 1.90 1.91 1.92 1.93 1.94 1.95 1.96 1.97 1.98 1.99 2.00 2.01 2.02 2.03 2.04
## [46] 2.05 2.06 2.07 2.08 2.09 2.10 2.11 2.12 2.13 2.14 2.15

Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 1.88 dan pada selang kepercayaan 95% nilai memiliki batas bawah 1.60 dan batas atas 2.15. Selang tersebut tidak memuat nilai satu sehingga dapat dikatakan bahwa data train tidak stasioner dalam ragam.

Penanganan Ketidakstasioneran Data

Penanganan ketidakstasioneran ragam

Transformasi Data

train.ts.new <- (train.ts)^2

Plot Box-Cox

index <- seq(1:244)
bc = boxcox(train.ts.new~index, lambda = seq(-2, 4, by=.01))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda
## [1] 0.94
#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]
##  [1] 0.80 0.81 0.82 0.83 0.84 0.85 0.86 0.87 0.88 0.89 0.90 0.91 0.92 0.93 0.94
## [16] 0.95 0.96 0.97 0.98 0.99 1.00 1.01 1.02 1.03 1.04 1.05 1.06 1.07
plot.ts(train.ts.new, lty=1, xlab="waktu", ylab="Data Transformasi", main="Plot Transformasi Data")

Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 0.94 dan pada selang kepercayaan 95% nilai memiliki batas bawah 0.80 dan batas atas 1.07. Selang tersebut memuat nilai satu sehingga dapat dikatakan bahwa data train stasioner dalam ragam.

Penanganan ketidakstasioneran rataan

train.diff <- diff(train.ts.new, differences = 1) 
plot.ts(train.diff, lty=1, xlab="Periode", ylab="Data Difference 1 saham", main="Plot Difference Persentase saham di Indonesia tahun 2003-2023")

Berdasarkan plot data deret waktu, terlihat bahwa data sudah stasioner dalam rataan ditandai dengan data bergerak pada nilai tengah tertentu (tidak terdapat trend ataupun musiman pada data).

Uji Stasioneritas Ulang

Plot ACF (rataan)

acf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot ACF cuts off pada lag ke 1. Hal ini menandakan data sudah stasioner dalam rataan dan ketidakstasioneran data telah berhasil tertangani.

Uji ADF (rataan)

tseries::adf.test(train.diff)
## Warning in tseries::adf.test(train.diff): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.diff
## Dickey-Fuller = -5.5483, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary

\(H_0\) : Data tidak stasioner dalam rataan

\(H_1\) : Data stasioner dalam rataan

Berdasarkan uji ADF tersebut, didapat p-value sebesar 0.01 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) atau data stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga dalam hal ini ketidakstasioneran data sudah berhasil ditangani dan dapat dilanjutkan ke pemodelan.

Identifikasi Model

Plot ACF

acf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot ACF cenderung cuts off pada lag ke 1, sehingga jika plot PACF dianggap tails of, maka model tentatifnya adalah ARIMA(0,1,1).

Plot PACF

pacf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot PACF cenderung cuts off pada lag ke 1, sehingga jika plot ACF dianggap tails of, maka model tentatifnya adalah ARIMA(1,1,0).

Plot EACF

eacf(train.diff)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o x x o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x o o o o o o x o o o  o  o  o 
## 3 x x x o o o o x o o o  o  o  o 
## 4 o x x o o o o x o o o  o  o  o 
## 5 x x x o o o o x o o o  o  o  o 
## 6 x o x o x o o o o o o  o  o  o 
## 7 x o x o x x x o o o o  o  o  o

Identifikasi model menggunakan plot EACF dilakukan dengan melihat ujung segitiga pada pola segitiga nol. Dalam hal ini model tentatif yang terbentuk adalah ARIMA(0,1,0), ARIMA(0,1,1), ARIMA(1,1,1), ARIMA(1,1,2), ARIMA(2,1,2), ARIMA(2,1,3).

Pendugaan Parameter Model Tentatif

ARIMA(0,1,0)

model1.da=Arima(train.diff, order=c(0,1,0),method="ML")
summary(model1.da) 
## Series: train.diff 
## ARIMA(0,1,0) 
## 
## sigma^2 = 4.454e+13:  log likelihood = -4146.11
## AIC=8294.21   AICc=8294.23   BIC=8297.7
## 
## Training set error measures:
##                     ME    RMSE     MAE MPE MAPE      MASE       ACF1
## Training set -6712.006 6660242 4715314 NaN  Inf 0.9958857 -0.5663285
#lmtest::coeftest(model1.da)

Didapatkan nilai AIC sebesar 8294.21 dan seluruh parameter tidak signifikan.

ARIMA(0,1,1)

model2.da=Arima(train.diff, order=c(0,1,1),method="ML")
summary(model2.da) 
## Series: train.diff 
## ARIMA(0,1,1) 
## 
## Coefficients:
##           ma1
##       -1.0000
## s.e.   0.0118
## 
## sigma^2 = 1.886e+13:  log likelihood = -4044.37
## AIC=8092.73   AICc=8092.78   BIC=8099.71
## 
## Training set error measures:
##                    ME    RMSE     MAE MPE MAPE     MASE       ACF1
## Training set 135060.8 4324919 3056400 Inf  Inf 0.645519 -0.1862011
lmtest::coeftest(model2.da) 
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.999999   0.011787 -84.838 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Didapatkan nilai AIC sebesar 8092.73 dan seluruh parameter signifikan.

ARIMA(1,1,1)

model3.da=Arima(train.diff, order=c(1,1,1),method="ML")
summary(model3.da)
## Series: train.diff 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.1820  -1.0000
## s.e.   0.0631   0.0129
## 
## sigma^2 = 1.828e+13:  log likelihood = -4040.29
## AIC=8086.58   AICc=8086.68   BIC=8097.04
## 
## Training set error measures:
##                    ME    RMSE     MAE MPE MAPE      MASE        ACF1
## Training set 165358.7 4249451 2985626 Inf  Inf 0.6305715 -0.01640867
lmtest::coeftest(model3.da)
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ar1 -0.181974   0.063139  -2.8821  0.00395 ** 
## ma1 -0.999999   0.012881 -77.6340  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Didapatkan nilai AIC sebesar 8086.58 dan seluruh parameter signifikan.

ARIMA(1,1,2)

model4.da=Arima(train.diff, order=c(1,1,2),method="ML")
summary(model4.da)
## Series: train.diff 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1     ma2
##       0.0101  -1.2017  0.2017
## s.e.  0.2608   0.2528  0.2525
## 
## sigma^2 = 1.831e+13:  log likelihood = -4040.01
## AIC=8088.02   AICc=8088.19   BIC=8101.97
## 
## Training set error measures:
##                    ME    RMSE     MAE MPE MAPE   MASE        ACF1
## Training set 176891.6 4243691 2984340 NaN  Inf 0.6303 -0.00534698
lmtest::coeftest(model4.da)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.010056   0.260820  0.0386    0.9692    
## ma1 -1.201729   0.252797 -4.7537 1.997e-06 ***
## ma2  0.201734   0.252452  0.7991    0.4242    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Didapatkan nilai AIC sebesar 8088.02 dan hanya parameter ma1 signifikan.

ARIMA(2,1,2)

model5.da=Arima(train.diff, order=c(2,1,2),method="ML")
summary(model5.da)
## Series: train.diff 
## ARIMA(2,1,2) 
## 
## Coefficients:
##           ar1      ar2      ma1      ma2
##       -0.9764  -0.2263  -0.2061  -0.7939
## s.e.   0.1143   0.0624   0.1018   0.1014
## 
## sigma^2 = 1.808e+13:  log likelihood = -4037.98
## AIC=8085.97   AICc=8086.22   BIC=8103.41
## 
## Training set error measures:
##                  ME    RMSE     MAE MPE MAPE      MASE        ACF1
## Training set 174139 4207933 2999397 NaN  Inf 0.6334799 0.001234364
lmtest::coeftest(model5.da)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.976416   0.114276 -8.5443 < 2.2e-16 ***
## ar2 -0.226283   0.062436 -3.6242 0.0002898 ***
## ma1 -0.206053   0.101756 -2.0250 0.0428695 *  
## ma2 -0.793942   0.101433 -7.8273 4.985e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Didapatkan nilai AIC sebesar 8085.97 dan seluruh parameter signifikan.

ARIMA(2,1,3)

model6.da=Arima(train.diff, order=c(2,1,3),method="ML")
summary(model6.da)  
## Series: train.diff 
## ARIMA(2,1,3) 
## 
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): NaNs produced
##           ar1      ar2     ma1      ma2      ma3
##       -1.8311  -0.9794  0.8465  -0.8702  -0.9672
## s.e.      NaN      NaN     NaN   0.0225      NaN
## 
## sigma^2 = 1.762e+13:  log likelihood = -4034.81
## AIC=8081.61   AICc=8081.97   BIC=8102.54
## 
## Training set error measures:
##                    ME    RMSE     MAE MPE MAPE      MASE       ACF1
## Training set 137023.4 4145921 2934886 NaN  Inf 0.6198551 -0.1339302
lmtest::coeftest(model6.da) 
## Warning in sqrt(diag(se)): NaNs produced
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -1.831055        NaN     NaN       NaN    
## ar2 -0.979383        NaN     NaN       NaN    
## ma1  0.846472        NaN     NaN       NaN    
## ma2 -0.870232   0.022489 -38.697 < 2.2e-16 ***
## ma3 -0.967179        NaN     NaN       NaN    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Didapatkan nilai AIC sebesar 8081.61 dan hanya parameter ma2 signifikan.

Berdasarkan pendugaan parameter di atas, nilai AIC terkecil 8085.97 dimiliki oleh model ARIMA(2,1,2) dan parameter model ARIMA(2,1,2) juga seluruhnya signifikan sehingga model yang dipilih adalah model ARIMA(2,1,2)

Overfitting

Tahapan selanjutnya adalah overfitting dilakukan dengan menaikkan orde AR(p) dan MA(q) dari model ARIMA(2,1,2) untuk melihat apakah terdapat model lain yang lebih baik dari model saat ini. Kandidat model overfitting adalah ARIMA(3,1,2) dan ARIMA(2,1,3).

ARIMA(3,1,2)

model5.1da=Arima(train.diff, order=c(3,1,2),method="ML")
summary(model5.1da) 
## Series: train.diff 
## ARIMA(3,1,2) 
## 
## Coefficients:
##           ar1      ar2     ar3      ma1      ma2
##       -0.9512  -0.1930  0.0338  -0.2245  -0.7755
## s.e.   0.1265   0.0913  0.0681   0.1104   0.1100
## 
## sigma^2 = 1.814e+13:  log likelihood = -4037.86
## AIC=8087.73   AICc=8088.08   BIC=8108.66
## 
## Training set error measures:
##                    ME    RMSE     MAE MPE MAPE      MASE         ACF1
## Training set 166545.9 4206294 2998926 NaN  Inf 0.6333805 -0.003975975
lmtest::coeftest(model5.1da)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.951159   0.126508 -7.5186 5.539e-14 ***
## ar2 -0.192956   0.091334 -2.1126   0.03463 *  
## ar3  0.033800   0.068135  0.4961   0.61984    
## ma1 -0.224463   0.110354 -2.0340   0.04195 *  
## ma2 -0.775534   0.110048 -7.0473 1.825e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Didapatkan nilai AIC sebesar 8087.73 dan hanya parameter ar3 tidak signifikan.

ARIMA(2,1,3)

model5.2da=Arima(train.diff, order=c(2,1,3),method="ML")
summary(model5.2da)
## Series: train.diff 
## ARIMA(2,1,3) 
## 
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): NaNs produced
##           ar1      ar2     ma1      ma2      ma3
##       -1.8311  -0.9794  0.8465  -0.8702  -0.9672
## s.e.      NaN      NaN     NaN   0.0225      NaN
## 
## sigma^2 = 1.762e+13:  log likelihood = -4034.81
## AIC=8081.61   AICc=8081.97   BIC=8102.54
## 
## Training set error measures:
##                    ME    RMSE     MAE MPE MAPE      MASE       ACF1
## Training set 137023.4 4145921 2934886 NaN  Inf 0.6198551 -0.1339302
lmtest::coeftest(model5.2da) 
## Warning in sqrt(diag(se)): NaNs produced
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -1.831055        NaN     NaN       NaN    
## ar2 -0.979383        NaN     NaN       NaN    
## ma1  0.846472        NaN     NaN       NaN    
## ma2 -0.870232   0.022489 -38.697 < 2.2e-16 ***
## ma3 -0.967179        NaN     NaN       NaN    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Didapatkan nilai AIC sebesar 8081.61 dan hanya parameter ma2 yang signifikan.

Setelah dilakukan overfitting pada model ARIMA(2,1,2) didapatkan hasil bahwa nilai AIC terkecil 8085.97 tetap dimiliki oleh model ARIMA(2,1,2) dengan seluruh parameternya signifikan sehingga model yang dipilih adalah model ARIMA(2,1,2)

Analisis Sisaan

Model terbaik hasil identifikasi kemudian dicek asumsi sisaannya. Sisaan model ARIMA harus memenuhi asumsi normalitas, kebebasan sisaan, dan kehomogenan ragam. Diagnostik model dilakukan secara eksplorasi dan uji formal.

Eksplorasi Sisaan

sisaan.da <- model5.da$residuals 
par(mfrow=c(2,2)) 
qqnorm(sisaan.da) 
qqline(sisaan.da, col = "blue", lwd = 2) 
plot(c(1:length(sisaan.da)),sisaan.da) 
acf(sisaan.da) 
pacf(sisaan.da) 

par(mfrow = c(1,1))

Berdasarkan plot kuantil-kuantil normal, secara eksplorasi ditunjukkan sisaan tidak menyebar normal ditandai dengan titik titik yang cenderung tidak mengikuti garis \(45^{\circ}\) pada kedua ujungnya. Kemudian dapat dilihat juga lebar pita sisaan yang cenderung tidak sama menandakan bahwa sisaan memiliki ragam yang tidak homogen. Plot ACF dan PACF sisaan ARIMA(1,1,2) juga tidak signifikan pada 20 lag awal yang menandakan saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.

Uji Formal

1. Uji kenormalanan sisaan

ks.test(sisaan.da,"pnorm")  #tak tolak H0 > sisaan menyebar normal
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan.da
## D = 0.50206, p-value < 2.2e-16
## alternative hypothesis: two-sided

Uji formal yang digunakan untuk normalitas adalah uji Kolmogorov-Smirnov (KS).

Hipotesis

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

Berdasarkan uji KS tersebut, didapat p-value sebesar 2.2e-16 yang kurang dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal. Hal ini sesuai dengan hasil eksplorasi menggunakan plot kuantil-kuantil normal.

2. Uji kebebasan sisaan (tidak ada autokorelasi)

Box.test(sisaan.da, type = "Ljung")  #tak tolak H0 > sisaan saling bebas
## 
##  Box-Ljung test
## 
## data:  sisaan.da
## X-squared = 0.00037484, df = 1, p-value = 0.9846

Uji formal untuk kebebasan sisaan menggunakan uji Ljung-Box.

Hipotesis

\(H_0\) : Sisaan saling bebas

\(H_1\) : Sisaan tidak tidak saling bebas

Berdasarkan uji Ljung-Box tersebut, didapat p-value sebesar 0.9846 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Hal ini sesuai dengan hasil eksplorasi menggunakan plot ACF dan PACF.

3. Uji kehomogenan sisaan

Box.test((sisaan.da)^2, type = "Ljung")  #tak tolak H0 > sisaan homogen
## 
##  Box-Ljung test
## 
## data:  (sisaan.da)^2
## X-squared = 5.6939, df = 1, p-value = 0.01702

Hipotesis

\(H_0\) : Ragam sisaan homogen

\(H_1\) : Ragam sisaan tidak homogen

Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, didapat p-value sebesar 0.01702 yang kurang dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa ragam sisaan tidak homogen sehingga perlu dilakukan penanganan.

4. Nilai tengah sisaan sama dengan nol

t.test(sisaan.da, mu = 0, conf.level = 0.95)  #tak tolak h0 > nilai tengah sisaan sama dengan 0
## 
##  One Sample t-test
## 
## data:  sisaan.da
## t = 0.64433, df = 242, p-value = 0.52
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -358232  706510
## sample estimates:
## mean of x 
##    174139

Terakhir, dengan uji-t, akan dicek apakah nilai tengah sisaan sama dengan nol.

Hipotesis

\(H_0\) : nilai tengah sisaan sama dengan 0

\(H_1\) : nilai tengah sisaan tidak sama dengan 0

Berdasarkan uji-ttersebut, didapat p-value sebesar 0.52 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol.

Peramalan ARIMA

Peramalan dilakukan menggunakan fungsi forecast() . Contoh peramalan berikut ini dilakukan untuk 61 hari ke depan.

ramalan.da <- forecast::forecast(model2.da, h = 61) 
ramalan.da
##     Point Forecast    Lo 80   Hi 80    Lo 95   Hi 95
## 245      -394804.5 -5971800 5182191 -8924083 8134474
## 246      -394804.5 -5971800 5182191 -8924083 8134474
## 247      -394804.5 -5971800 5182191 -8924083 8134474
## 248      -394804.5 -5971800 5182191 -8924083 8134474
## 249      -394804.5 -5971800 5182191 -8924083 8134474
## 250      -394804.5 -5971800 5182191 -8924083 8134474
## 251      -394804.5 -5971800 5182191 -8924083 8134474
## 252      -394804.5 -5971800 5182191 -8924083 8134474
## 253      -394804.5 -5971800 5182191 -8924083 8134474
## 254      -394804.5 -5971800 5182191 -8924083 8134474
## 255      -394804.5 -5971800 5182191 -8924083 8134474
## 256      -394804.5 -5971800 5182191 -8924083 8134474
## 257      -394804.5 -5971800 5182191 -8924083 8134474
## 258      -394804.5 -5971800 5182191 -8924083 8134474
## 259      -394804.5 -5971800 5182191 -8924083 8134474
## 260      -394804.5 -5971800 5182191 -8924083 8134474
## 261      -394804.5 -5971800 5182191 -8924083 8134474
## 262      -394804.5 -5971800 5182191 -8924083 8134474
## 263      -394804.5 -5971800 5182191 -8924083 8134474
## 264      -394804.5 -5971800 5182191 -8924083 8134474
## 265      -394804.5 -5971800 5182191 -8924083 8134474
## 266      -394804.5 -5971800 5182191 -8924083 8134474
## 267      -394804.5 -5971800 5182191 -8924083 8134474
## 268      -394804.5 -5971800 5182191 -8924083 8134474
## 269      -394804.5 -5971800 5182191 -8924083 8134474
## 270      -394804.5 -5971800 5182191 -8924083 8134474
## 271      -394804.5 -5971800 5182191 -8924083 8134474
## 272      -394804.5 -5971800 5182191 -8924083 8134474
## 273      -394804.5 -5971800 5182191 -8924083 8134474
## 274      -394804.5 -5971800 5182191 -8924083 8134474
## 275      -394804.5 -5971800 5182191 -8924083 8134474
## 276      -394804.5 -5971800 5182191 -8924083 8134474
## 277      -394804.5 -5971800 5182191 -8924083 8134474
## 278      -394804.5 -5971800 5182191 -8924083 8134474
## 279      -394804.5 -5971800 5182191 -8924083 8134474
## 280      -394804.5 -5971800 5182191 -8924083 8134474
## 281      -394804.5 -5971800 5182191 -8924083 8134474
## 282      -394804.5 -5971800 5182191 -8924083 8134474
## 283      -394804.5 -5971800 5182191 -8924083 8134474
## 284      -394804.5 -5971800 5182191 -8924083 8134474
## 285      -394804.5 -5971800 5182191 -8924083 8134474
## 286      -394804.5 -5971800 5182191 -8924083 8134474
## 287      -394804.5 -5971800 5182191 -8924083 8134474
## 288      -394804.5 -5971800 5182191 -8924083 8134474
## 289      -394804.5 -5971800 5182191 -8924083 8134474
## 290      -394804.5 -5971800 5182191 -8924083 8134474
## 291      -394804.5 -5971800 5182191 -8924083 8134474
## 292      -394804.5 -5971800 5182191 -8924083 8134474
## 293      -394804.5 -5971800 5182191 -8924083 8134474
## 294      -394804.5 -5971800 5182191 -8924083 8134474
## 295      -394804.5 -5971800 5182191 -8924083 8134474
## 296      -394804.5 -5971800 5182191 -8924083 8134474
## 297      -394804.5 -5971800 5182191 -8924083 8134474
## 298      -394804.5 -5971800 5182191 -8924083 8134474
## 299      -394804.5 -5971800 5182191 -8924083 8134474
## 300      -394804.5 -5971800 5182191 -8924083 8134474
## 301      -394804.5 -5971800 5182191 -8924083 8134474
## 302      -394804.5 -5971800 5182191 -8924083 8134474
## 303      -394804.5 -5971800 5182191 -8924083 8134474
## 304      -394804.5 -5971800 5182191 -8924083 8134474
## 305      -394804.5 -5971800 5182191 -8924083 8134474

Plot Ramalan

data.ramalan.da <- ramalan.da$mean
plot(ramalan.da)

Dapat dilihat bahwa ramalan ARIMA(2,1,2) cenderung stabil hingga akhir periode. Selanjutnya, dapat dicari nilai akurasi antara hasil ramalan dengan data uji sebagai berikut.

Transformasi Balik

pt_1 <- train.ts.new[244] #nilai akhir data latih
hasil.forc.Diff <- data.ramalan.da
hasil <- diffinv(hasil.forc.Diff, differences = 1) + pt_1
hasil1 <- abs(hasil)^(1/2)

Akurasi

Selanjutnya, dapat dicari nilai akurasi antara hasil ramalan yang sudah ditransformasi balik dengan data uji sebagai berikut.

perbandingan<-matrix(data=c(head(test.ts, n=61),  hasil1[-1]),
                     nrow = 61, ncol = 2)
colnames(perbandingan)<-c("Aktual","Hasil Forecast")
head(perbandingan)
##      Aktual Hasil Forecast
## [1,]   4500       4506.406
## [2,]   4630       4462.386
## [3,]   4810       4417.928
## [4,]   4830       4373.017
## [5,]   4840       4327.641
## [6,]   4810       4281.784
accuracy(ts(hasil1[-1]), head(test.ts, n=61))
##               ME     RMSE    MAE      MPE     MAPE      ACF1 Theil's U
## Test set 1574.49 1783.637 1574.7 38.11491 38.11958 0.9455622  11.73042

Diperoleh nilai MAPE sebesar 38.11958 yang menandakan bahwa hasil peramalan dari model ARIMA(2,1,2) masih kurang baik.

Penanganan Ketidakhomogenan Ragam Sisaan

Uji Efek ARCH

ARCH(1)

bydArchTest1 <- ArchTest(sisaan.da, lags=1, demean=TRUE)
bydArchTest1
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  sisaan.da
## Chi-squared = 5.5368, df = 1, p-value = 0.01862

Didapatkan p-value sebesar 0.01862 kurang dari taraf nyata 5%.

ARCH(3)

bydArchTest3 <- ArchTest(sisaan.da, lags=3, demean=TRUE)
bydArchTest3
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  sisaan.da
## Chi-squared = 14.173, df = 3, p-value = 0.002679

Didapatkan p-value sebesar 0.002679 kurang dari taraf nyata 5%

ARCH(5)

bydArchTest5 <- ArchTest(sisaan.da, lags=5, demean=TRUE)
bydArchTest5
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  sisaan.da
## Chi-squared = 13.823, df = 5, p-value = 0.01678

Didapatkan p-value sebesar 0.01678 kurang dari taraf nyata 5%

ARCH(7)

bydArchTest7 <- ArchTest(sisaan.da, lags=7, demean=TRUE)
train.diff <- train.ts
bydArchTest7
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  sisaan.da
## Chi-squared = 15.167, df = 7, p-value = 0.03391

Didapatkan p-value sebesar 0.03391 kurang dari taraf nyata 5%

ARCH(10)

bydArchTest10 <- ArchTest(sisaan.da, lags=10, demean=TRUE)
bydArchTest10
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  sisaan.da
## Chi-squared = 20.658, df = 10, p-value = 0.02361

Didapatkan p-value sebesar 0.02361 kurang dari taraf nyata 5% ` ### ARCH(12)

bydArchTest12 <- ArchTest(sisaan.da, lags=12, demean=TRUE)
bydArchTest12
## 
##  ARCH LM-test; Null hypothesis: no ARCH effects
## 
## data:  sisaan.da
## Chi-squared = 21.827, df = 12, p-value = 0.03951

Pada lag ke-12, didapatkan p-value sebesar 0.03951 kurang dari taraf nyata 5%. Hal ini menunjukkan bahwa hingga lag ke-12 masih signifikan sehingga akan dicoba dengan GARCH.

Uji Efek GARCH

GARCH(0,1)

garch11 <- rugarch::ugarchspec(mean.model = list(armaOrder = c(2,2), include.mean = FALSE),
                     variance.model = list(model = "sGARCH",garchOrder=c(0,1)))
garch11 <- rugarch::ugarchfit(data = train.diff, spec = garch11)
garch11
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(0,1)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##          Estimate  Std. Error  t value Pr(>|t|)
## ar1      0.782702    0.215122  3.63840 0.000274
## ar2      0.213048    0.214374  0.99382 0.320313
## ma1      0.082917    0.206947  0.40067 0.688663
## ma2     -0.075122    0.078129 -0.96151 0.336296
## omega 4489.601547  929.500777  4.83012 0.000001
## beta1    0.927337    0.013059 71.01045 0.000000
## 
## Robust Standard Errors:
##          Estimate  Std. Error  t value Pr(>|t|)
## ar1      0.782702  1.5300e-01  5.11575 0.000000
## ar2      0.213048  1.5209e-01  1.40081 0.161271
## ma1      0.082917  1.3804e-01  0.60066 0.548064
## ma2     -0.075122  9.6195e-02 -0.78094 0.434840
## omega 4489.601547  1.1489e+03  3.90775 0.000093
## beta1    0.927337  1.3927e-02 66.58419 0.000000
## 
## LogLikelihood : -1842.312 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       15.150
## Bayes        15.236
## Shibata      15.149
## Hannan-Quinn 15.185
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       13.81 0.0002025
## Lag[2*(p+q)+(p+q)-1][11]     15.72 0.0000000
## Lag[4*(p+q)+(p+q)-1][19]     17.88 0.0027819
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      57.19 3.963e-14
## Lag[2*(p+q)+(p+q)-1][2]     57.19 7.772e-16
## Lag[4*(p+q)+(p+q)-1][5]     57.20 5.551e-16
## d.o.f=1
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[2] 0.0004003 0.500 2.000  0.9840
## ARCH Lag[4] 0.0097585 1.397 1.611  0.9994
## ARCH Lag[6] 0.0159181 2.222 1.500  1.0000
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  4.0074
## Individual Statistics:             
## ar1   0.10969
## ar2   0.11280
## ma1   0.04891
## ma2   0.04289
## omega 0.23320
## beta1 0.20630
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.49 1.68 2.12
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                     t-value      prob sig
## Sign Bias            1.1070 2.694e-01    
## Negative Sign Bias   0.2227 8.240e-01    
## Positive Sign Bias  13.7289 5.150e-32 ***
## Joint Effect       190.1269 5.731e-41 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     31.41      0.03638
## 2    30     39.28      0.09644
## 3    40     42.89      0.30816
## 4    50     69.52      0.02844
## 
## 
## Elapsed time : 0.09767413

Didapatkan nilai AIC sebesar 15.150 dengan beberapa parameter yang tidak signifikan (p-value < 0.05) dan semua ARCH LM Test signifikan (p-value > 0.05)

GARCH(0,2)

garch11 <- rugarch::ugarchspec(mean.model = list(armaOrder = c(2,2), include.mean = FALSE),
                     variance.model = list(model = "sGARCH",garchOrder=c(0,2)))
garch11 <- rugarch::ugarchfit(data = train.diff, spec = garch11)
garch11
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(0,2)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##          Estimate  Std. Error  t value Pr(>|t|)
## ar1    8.0202e-01  1.4792e-01  5.42205  0.00000
## ar2    1.9340e-01  1.4734e-01  1.31257  0.18933
## ma1    3.2344e-02  1.3648e-01  0.23699  0.81266
## ma2   -3.9306e-02  6.2043e-02 -0.63353  0.52638
## omega  2.5976e+04  2.0308e+04  1.27906  0.20087
## beta1  6.5937e-01  6.3011e-01  1.04644  0.29536
## beta2  0.0000e+00  2.8385e-01  0.00000  1.00000
## 
## Robust Standard Errors:
##          Estimate  Std. Error  t value Pr(>|t|)
## ar1    8.0202e-01  7.2165e-02 11.11375 0.000000
## ar2    1.9340e-01  7.1043e-02  2.72224 0.006484
## ma1    3.2344e-02  4.4166e-02  0.73233 0.463965
## ma2   -3.9306e-02  5.8288e-02 -0.67434 0.500093
## omega  2.5976e+04  1.7424e+04  1.49083 0.136006
## beta1  6.5937e-01  1.0502e+00  0.62787 0.530091
## beta2  0.0000e+00  8.2034e-01  0.00000 1.000000
## 
## LogLikelihood : -1836.126 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       15.108
## Bayes        15.208
## Shibata      15.106
## Hannan-Quinn 15.148
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       14.65 0.0001297
## Lag[2*(p+q)+(p+q)-1][11]     17.51 0.0000000
## Lag[4*(p+q)+(p+q)-1][19]     20.85 0.0001696
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      57.82 2.875e-14
## Lag[2*(p+q)+(p+q)-1][5]     57.83 3.331e-16
## Lag[4*(p+q)+(p+q)-1][9]     57.86 5.551e-15
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]  0.001172 0.500 2.000  0.9727
## ARCH Lag[5]  0.023178 1.440 1.667  0.9984
## ARCH Lag[7]  0.040755 2.315 1.543  0.9999
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  7.008
## Individual Statistics:             
## ar1   0.06847
## ar2   0.07070
## ma1   0.08334
## ma2   0.08418
## omega 0.20811
## beta1 0.12297
## beta2 0.10427
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.69 1.9 2.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                     t-value      prob sig
## Sign Bias            0.9876 3.243e-01    
## Negative Sign Bias   0.3855 7.002e-01    
## Positive Sign Bias  13.9774 7.553e-33 ***
## Joint Effect       197.0475 1.833e-42 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     31.41      0.03638
## 2    30     42.23      0.05355
## 3    40     47.48      0.16550
## 4    50     73.21      0.01408
## 
## 
## Elapsed time : 0.1932759

Didapatkan nilai AIC sebesar 15.108 dengan beberapa parameter yang tidak signifikan (p-value < 0.05) dan semua ARCH LM Test signifikan (p-value > 0.05)

GARCH(0,3)

garch11 <- rugarch::ugarchspec(mean.model = list(armaOrder = c(2,2), include.mean = FALSE),
                     variance.model = list(model = "sGARCH",garchOrder=c(0,3)))
garch11 <- rugarch::ugarchfit(data = train.diff, spec = garch11)
garch11
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(0,3)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##          Estimate  Std. Error  t value Pr(>|t|)
## ar1    7.8386e-01  2.7260e-01  2.87551 0.004034
## ar2    2.1193e-01  2.7172e-01  0.77995 0.435420
## ma1    5.8515e-02  2.7060e-01  0.21625 0.828796
## ma2   -6.5732e-02  8.4599e-02 -0.77698 0.437168
## omega  1.4566e+04  2.8259e+03  5.15442 0.000000
## beta1  8.0160e-01  3.7989e+00  0.21100 0.832883
## beta2  0.0000e+00  5.9028e+00  0.00000 1.000000
## beta3  0.0000e+00  2.2281e+00  0.00000 1.000000
## 
## Robust Standard Errors:
##          Estimate  Std. Error   t value Pr(>|t|)
## ar1    7.8386e-01  2.1352e-01  3.671210 0.000241
## ar2    2.1193e-01  2.1241e-01  0.997726 0.318412
## ma1    5.8515e-02  1.9560e-01  0.299157 0.764821
## ma2   -6.5732e-02  1.0982e-01 -0.598525 0.549490
## omega  1.4566e+04  3.9702e+04  0.366879 0.713709
## beta1  8.0160e-01  8.4002e+00  0.095426 0.923977
## beta2  0.0000e+00  1.2550e+01  0.000000 1.000000
## beta3  0.0000e+00  4.6899e+00  0.000000 1.000000
## 
## LogLikelihood : -1837.29 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       15.125
## Bayes        15.240
## Shibata      15.123
## Hannan-Quinn 15.172
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       14.75 0.0001229
## Lag[2*(p+q)+(p+q)-1][11]     17.03 0.0000000
## Lag[4*(p+q)+(p+q)-1][19]     19.55 0.0006005
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       57.59 3.231e-14
## Lag[2*(p+q)+(p+q)-1][8]      57.62 3.109e-15
## Lag[4*(p+q)+(p+q)-1][14]     57.65 2.677e-13
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]   0.02087 0.500 2.000  0.8851
## ARCH Lag[6]   0.03285 1.461 1.711  0.9976
## ARCH Lag[8]   0.05196 2.368 1.583  0.9999
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  -0.8131
## Individual Statistics:             
## ar1   0.08624
## ar2   0.09022
## ma1   0.08901
## ma2   0.04663
## omega 0.13935
## beta1 0.11061
## beta2 0.10692
## beta3 0.10563
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.89 2.11 2.59
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                     t-value      prob sig
## Sign Bias            0.9625 3.368e-01    
## Negative Sign Bias   0.3537 7.238e-01    
## Positive Sign Bias  13.9097 1.274e-32 ***
## Joint Effect       195.5345 3.890e-42 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     29.77     0.054793
## 2    30     41.98     0.056350
## 3    40     49.77     0.115762
## 4    50     75.26     0.009324
## 
## 
## Elapsed time : 0.240587

Didapatkan nilai AIC sebesar 15.125 dengan beberapa parameter yang tidak signifikan (p-value < 0.05) dan semua ARCH LM Test signifikan (p-value > 0.05)

GARCH(1,1)

garch11 <- rugarch::ugarchspec(mean.model = list(armaOrder = c(2,2), include.mean = FALSE),
                     variance.model = list(model = "sGARCH",garchOrder=c(1,1)))
garch11 <- rugarch::ugarchfit(data = train.diff, spec = garch11)
garch11
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,1)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##          Estimate  Std. Error  t value Pr(>|t|)
## ar1    8.5422e-02  2.7219e-02   3.1384 0.001699
## ar2    9.0265e-01  2.6904e-02  33.5510 0.000000
## ma1    1.0411e+00  5.6628e-02  18.3848 0.000000
## ma2    1.9940e-01  5.4570e-02   3.6540 0.000258
## omega  4.4896e+03  1.8959e+03   2.3681 0.017881
## alpha1 3.1253e-01  1.0247e-01   3.0500 0.002289
## beta1  6.7078e-01  4.8102e-02  13.9449 0.000000
## 
## Robust Standard Errors:
##          Estimate  Std. Error  t value Pr(>|t|)
## ar1    8.5422e-02  4.4195e-02  1.93284 0.053256
## ar2    9.0265e-01  4.3809e-02 20.60400 0.000000
## ma1    1.0411e+00  2.6112e-01  3.98703 0.000067
## ma2    1.9940e-01  2.0546e-01  0.97048 0.331808
## omega  4.4896e+03  2.8451e+03  1.57801 0.114564
## alpha1 3.1253e-01  2.3762e-01  1.31524 0.188428
## beta1  6.7078e-01  1.1516e-01  5.82475 0.000000
## 
## LogLikelihood : -1812.97 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       14.918
## Bayes        15.018
## Shibata      14.916
## Hannan-Quinn 14.958
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic p-value
## Lag[1]                       2.686  0.1013
## Lag[2*(p+q)+(p+q)-1][11]     5.707  0.6773
## Lag[4*(p+q)+(p+q)-1][19]     8.301  0.7508
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic   p-value
## Lag[1]                      24.77 6.465e-07
## Lag[2*(p+q)+(p+q)-1][5]     25.17 7.164e-07
## Lag[4*(p+q)+(p+q)-1][9]     25.54 7.035e-06
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]  0.001074 0.500 2.000  0.9739
## ARCH Lag[5]  0.212553 1.440 1.667  0.9625
## ARCH Lag[7]  0.515738 2.315 1.543  0.9769
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  11.0502
## Individual Statistics:              
## ar1    0.02157
## ar2    0.02330
## ma1    0.08116
## ma2    0.17956
## omega  0.55381
## alpha1 0.07554
## beta1  0.21605
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.69 1.9 2.35
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias           0.8967 3.708e-01    
## Negative Sign Bias  0.2514 8.017e-01    
## Positive Sign Bias  5.5385 8.024e-08 ***
## Joint Effect       30.8762 9.026e-07 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     23.87       0.2012
## 2    30     34.36       0.2263
## 3    40     41.25       0.3727
## 4    50     54.36       0.2777
## 
## 
## Elapsed time : 0.09920788
res1<-garch11@fit$residuals
res1.ts<-ts(res1)
plot.ts(res1.ts)

Didapatkan nilai AIC sebesar 14.918 dengan semua parameter signifikan (p-value < 0.05) dan semua ARCH LM Test signifikan (p-value > 0.05)

GARCH(1,2)

garch11 <- rugarch::ugarchspec(mean.model = list(armaOrder = c(2,2), include.mean = FALSE),
                     variance.model = list(model = "sGARCH",garchOrder=c(1,2)))
garch11 <- rugarch::ugarchfit(data = train.diff, spec = garch11)
garch11
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,2)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##          Estimate  Std. Error  t value Pr(>|t|)
## ar1    1.1236e-01  3.5387e-02   3.1753 0.001497
## ar2    8.7815e-01  3.4727e-02  25.2869 0.000000
## ma1    1.0581e+00  7.5801e-02  13.9589 0.000000
## ma2    2.1352e-01  6.5974e-02   3.2365 0.001210
## omega  2.2926e+04  7.2390e+03   3.1670 0.001540
## alpha1 2.7909e-01  1.0074e-01   2.7705 0.005597
## beta1  4.5706e-01  3.5635e-01   1.2826 0.199635
## beta2  0.0000e+00  2.5258e-01   0.0000 1.000000
## 
## Robust Standard Errors:
##          Estimate  Std. Error  t value Pr(>|t|)
## ar1    1.1236e-01  7.7560e-02  1.44875 0.147408
## ar2    8.7815e-01  7.5759e-02 11.59136 0.000000
## ma1    1.0581e+00  4.7489e-01  2.22813 0.025872
## ma2    2.1352e-01  3.8175e-01  0.55934 0.575932
## omega  2.2926e+04  1.8829e+04  1.21759 0.223380
## alpha1 2.7909e-01  3.7352e-01  0.74718 0.454954
## beta1  4.5706e-01  8.2399e-01  0.55469 0.579109
## beta2  0.0000e+00  4.6487e-01  0.00000 1.000000
## 
## LogLikelihood : -1797.321 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       14.798
## Bayes        14.912
## Shibata      14.796
## Hannan-Quinn 14.844
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       7.404 6.509e-03
## Lag[2*(p+q)+(p+q)-1][11]     9.747 6.343e-08
## Lag[4*(p+q)+(p+q)-1][19]    12.774 1.254e-01
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       34.57 4.120e-09
## Lag[2*(p+q)+(p+q)-1][8]      35.49 6.723e-09
## Lag[4*(p+q)+(p+q)-1][14]     35.73 2.702e-07
## d.o.f=3
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[4]   0.04104 0.500 2.000  0.8395
## ARCH Lag[6]   0.04710 1.461 1.711  0.9960
## ARCH Lag[8]   0.29617 2.368 1.583  0.9943
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  3.24
## Individual Statistics:              
## ar1    0.04115
## ar2    0.02869
## ma1    0.05864
## ma2    0.15408
## omega  0.18709
## alpha1 0.10686
## beta1  0.06329
## beta2  0.05850
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          1.89 2.11 2.59
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value      prob sig
## Sign Bias           0.8740 3.830e-01    
## Negative Sign Bias  0.2145 8.303e-01    
## Positive Sign Bias  6.5600 3.301e-10 ***
## Joint Effect       43.3776 2.046e-09 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     17.15       0.5799
## 2    30     43.95       0.0371
## 3    40     37.64       0.5319
## 4    50     50.26       0.4232
## 
## 
## Elapsed time : 0.381856

Didapatkan nilai AIC sebesar 14.798 dengan beberapa parameter yang tidak signifikan (p-value < 0.05) dan semua ARCH LM Test signifikan (p-value > 0.05)

GARCH(1,3)

garch11 <- rugarch::ugarchspec(mean.model = list(armaOrder = c(2,2), include.mean = FALSE),
                     variance.model = list(model = "sGARCH",garchOrder=c(1,3)))
garch11 <- rugarch::ugarchfit(data = train.diff, spec = garch11)
garch11
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(1,3)
## Mean Model   : ARFIMA(2,0,2)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##           Estimate  Std. Error  t value Pr(>|t|)
## ar1     8.0997e-01  2.0537e-02 39.43931 0.000000
## ar2     1.8320e-01  2.0411e-02  8.97545 0.000000
## ma1     1.5122e-02  6.2524e-02  0.24187 0.808883
## ma2    -8.9870e-03  4.3041e-02 -0.20880 0.834608
## omega   2.5829e+04  5.3677e+03  4.81185 0.000001
## alpha1  1.9445e-01  7.4672e-02  2.60410 0.009212
## beta1   4.8010e-01  2.8530e-01  1.68277 0.092420
## beta2   0.0000e+00  4.2823e-01  0.00000 1.000000
## beta3   0.0000e+00  1.8359e-01  0.00000 1.000000
## 
## Robust Standard Errors:
##           Estimate  Std. Error   t value Pr(>|t|)
## ar1     8.0997e-01  3.4270e-03 236.35483 0.000000
## ar2     1.8320e-01  2.3240e-03  78.82647 0.000000
## ma1     1.5122e-02  4.6900e-02   0.32244 0.747119
## ma2    -8.9870e-03  2.9808e-02  -0.30149 0.763043
## omega   2.5829e+04  1.1254e+04   2.29504 0.021731
## alpha1  1.9445e-01  9.0796e-02   2.14164 0.032222
## beta1   4.8010e-01  4.6106e-01   1.04130 0.297738
## beta2   0.0000e+00  6.1991e-01   0.00000 1.000000
## beta3   0.0000e+00  3.0254e-01   0.00000 1.000000
## 
## LogLikelihood : -1829.232 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       15.067
## Bayes        15.196
## Shibata      15.065
## Hannan-Quinn 15.119
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       14.97 1.092e-04
## Lag[2*(p+q)+(p+q)-1][11]     17.82 0.000e+00
## Lag[4*(p+q)+(p+q)-1][19]     21.54 8.437e-05
## d.o.f=4
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       56.20 6.539e-14
## Lag[2*(p+q)+(p+q)-1][11]     56.35 6.917e-14
## Lag[4*(p+q)+(p+q)-1][19]     56.49 1.671e-11
## d.o.f=4
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[5]   0.01370 0.500 2.000  0.9068
## ARCH Lag[7]   0.04489 1.473 1.746  0.9966
## ARCH Lag[9]   0.18047 2.402 1.619  0.9984
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  1.4328
## Individual Statistics:              
## ar1    0.10597
## ar2    0.10597
## ma1    0.07169
## ma2    0.08839
## omega  0.14660
## alpha1 0.07529
## beta1  0.07726
## beta2  0.06010
## beta3  0.04698
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.1 2.32 2.82
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                     t-value      prob sig
## Sign Bias            1.2576 2.098e-01    
## Negative Sign Bias   0.1474 8.829e-01    
## Positive Sign Bias  13.7288 5.153e-32 ***
## Joint Effect       189.8988 6.419e-41 ***
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     27.15      0.10126
## 2    30     39.52      0.09202
## 3    40     45.84      0.20970
## 4    50     59.69      0.14093
## 
## 
## Elapsed time : 0.248148

Didapatkan nilai AIC sebesar 15.067 dengan beberapa parameter yang tidak signifikan (p-value < 0.05) dan semua ARCH LM Test signifikan (p-value > 0.05)

Setelah dilakukan uji efek GARCH, didapatkan nilai AIC terkecil 14.918 yang dimiliki oleh model GARCH(1,1) dengan nilai p-value yang lebih besar dari taraf nyata 5% pada ARCH LM Test dan seluruh parameter model GARCH(1,1) signifikan sehingga model yang dipilih adalah model GARCH(1,1).

Peramalan ARIMA-GARCH

Peramalan dilakukan untuk 61 periode ke depan.

Plot Peramalan Pada Model

forc<- ugarchforecast(fitORspec = garch11, n.ahead = 61, n.roll = 0)
plot(forc, which= 1)

## Akurasi

Selanjutnya, dapat dicari nilai akurasi antara hasil ramalan yang sudah ditransformasi balik dengan data uji sebagai berikut.

pt_1 <- train.ts.new[244] #nilai akhlir data latih
hasil.forc.Diff <- forc@forecast$seriesFor[,1]
hasilgarch <- diffinv(hasil.forc.Diff, differences = 1) + pt_1
hasilakhir <- abs(hasilgarch)^(1/2) 
perbandingan <- data.frame("Aktual"= test.ts,
                           "Ramalan" = hasilakhir[-1])
head(perbandingan,10)
##    Aktual  Ramalan
## 1    4500 4550.498
## 2    4630 4550.993
## 3    4810 4551.484
## 4    4830 4551.973
## 5    4840 4552.460
## 6    4810 4552.943
## 7    5300 4553.423
## 8    4980 4553.901
## 9    4540 4554.376
## 10   4540 4554.848

Nilai MAPE

Periode  <-c(test.ts[245:305]) 
dataframe <- data.frame(Periode, perbandingan) 
T <- nrow(dataframe) 
MAPE <- 1/T*sum(abs((dataframe$Aktual-dataframe$Ramalan)/dataframe$Aktual)*100)
MAPE
## [1] 9.316642

Diperoleh nilai MAPE sebesar 9.316642 yang menandakan bahwa hasil peramalan dari model ARIMA(2,1,2)-GARCH(1,1) sudah sangat baik.

Pemodelan

ARIMA(2,1,2)-GARCH(1,1)

coef(garch11)
##           ar1           ar2           ma1           ma2         omega 
##  8.099724e-01  1.832018e-01  1.512248e-02 -8.986756e-03  2.582864e+04 
##        alpha1         beta1         beta2         beta3 
##  1.944534e-01  4.800979e-01  1.478834e-09  6.002961e-15

Model Mean ARIMA (2,1,2)

\[ Y_t = -0.8099Y_{t-1}-0.1832Y_{t-2}+ e_t -0.0151e_{t-1} -0.0089e_{t-2} \]

Model Varian GARCH (1,1)

\[ σ^2 = 2.5828×10^4+ 0.19444e^2_{t-1} + 0.4801σ^2_{t-1} + 1.478×10^-9σ^2_{t-2} + 6.003×10^-15σ^2_{t-3} \]

Kesimpulan

Peramalan data harga saham PT Unilever Indonesia Tbk untuk 61 periode (pekan) ke depan dengan menggunakan model ARIMA(2,1,2) tidak dapat mengatasi adanya efek heteroskedastisitas dengan nilai MAPE sebesar 38,11%. Nilai MAPE tersebut menunjukkan bahwa peramalan masih belum cukup baik. Efek heteroskedastisitas tersebut ditangani dengan model ARCH-GARCH.Lag model ARCH masih signifikan hingga lag 12 sehingga selanjutnya ditangani dengan model GARCH dan dibandingkan nilai AIC dan signifikansi parameter setiap model. Didapatkan model GARCH(1,1) yang memiliki nilai AIC terkecil dan semua parameter signifikan. Model ARIMA(2,1,2)-GARCH(1,1) sebagai model terbaik untuk meramalkan harga saham PT Unilever Indonesia Tbk dengan nilai MAPE sebesar 9,31%.