Packages yang digunakakan

library(ggplot2)
library(tsibble)
## 
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, union
library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(MASS)
library(forecast)
library(TSA)
## 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(TTR)
library(aTSA)
## 
## Attaching package: 'aTSA'
## 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(graphics)

Penyiapan Data

#Input data
meta = rio::import("https://raw.githubusercontent.com/mrnabilnaufal07/mpdw/main/KELOMPOK%206/saham_meta%20-%20Close.csv")
n = NROW(meta)

#Menghapus data Date dan Periode
meta.close <- meta$Close
View(meta)

#Mengubah tipe data menjadi time series
meta.ts<-ts(meta.close)

Eksplorasi Data

Plot Data Penuh

plot.ts(meta.ts, xlab="waktu", ylab="Harga Close Saham Meta (USD)")

ggsave("Plot harga closing.png", dpi=1000)
## Saving 7 x 5 in image

Berdasarkan plot data deret waktu, terlihat bahwa data memiliki pola tren naik pada periode 1 hingga 190 an. Namun, memasuki periode 200 tren data mulai menurun Berdasarkan pola data, pembagian data latih dan data uji ditetapkan dengan proporsi 90%:10%.

Plot Data Latih

meta.train<-meta$Close[1:275]
train.ts<-ts(meta.train)
plot.ts(train.ts, lty=1, xlab="waktu", ylab="Harga Close Saham Meta (USD)", main="Plot Harga Close Saham Meta (USD)")

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

Plot Data Uji

meta.test<-meta$Close[276:305]
test.ts<-ts(meta.test)
plot.ts(test.ts, lty=1, xlab="waktu", ylab="Harga Close Saham Meta (USD)", main="Plot Harga Close Saham Meta (USD)")

Kestasioneran Dalam Ragam

Plot Box-Cox

index <- seq(1:275)
bc = boxcox(train.ts~index, lambda = seq(-0.7,0.4,by=0.01))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda
## [1] -0.17
#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]
##  [1] -5.200000e-01 -5.100000e-01 -5.000000e-01 -4.900000e-01 -4.800000e-01
##  [6] -4.700000e-01 -4.600000e-01 -4.500000e-01 -4.400000e-01 -4.300000e-01
## [11] -4.200000e-01 -4.100000e-01 -4.000000e-01 -3.900000e-01 -3.800000e-01
## [16] -3.700000e-01 -3.600000e-01 -3.500000e-01 -3.400000e-01 -3.300000e-01
## [21] -3.200000e-01 -3.100000e-01 -3.000000e-01 -2.900000e-01 -2.800000e-01
## [26] -2.700000e-01 -2.600000e-01 -2.500000e-01 -2.400000e-01 -2.300000e-01
## [31] -2.200000e-01 -2.100000e-01 -2.000000e-01 -1.900000e-01 -1.800000e-01
## [36] -1.700000e-01 -1.600000e-01 -1.500000e-01 -1.400000e-01 -1.300000e-01
## [41] -1.200000e-01 -1.100000e-01 -1.000000e-01 -9.000000e-02 -8.000000e-02
## [46] -7.000000e-02 -6.000000e-02 -5.000000e-02 -4.000000e-02 -3.000000e-02
## [51] -2.000000e-02 -1.000000e-02  1.110223e-16  1.000000e-02  2.000000e-02
## [56]  3.000000e-02  4.000000e-02  5.000000e-02  6.000000e-02  7.000000e-02
## [61]  8.000000e-02  9.000000e-02  1.000000e-01  1.100000e-01  1.200000e-01
## [66]  1.300000e-01  1.400000e-01  1.500000e-01  1.600000e-01  1.700000e-01
## [71]  1.800000e-01  1.900000e-01  2.000000e-01
ggsave("Boxcox sebelum trans.png", dpi=1000)
## Saving 7 x 5 in image

Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar -0.17 dan pada selang kepercayaan 95% nilai memiliki batas bawah -0.52 dan batas atas 0.20. Selang tersebut tidak memuat nilai satu sehingga dapat dikatakan bahwa data yang digunakan tidak stasioner dalam ragam. Nilai lambda = 0 dipilih untuk melakukan transformasi,

Penanganan Ketidakstasioneran dalam ragam (Lambda = 0)

# transformasi dengan nilai lambda = 0
data_boxcox = log(train.ts)

par(mfrow=c(1,2))
plot.ts(train.ts, lty=1, xlab="waktu", ylab="Harga Awal", main="Plot Data Train")
plot.ts(data_boxcox, lty=1, xlab="waktu", ylab="Harga Transformasi", main="Plot Data Train (Transformasi)")

Plot Box-Cox

index2 <- seq(1:275)
bc2 = boxcox(data_boxcox~index, lambda = seq(-2,2.5,by=0.01))

#Nilai Rounded Lambda
lambda2 <- bc2$x[which.max(bc2$y)]
lambda2
## [1] 0.08
#SK
bc2$x[bc2$y > max(bc2$y) - 1/2 * qchisq(.95,1)]
##   [1] -1.81 -1.80 -1.79 -1.78 -1.77 -1.76 -1.75 -1.74 -1.73 -1.72 -1.71 -1.70
##  [13] -1.69 -1.68 -1.67 -1.66 -1.65 -1.64 -1.63 -1.62 -1.61 -1.60 -1.59 -1.58
##  [25] -1.57 -1.56 -1.55 -1.54 -1.53 -1.52 -1.51 -1.50 -1.49 -1.48 -1.47 -1.46
##  [37] -1.45 -1.44 -1.43 -1.42 -1.41 -1.40 -1.39 -1.38 -1.37 -1.36 -1.35 -1.34
##  [49] -1.33 -1.32 -1.31 -1.30 -1.29 -1.28 -1.27 -1.26 -1.25 -1.24 -1.23 -1.22
##  [61] -1.21 -1.20 -1.19 -1.18 -1.17 -1.16 -1.15 -1.14 -1.13 -1.12 -1.11 -1.10
##  [73] -1.09 -1.08 -1.07 -1.06 -1.05 -1.04 -1.03 -1.02 -1.01 -1.00 -0.99 -0.98
##  [85] -0.97 -0.96 -0.95 -0.94 -0.93 -0.92 -0.91 -0.90 -0.89 -0.88 -0.87 -0.86
##  [97] -0.85 -0.84 -0.83 -0.82 -0.81 -0.80 -0.79 -0.78 -0.77 -0.76 -0.75 -0.74
## [109] -0.73 -0.72 -0.71 -0.70 -0.69 -0.68 -0.67 -0.66 -0.65 -0.64 -0.63 -0.62
## [121] -0.61 -0.60 -0.59 -0.58 -0.57 -0.56 -0.55 -0.54 -0.53 -0.52 -0.51 -0.50
## [133] -0.49 -0.48 -0.47 -0.46 -0.45 -0.44 -0.43 -0.42 -0.41 -0.40 -0.39 -0.38
## [145] -0.37 -0.36 -0.35 -0.34 -0.33 -0.32 -0.31 -0.30 -0.29 -0.28 -0.27 -0.26
## [157] -0.25 -0.24 -0.23 -0.22 -0.21 -0.20 -0.19 -0.18 -0.17 -0.16 -0.15 -0.14
## [169] -0.13 -0.12 -0.11 -0.10 -0.09 -0.08 -0.07 -0.06 -0.05 -0.04 -0.03 -0.02
## [181] -0.01  0.00  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09  0.10
## [193]  0.11  0.12  0.13  0.14  0.15  0.16  0.17  0.18  0.19  0.20  0.21  0.22
## [205]  0.23  0.24  0.25  0.26  0.27  0.28  0.29  0.30  0.31  0.32  0.33  0.34
## [217]  0.35  0.36  0.37  0.38  0.39  0.40  0.41  0.42  0.43  0.44  0.45  0.46
## [229]  0.47  0.48  0.49  0.50  0.51  0.52  0.53  0.54  0.55  0.56  0.57  0.58
## [241]  0.59  0.60  0.61  0.62  0.63  0.64  0.65  0.66  0.67  0.68  0.69  0.70
## [253]  0.71  0.72  0.73  0.74  0.75  0.76  0.77  0.78  0.79  0.80  0.81  0.82
## [265]  0.83  0.84  0.85  0.86  0.87  0.88  0.89  0.90  0.91  0.92  0.93  0.94
## [277]  0.95  0.96  0.97  0.98  0.99  1.00  1.01  1.02  1.03  1.04  1.05  1.06
## [289]  1.07  1.08  1.09  1.10  1.11  1.12  1.13  1.14  1.15  1.16  1.17  1.18
## [301]  1.19  1.20  1.21  1.22  1.23  1.24  1.25  1.26  1.27  1.28  1.29  1.30
## [313]  1.31  1.32  1.33  1.34  1.35  1.36  1.37  1.38  1.39  1.40  1.41  1.42
## [325]  1.43  1.44  1.45  1.46  1.47  1.48  1.49  1.50  1.51  1.52  1.53  1.54
## [337]  1.55  1.56  1.57  1.58  1.59  1.60  1.61  1.62  1.63  1.64  1.65  1.66
## [349]  1.67  1.68  1.69  1.70  1.71  1.72  1.73  1.74  1.75  1.76  1.77  1.78
## [361]  1.79  1.80  1.81  1.82  1.83  1.84  1.85  1.86  1.87  1.88  1.89  1.90
## [373]  1.91  1.92  1.93  1.94  1.95  1.96  1.97  1.98  1.99  2.00  2.01

Plot Boxcox memuat nilai satu sehingga dapat dikatakan bahwa data yang digunakan sudah stasioner dalam ragam.

Kestasioneran Dalam Rataan

Plot ACF

acf(data_boxcox)

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

Uji ADF

tseries::adf.test(data_boxcox)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_boxcox
## Dickey-Fuller = -1.7717, Lag order = 6, p-value = 0.6722
## 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.6722 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa data tidak stasioner dalam rataan. Hal ini sesuai dengan hasil eksplorasi menggunakan plot time series dan plot ACF, sehingga ketidakstasioneran model kedepannya harus ditangani.

Penanganan Ketidakstasioneran dalam rataan

train.diff<-diff(data_boxcox,differences = 1) 
plot.ts(train.diff, lty=1, xlab="waktu", ylab="Data Difference 1 Saham Meta (USD)", main="Plot Harga Close Saham Meta (USD)")

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)

Plot ACF

acf(train.diff)

Berdasarkan plot tersebut, tidak terlihat adanya cut off pada lag 1. Artinya data sudah stasioner dalam rataan. Namun, akan dipastikan kembali dengan pengujian formal.

Uji ADF

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.6838, 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. Dalam hal ini ketidakstasioneran data dalam rataan sudah berhasil ditangani dan dapat dilanjutkan ke identifikasi model.

Identifikasi Model

Plot ACF

acf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot ACF tidak cuts off.

Plot PACF

pacf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot PACF tidak cuts off.

Tidak dapat mengidentifikasi model menggunakan plot ACF dan PACF

Plot EACF

eacf(data_boxcox)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x x x x x x x x x  x  x  x 
## 1 o o o o o o o o o o o  o  o  o 
## 2 x o o x o o o o o o o  o  o  o 
## 3 x x o x o o o o o o o  o  o  o 
## 4 x x x x o o o o o o o  o  o  o 
## 5 x x o o o o o o o o o  o  o  o 
## 6 x x o o o o o o o o o  o  o  o 
## 7 x o o o o o o 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(1,1,0), ARIMA(1,1,1) ARIMA(1,1,3), ARIMA(2,1,2).

Pendugaan Parameter Model Tentatif

ARIMA(1,1,0)

model1.da=Arima(data_boxcox, order=c(1,1,0),method="ML")
summary(model1.da) #AIC=-792.34    
## Series: data_boxcox 
## ARIMA(1,1,0) 
## 
## Coefficients:
##           ar1
##       -0.0580
## s.e.   0.0603
## 
## sigma^2 = 0.003213:  log likelihood = 398.17
## AIC=-792.34   AICc=-792.29   BIC=-785.11
## 
## Training set error measures:
##                        ME       RMSE        MAE        MPE     MAPE      MASE
## Training set 0.0005747142 0.05647723 0.04057241 0.00454333 0.772686 0.9963321
##                      ACF1
## Training set 0.0007095721
lmtest::coeftest(model1.da) #seluruh parameter tidak signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 -0.058032   0.060273 -0.9628   0.3356

ARIMA(1,1,1)

model2.da=Arima(data_boxcox, order=c(1,1,1),method="ML")
summary(model2.da) #AIC=-790.36
## Series: data_boxcox 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1     ma1
##       -0.1447  0.0876
## s.e.   0.5950  0.5961
## 
## sigma^2 = 0.003225:  log likelihood = 398.18
## AIC=-790.36   AICc=-790.27   BIC=-779.52
## 
## Training set error measures:
##                        ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.0005725358 0.05647433 0.04057693 0.00453167 0.7727692 0.9964432
##                       ACF1
## Training set -0.0006388511
lmtest::coeftest(model2.da) #seluruh parameter tidak signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 -0.144711   0.595013 -0.2432   0.8078
## ma1  0.087634   0.596117  0.1470   0.8831

ARIMA(1,1,3)

model3.da=Arima(data_boxcox, order=c(1,1,3),method="ML")
summary(model3.da) #AIC=-790.25
## Series: data_boxcox 
## ARIMA(1,1,3) 
## 
## Coefficients:
##          ar1      ma1     ma2     ma3
##       0.2962  -0.3664  0.0633  0.0970
## s.e.  0.3159   0.3107  0.0649  0.0637
## 
## sigma^2 = 0.003202:  log likelihood = 400.13
## AIC=-790.25   AICc=-790.03   BIC=-772.19
## 
## Training set error measures:
##                        ME       RMSE        MAE         MPE      MAPE      MASE
## Training set 0.0005472865 0.05606997 0.04070347 0.005139251 0.7752716 0.9995507
##                     ACF1
## Training set 0.002435993
lmtest::coeftest(model3.da) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1  0.296222   0.315918  0.9377   0.3484
## ma1 -0.366376   0.310693 -1.1792   0.2383
## ma2  0.063296   0.064893  0.9754   0.3294
## ma3  0.097021   0.063692  1.5233   0.1277

ARIMA(2,1,2)

model4.da=Arima(data_boxcox, order=c(2,1,2),method="ML")
summary(model4.da) #AIC=-792.21
## Series: data_boxcox 
## ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       0.8845  -0.6766  -0.9619  0.7931
## s.e.  0.1782   0.1851   0.1490  0.1545
## 
## sigma^2 = 0.003178:  log likelihood = 401.11
## AIC=-792.21   AICc=-791.99   BIC=-774.14
## 
## Training set error measures:
##                        ME       RMSE        MAE         MPE     MAPE      MASE
## Training set 0.0005357446 0.05586216 0.04043367 0.004527275 0.770401 0.9929253
##                     ACF1
## Training set 0.008159536
lmtest::coeftest(model4.da) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.88450    0.17825  4.9622 6.969e-07 ***
## ar2 -0.67665    0.18508 -3.6559 0.0002563 ***
## ma1 -0.96193    0.14900 -6.4561 1.074e-10 ***
## ma2  0.79313    0.15445  5.1351 2.819e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
knitr::kable(
  cbind(c("ARIMA (1,1,0)","ARIMA (1,1,1)","ARIMA (1,1,3)","ARIMA (2,1,2)"),
        c(model1.da$aic,model2.da$aic,model3.da$aic,model4.da$aic)),
  col.names=c("Model","AIC")
)
Model AIC
ARIMA (1,1,0) -792.335565123414
ARIMA (1,1,1) -790.363803598181
ARIMA (1,1,3) -790.252191907801
ARIMA (2,1,2) -792.210272086609

Berdasarkan pendugaan parameter di atas, model ARIMA(2,1,2) dengan seluruh parameternya signifikan pada taraf 5% dan AIC terkecil. Oleh karena itu model ARIMA(2,1,2) dipilih untuk proses berikutnya

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

#Eksplorasi 
sisaan.da <- model4.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}\). 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(0,1,1) juga ada yang signifikan pada 20 lag awal yang menandakan tidak saling bebas. Kondisi ini akan diuji lebih lanjut dengan uji formal.

Uji Formal

#1) Sisaan Menyebar Normal 
tseries::jarque.bera.test(sisaan.da);  #tak tolak H0 > sisaan menyebar normal
## 
##  Jarque Bera Test
## 
## data:  sisaan.da
## X-squared = 165.7, df = 2, p-value < 2.2e-16
nortest::ad.test(sisaan.da)
## 
##  Anderson-Darling normality test
## 
## data:  sisaan.da
## A = 2.3549, p-value = 5.651e-06

Selain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Jarque Bera . Hipotesis pada uji Jarque Bera adalah sebagai berikut.

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

Berdasarkan uji Jarque-Bera tersebut, didapat p-value kurang dari 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) Sisaan saling bebas/tidak ada autokorelasi 
Box.test(sisaan.da, type = "Ljung")  #tak tolak H0 > sisaan saling bebas
## 
##  Box-Ljung test
## 
## data:  sisaan.da
## X-squared = 0.018509, df = 1, p-value = 0.8918

Selanjutnya akan dilakukan uji formal untuk kebebasan sisaan menggunakan uji Ljung-Box. Hipotesis yang digunakan adalah sebagai berikut.

\(H_0\) : Sisaan saling bebas

\(H_1\) : Sisaan tidak tidak saling bebas

Berdasarkan uji Ljung-Box tersebut, didapat p-value sebesar 0.8918 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas/tidak ada autokorelasi. Hal ini tidak sesuai dengan hasil eksplorasi.

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

Hipotesis yang digunakan untuk uji kehomogenan ragam adalah sebagai berikut.

\(H_0\) : Ragam sisaan homogen

\(H_1\) : Ragam sisaan tidak homogen

Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, didapat p-value sebesar 0.2504 yang lebih dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa ragam sisaan homogen.

#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.15876, df = 274, p-value = 0.874
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.006107697  0.007179186
## sample estimates:
##    mean of x 
## 0.0005357446

Terakhir, dengan uji-t, akan dicek apakah nilai tengah sisaan sama dengan nol. Hipotesis yang diujikan sebagai berikut.

\(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.874 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol. Hal ini sama dengan hasil eksplorasi.

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(2,1,3) dan ARIMA(3,1,2).

Model ARIMA(2,1,3)

model.overfit1=Arima(data_boxcox, order=c(2,1,3),method="ML")
summary(model.overfit1) #-790.28  
## Series: data_boxcox 
## ARIMA(2,1,3) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     ma3
##       0.8273  -0.6387  -0.8960  0.7367  0.0253
## s.e.  0.3046   0.2659   0.3111  0.2932  0.0954
## 
## sigma^2 = 0.003189:  log likelihood = 401.14
## AIC=-790.28   AICc=-789.97   BIC=-768.6
## 
## Training set error measures:
##                        ME      RMSE        MAE         MPE      MAPE      MASE
## Training set 0.0005298075 0.0558551 0.04042824 0.004519325 0.7703156 0.9927918
##                       ACF1
## Training set -0.0003483687
lmtest::coeftest(model.overfit1) #ma3 tidak signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)   
## ar1  0.827344   0.304625  2.7159 0.006609 **
## ar2 -0.638701   0.265866 -2.4023 0.016291 * 
## ma1 -0.896035   0.311092 -2.8803 0.003973 **
## ma2  0.736692   0.293205  2.5125 0.011986 * 
## ma3  0.025332   0.095443  0.2654 0.790692   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model ARIMA(3,1,2)

model.overfit2=Arima(data_boxcox, order=c(3,1,2),method="ML")
summary(model.overfit2) #-790.28 
## Series: data_boxcox 
## ARIMA(3,1,2) 
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2
##       0.8611  -0.6674  0.0202  -0.9303  0.7690
## s.e.  0.2191   0.2103  0.0789   0.2112  0.2053
## 
## sigma^2 = 0.003189:  log likelihood = 401.14
## AIC=-790.28   AICc=-789.96   BIC=-768.6
## 
## Training set error measures:
##                        ME       RMSE        MAE         MPE      MAPE      MASE
## Training set 0.0005300233 0.05585548 0.04042884 0.004520478 0.7703243 0.9928066
##                      ACF1
## Training set 0.0001175852
lmtest::coeftest(model.overfit2) #ar3 tidak signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.861141   0.219081  3.9307 8.470e-05 ***
## ar2 -0.667431   0.210323 -3.1734 0.0015068 ** 
## ar3  0.020182   0.078934  0.2557 0.7981980    
## ma1 -0.930307   0.211155 -4.4058 1.054e-05 ***
## ma2  0.769010   0.205300  3.7458 0.0001798 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model yang dipilih tetap model ARIMA(2,1,2)

Peramalan

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

#---FORECAST---#
ramalan.da <- forecast::forecast(model4.da, h = 30) 
ramalan.da
##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 276       5.386578 5.314328 5.458828 5.276081 5.497075
## 277       5.388176 5.289875 5.486478 5.237838 5.538515
## 278       5.382257 5.261508 5.503007 5.197587 5.566928
## 279       5.375940 5.232742 5.519139 5.156937 5.594943
## 280       5.374358 5.209990 5.538726 5.122979 5.625737
## 281       5.377233 5.194719 5.559748 5.098101 5.656365
## 282       5.380847 5.183282 5.578411 5.078698 5.682995
## 283       5.382097 5.171386 5.592809 5.059842 5.704353
## 284       5.380759 5.157536 5.603981 5.039369 5.722148
## 285       5.378728 5.143035 5.614422 5.018266 5.739191
## 286       5.377838 5.129851 5.625825 4.998574 5.757102
## 287       5.378425 5.118738 5.638112 4.981267 5.775582
## 288       5.379546 5.108953 5.650139 4.965710 5.793382
## 289       5.380141 5.099302 5.660980 4.950634 5.809647
## 290       5.379908 5.089205 5.670611 4.935316 5.824500
## 291       5.379300 5.078927 5.679673 4.919918 5.838681
## 292       5.378919 5.069051 5.688788 4.905016 5.852823
## 293       5.378994 5.059888 5.698101 4.890963 5.867026
## 294       5.379318 5.051293 5.707343 4.877647 5.880989
## 295       5.379554 5.042914 5.716194 4.864708 5.894400
## 296       5.379543 5.034523 5.724563 4.851880 5.907206
## 297       5.379374 5.026146 5.732602 4.839159 5.919589
## 298       5.379232 5.017948 5.740516 4.826695 5.931769
## 299       5.379221 5.010042 5.748400 4.814610 5.943831
## 300       5.379307 5.002409 5.756204 4.802892 5.955721
## 301       5.379390 4.994949 5.763832 4.791437 5.967344
## 302       5.379406 4.987574 5.771239 4.780150 5.978662
## 303       5.379364 4.980273 5.778454 4.769007 5.989720
## 304       5.379315 4.973087 5.785544 4.758042 6.000588
## 305       5.379301 4.966052 5.792550 4.747292 6.011311
data.ramalan.da <- ramalan.da$mean
plot(ramalan.da, xlab="Periode (Tahun)", col = "blue", lwd=1.5)

Berdasarkan hasil plot ramalan di atas, 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.

pt_1 <- (exp(data_boxcox[length(data_boxcox)])) #nilai akhir data latih
hasil.forc.Diff <- data.ramalan.da
hasil <- diffinv(hasil.forc.Diff, differences = 1) + pt_1

#has.1 sama hasilnta dengan: cumsum(c(pt_1,hasil.forc.Diff))
ts.plot(train.ts,hasil, ylab="Jumlah kematian",
        xlab = "Periode (Minggu)",
        col = "blue", lwd=1.5)

perbandingan.da<-matrix(data=c(head(test.ts, n=30), hasil[-1]),
                     nrow = 30, ncol = 2)
colnames(perbandingan.da)<-c("Aktual","Hasil Forecast")
perbandingan.da
##       Aktual Hasil Forecast
##  [1,] 221.49       221.4866
##  [2,] 212.89       226.8748
##  [3,] 240.32       232.2570
##  [4,] 232.78       237.6330
##  [5,] 233.81       243.0073
##  [6,] 245.64       248.3845
##  [7,] 262.04       253.7654
##  [8,] 272.61       259.1475
##  [9,] 264.95       264.5282
## [10,] 281.00       269.9070
## [11,] 288.73       275.2848
## [12,] 286.98       280.6632
## [13,] 290.53       286.0428
## [14,] 308.87       291.4229
## [15,] 294.26       296.8028
## [16,] 325.48       302.1821
## [17,] 310.73       307.5611
## [18,] 301.64       312.9400
## [19,] 283.25       318.3194
## [20,] 285.50       323.6989
## [21,] 296.38       329.0785
## [22,] 297.89       334.4578
## [23,] 300.31       339.8371
## [24,] 299.08       345.2163
## [25,] 300.21       350.5956
## [26,] 315.43       355.9750
## [27,] 314.69       361.3544
## [28,] 308.65       366.7338
## [29,] 296.73       372.1131
## [30,] 302.66       377.4924
accuracy(ts(hasil[-1]), head(test.ts, n=30))
##                 ME     RMSE     MAE      MPE     MAPE      ACF1 Theil's U
## Test set -16.97444 32.54952 24.2732 -5.71993 8.239161 0.8240846  2.358312

Diperoleh nilai MAPE sebesar 8.239161%. Nilai ini cukup baik dan menunjukkan kalau model layak untuk digunakan.

Kesimpulan: Model ARIMA(2,1,2) terpilih.