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
datadeath<-rio::import("https://raw.githubusercontent.com/mrnabilnaufal07/mpdw/main/Data/Data-Nabil.csv")

#Menghapus data periode
datadeath <- datadeath$Meninggal

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

Eksplorasi Data

Plot Data Penuh

plot.ts(datadeath.ts, xlab="waktu", ylab="Jumlah Kematian", main="Plot Kasus Meninggal Karena Diabetes")

Berdasarkan plot data deret waktu, terlihat bahwa data memiliki pola siklik. Berdasarkan pola data, pembagian data latih dan data uji ditetapkan dengan proporsi 70%:30%.

Plot Data Latih

deathtrain<-datadeath[1:120]
train.ts<-ts(deathtrain)
plot.ts(train.ts, lty=1, xlab="waktu", ylab="death", main="Plot Death Train")

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

deathtest<-datadeath[121:173]
test.ts<-ts(deathtest)
plot.ts(test.ts, lty=1, xlab="waktu", ylab="death", main="Plot death Test")

Uji Stasioneritas Data

Plot ACF

acf(train.ts)

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(train.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.ts
## Dickey-Fuller = -2.6413, Lag order = 4, p-value = 0.3102
## 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.3102 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

Plot Box-Cox

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

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda
## [1] -0.15
#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]
##   [1] -0.86 -0.85 -0.84 -0.83 -0.82 -0.81 -0.80 -0.79 -0.78 -0.77 -0.76 -0.75
##  [13] -0.74 -0.73 -0.72 -0.71 -0.70 -0.69 -0.68 -0.67 -0.66 -0.65 -0.64 -0.63
##  [25] -0.62 -0.61 -0.60 -0.59 -0.58 -0.57 -0.56 -0.55 -0.54 -0.53 -0.52 -0.51
##  [37] -0.50 -0.49 -0.48 -0.47 -0.46 -0.45 -0.44 -0.43 -0.42 -0.41 -0.40 -0.39
##  [49] -0.38 -0.37 -0.36 -0.35 -0.34 -0.33 -0.32 -0.31 -0.30 -0.29 -0.28 -0.27
##  [61] -0.26 -0.25 -0.24 -0.23 -0.22 -0.21 -0.20 -0.19 -0.18 -0.17 -0.16 -0.15
##  [73] -0.14 -0.13 -0.12 -0.11 -0.10 -0.09 -0.08 -0.07 -0.06 -0.05 -0.04 -0.03
##  [85] -0.02 -0.01  0.00  0.01  0.02  0.03  0.04  0.05  0.06  0.07  0.08  0.09
##  [97]  0.10  0.11  0.12  0.13  0.14  0.15  0.16  0.17  0.18  0.19  0.20  0.21
## [109]  0.22  0.23  0.24  0.25  0.26  0.27  0.28  0.29  0.30  0.31  0.32  0.33
## [121]  0.34  0.35  0.36  0.37  0.38  0.39  0.40  0.41  0.42  0.43  0.44  0.45
## [133]  0.46  0.47  0.48  0.49  0.50  0.51  0.52  0.53  0.54  0.55  0.56

Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar -0.15 dan pada selang kepercayaan 95% nilai memiliki batas bawah -0.86 dan batas atas 0.56. Selang tersebut tidak memuat nilai satu sehingga dapat dikatakan bahwa data yang digunakan tidak stasioner dalam ragam.

Penanganan Ketidakstasioneran Data

train.diff<-diff(train.ts,differences = 1) 
plot.ts(train.diff, lty=1, xlab="waktu", ylab="Data Difference 1 Death", main="Plot Difference Death")

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, 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

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 = -6.9241, Lag order = 4, 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 dalam rataan 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 2, sehingga jika plot ACF dianggap tails of, maka model tentatifnya adalah ARIMA(2,1,0).

Jika baik plot ACF maupun plot PACF keduanya dianggap tails of, maka model yang terbentuk adalah ARIMA(1,1,1)

Plot EACF

eacf(train.diff)
## 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  x 
## 1 x x o o o o o o o o o  o  o  x 
## 2 x o x o o o o o o o o  o  o  o 
## 3 x o x o o o o o o o o  o  o  o 
## 4 x x x o o o o o o o o  o  o  o 
## 5 x o x x o o o o o o o  o  o  o 
## 6 x o x o o o o o o o o  o  o  o 
## 7 o x x 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(0,1,2), ARIMA(1,1,2), dan ARIMA(2,1,1)

Pendugaan Parameter Model Tentatif

ARIMA(0,1,1)

model1.da=Arima(train.diff, order=c(0,1,1),method="ML")
summary(model1.da) #AIC=861.26 
## Series: train.diff 
## ARIMA(0,1,1) 
## 
## Coefficients:
##          ma1
##       -1.000
## s.e.   0.021
## 
## sigma^2 = 81.05:  log likelihood = -428.63
## AIC=861.26   AICc=861.37   BIC=866.8
## 
## Training set error measures:
##                       ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set -0.02097883 8.926839 7.156795 NaN  Inf 0.5530464 -0.6061335
lmtest::coeftest(model1.da) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -1.000000   0.021007 -47.602 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(2,1,0)

model2.da=Arima(train.diff, order=c(2,1,0),method="ML")
summary(model2.da) #AIC=864.86 
## Series: train.diff 
## ARIMA(2,1,0) 
## 
## Coefficients:
##           ar1      ar2
##       -1.1229  -0.4893
## s.e.   0.0799   0.0797
## 
## sigma^2 = 85.28:  log likelihood = -429.43
## AIC=864.86   AICc=865.07   BIC=873.18
## 
## Training set error measures:
##                       ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set 0.007483314 9.117691 6.984824 NaN  Inf 0.5397572 -0.1757061
lmtest::coeftest(model2.da) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -1.122912   0.079927 -14.0492 < 2.2e-16 ***
## ar2 -0.489348   0.079750  -6.1361  8.46e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(1,1,1)

model3.da=Arima(train.diff, order=c(1,1,1),method="ML")
summary(model3.da) #AIC=810.27
## Series: train.diff 
## ARIMA(1,1,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.5995  -1.0000
## s.e.   0.0726   0.0219
## 
## sigma^2 = 51.57:  log likelihood = -402.13
## AIC=810.27   AICc=810.48   BIC=818.58
## 
## Training set error measures:
##                      ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set 0.03990198 7.089887 5.486836 NaN  Inf 0.4239991 -0.1514496
lmtest::coeftest(model3.da) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.599516   0.072619  -8.2556 < 2.2e-16 ***
## ma1 -0.999999   0.021878 -45.7086 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(0,1,2)

model4.da=Arima(train.diff, order=c(0,1,2),method="ML")
summary(model4.da) #AIC=799.27
## Series: train.diff 
## ARIMA(0,1,2) 
## 
## Coefficients:
##           ma1     ma2
##       -1.7523  0.7523
## s.e.   0.0661  0.0620
## 
## sigma^2 = 46.11:  log likelihood = -396.64
## AIC=799.27   AICc=799.48   BIC=807.59
## 
## Training set error measures:
##                      ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set -0.2056647 6.704611 5.374915 NaN  Inf 0.4153504 -0.1336757
lmtest::coeftest(model4.da) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -1.752257   0.066142 -26.492 < 2.2e-16 ***
## ma2  0.752266   0.061965  12.140 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(1,1,2)

model5.da=Arima(train.diff, order=c(1,1,2),method="ML")
summary(model5.da) #AIC=798.02 
## Series: train.diff 
## ARIMA(1,1,2) 
## 
## Coefficients:
##           ar1      ma1     ma2
##       -0.2213  -1.6456  0.6456
## s.e.   0.1226   0.1100  0.1077
## 
## sigma^2 = 45.36:  log likelihood = -395.01
## AIC=798.02   AICc=798.37   BIC=809.1
## 
## Training set error measures:
##                       ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set -0.09217524 6.620467 5.269615 NaN  Inf 0.4072132 0.00484327
lmtest::coeftest(model5.da) #terdapat parameter tidak signifikan
## 
## z test of coefficients:
## 
##     Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.22134    0.12265  -1.8047   0.07113 .  
## ma1 -1.64562    0.11001 -14.9585 < 2.2e-16 ***
## ma2  0.64562    0.10770   5.9948 2.037e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(2,1,1)

model6.da=Arima(train.diff, order=c(2,1,1),method="ML")
summary(model6.da) #AIC=805.48
## Series: train.diff 
## ARIMA(2,1,1) 
## 
## Coefficients:
##           ar1      ar2      ma1
##       -0.7445  -0.2354  -1.0000
## s.e.   0.0891   0.0888   0.0226
## 
## sigma^2 = 48.88:  log likelihood = -398.74
## AIC=805.48   AICc=805.83   BIC=816.56
## 
## Training set error measures:
##                      ME     RMSE     MAE MPE MAPE      MASE       ACF1
## Training set 0.07668715 6.873008 5.30659 NaN  Inf 0.4100705 -0.0506046
lmtest::coeftest(model6.da) #seluruh parameter tidak signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.744495   0.089120  -8.3538 < 2.2e-16 ***
## ar2 -0.235372   0.088753  -2.6520  0.008002 ** 
## ma1 -0.999999   0.022557 -44.3319 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan pendugaan parameter di atas, nilai AIC terkecil dimiliki oleh model ARIMA(1,1,2) namun ada parameter yang tidak signifikan. Namun, model ARIMA(0,1,2) memiliki nilai AIC kedua terkecil dan parameter model ARIMA(0,1,2) juga seluruhnya signifikan sehingga model yang dipilih adalah model ARIMA(0,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

#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 menyebar normal ditandai dengan titik titik yang cenderung mengikuti garis \(45^{\circ}\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung sama menandakan bahwa sisaan memiliki ragam yang homogen. Plot ACF dan PACF sisaan ARIMA(0,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) Sisaan Menyebar Normal 
tseries::jarque.bera.test(sisaan.da)  #tak tolak H0 > sisaan menyebar normal
## 
##  Jarque Bera Test
## 
## data:  sisaan.da
## X-squared = 0.69234, df = 2, p-value = 0.7074

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 KS tersebut, didapat p-value sebesar 0.7074 yang lebih dari taraf nyata 5% sehingga tidak tolak \(H_0\) dan menandakan bahwa sisaan 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 = 2.1805, df = 1, p-value = 0.1398

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.1398 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas. Hal ini 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 = 0.89689, df = 1, p-value = 0.3436

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.3436 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.33337, df = 118, p-value = 0.7394
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.427332  1.016003
## sample estimates:
##  mean of x 
## -0.2056647

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.7394 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(0,1,2) untuk melihat apakah terdapat model lain yang lebih baik dari model saat ini. Kandidat model overfitting adalah ARIMA(1,1,2) dan ARIMA(0,1,3).

Model ARIMA(1,1,2)

model.overfit1=Arima(train.diff, order=c(1,1,2),method="ML")
summary(model.overfit1) #798.02
## Series: train.diff 
## ARIMA(1,1,2) 
## 
## Coefficients:
##           ar1      ma1     ma2
##       -0.2213  -1.6456  0.6456
## s.e.   0.1226   0.1100  0.1077
## 
## sigma^2 = 45.36:  log likelihood = -395.01
## AIC=798.02   AICc=798.37   BIC=809.1
## 
## Training set error measures:
##                       ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set -0.09217524 6.620467 5.269615 NaN  Inf 0.4072132 0.00484327
lmtest::coeftest(model.overfit1) #ar1 tidak signifikan
## 
## z test of coefficients:
## 
##     Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.22134    0.12265  -1.8047   0.07113 .  
## ma1 -1.64562    0.11001 -14.9585 < 2.2e-16 ***
## ma2  0.64562    0.10770   5.9948 2.037e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model ARIMA(0,1,3)

model.overfit2=Arima(train.diff, order=c(0,1,3),method="ML")
summary(model.overfit2) #798.74
## Series: train.diff 
## ARIMA(0,1,3) 
## 
## Coefficients:
##           ma1     ma2      ma3
##       -1.8454  0.9753  -0.1299
## s.e.   0.0922  0.1549   0.0808
## 
## sigma^2 = 45.62:  log likelihood = -395.37
## AIC=798.74   AICc=799.1   BIC=809.83
## 
## Training set error measures:
##                      ME     RMSE      MAE MPE MAPE      MASE        ACF1
## Training set -0.1070217 6.639746 5.310216 NaN  Inf 0.4103507 -0.02325722
lmtest::coeftest(model.overfit2) #ma3 tidak signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ma1 -1.845431   0.092214 -20.0125 < 2.2e-16 ***
## ma2  0.975311   0.154948   6.2944 3.085e-10 ***
## ma3 -0.129878   0.080831  -1.6068    0.1081    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#model yang dipilih adalah model awal, yaitu ARIMA(0,1,2)

Peramalan

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

#---FORECAST---#
ramalan.da <- forecast::forecast(model4.da, h = 30) 
ramalan.da
##     Point Forecast      Lo 80     Hi 80     Lo 95    Hi 95
## 121    0.450222594  -8.290773  9.191218 -12.91798 13.81842
## 122   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 123   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 124   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 125   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 126   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 127   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 128   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 129   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 130   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 131   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 132   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 133   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 134   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 135   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 136   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 137   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 138   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 139   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 140   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 141   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 142   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 143   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 144   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 145   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 146   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 147   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 148   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 149   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
## 150   -0.005521503 -10.897576 10.886533 -16.66348 16.65244
data.ramalan.da <- ramalan.da$mean
plot(ramalan.da)

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

pt_1 <- train.ts[120] #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)

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,]     28       29.45022
##  [2,]     38       29.44470
##  [3,]     42       29.43918
##  [4,]     30       29.43366
##  [5,]     31       29.42814
##  [6,]     28       29.42262
##  [7,]     37       29.41709
##  [8,]     39       29.41157
##  [9,]     40       29.40605
## [10,]     34       29.40053
## [11,]     46       29.39501
## [12,]     43       29.38949
## [13,]     30       29.38396
## [14,]     26       29.37844
## [15,]     46       29.37292
## [16,]     46       29.36740
## [17,]     44       29.36188
## [18,]     33       29.35636
## [19,]     39       29.35084
## [20,]     48       29.34531
## [21,]     37       29.33979
## [22,]     54       29.33427
## [23,]     30       29.32875
## [24,]     32       29.32323
## [25,]     39       29.31771
## [26,]     46       29.31219
## [27,]     37       29.30666
## [28,]     35       29.30114
## [29,]     26       29.29562
## [30,]     32       29.29010
accuracy(ts(hasil[-1]), head(test.ts, n=30))
##                ME     RMSE      MAE      MPE     MAPE       ACF1 Theil's U
## Test set 7.829839 10.60193 8.466299 18.07318 20.46849 0.09501679  1.118426

Didapati bahwa nilai MAPE mencapai 20.47%, yang mengindikasikan bahwa model prediksi rata-rata mengalami deviasi sekitar 21% dari nilai sebenarnya. Saat dilakukan pengecekan stasioneritas, ditemukan bahwa data tidak menunjukkan stasioneritas baik dalam rataan maupun dalam ragam. Namun, dalam analisis ini, hanya masalah ketidakstasioneran dalam rataan yang telah diatasi. Sayangnya, ketidakstasioneran dalam ragam masih menjadi permasalahan yang belum terpecahkan. Hal ini dapat berpotensi menjadi salah satu penyebab mengapa nilai MAPE dari model masih relatif tinggi.