Pendugaan Parameter, Diagnostik Model, dan Peramalan data deret waktu

Packages dan Data

library(ggplot2)
library(tsibble)
library(tseries)
library(MASS)
library(forecast)
library(TSA)
library(TTR)
library(aTSA)
library(graphics)
library(rio)

Pada permasalahan kali ini akan digunakan data pribadi yang telah diupload pada github yang tertera. Data ini merupakan data nilai tukar rupiah terhadap mata uang asing (USD) per bulan pada abad ke-21 (terhitung sejak januari 2001).

df<- import("https://raw.githubusercontent.com/fax17/MPDW/main/data/data.csv")
data <- df[-c(1,2)] #mengambil nilai tukarnya saja
data.ts <- ts(df$USD)

Eksplorasi Data

Plot Data Penuh

plot.ts(data.ts, lty=1, xlab="Waktu", ylab="Nilai Tukar", main="Plot Data Nilai Tukar")

Berdasarkan plot data deret waktu, terlihat bahwa data cenderung memiliki trend yang naik meski sempat beberapa kali turun terlebih dahulu. Berdasarkan pola data ini, peneliti akan membagi data latih dan data uji sebanyak 240 dan 31 amatan. Pembagian ini dilandaskan bahwa data uji merupakan data 3 tahun (2 tahun 7 bulan) terakhir setelah tahun 2020 dan data latih adalah 20 tahun pertama di abad ke-21 ini. Pembagian ini subjektif berdasarkan keinginan peneliti.

Plot Data Latih

datatrain<-data[1:240,]
train.ts<-ts(datatrain)
plot.ts(train.ts, lty=1, xlab="Waktu", ylab="Data", main="Plot Data Train")

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

Plot Data Uji

datatest<-data[241:271,]
test.ts<-ts(datatest)
plot.ts(test.ts, lty=1, xlab="Waktu", ylab="Data", main="Plot Data 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.3793, Lag order = 6, p-value = 0.4163
## 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.4163 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:240)
bc = MASS::boxcox(train.ts~index, lambda = seq(-2,2,by=1))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda
## [1] -0.8282828
#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]
##  [1] -1.6363636 -1.5959596 -1.5555556 -1.5151515 -1.4747475 -1.4343434
##  [7] -1.3939394 -1.3535354 -1.3131313 -1.2727273 -1.2323232 -1.1919192
## [13] -1.1515152 -1.1111111 -1.0707071 -1.0303030 -0.9898990 -0.9494949
## [19] -0.9090909 -0.8686869 -0.8282828 -0.7878788 -0.7474747 -0.7070707
## [25] -0.6666667 -0.6262626 -0.5858586 -0.5454545 -0.5050505 -0.4646465
## [31] -0.4242424 -0.3838384 -0.3434343 -0.3030303 -0.2626263 -0.2222222
## [37] -0.1818182 -0.1414141 -0.1010101

Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar -0.8282 dan pada selang kepercayaan 95% nilai memiliki batas bawah -1.6363 dan batas atas -0.1010. Selang tersebut tidak memuat nilai satu sehingga dapat dikatakan bahwa data 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 Data", main="Plot Difference Data")

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 setelah lag ke 2. 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.541, 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 setelah lag ke 2, sehingga jika plot PACF dianggap tails of, maka model tentatifnya adalah ARIMA(0,1,2).

Plot PACF

pacf(train.diff)

Berdasarkan plot tersebut, terlihat bahwa plot PACF cenderung cuts off setelah 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

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

Pendugaan Parameter Model Tentatif

ARIMA(0,1,2)

model1.da <- forecast::Arima(train.diff, order=c(0,1,2),method="ML")
summary(model1.da) #AIC=3517.59
## Series: train.diff 
## ARIMA(0,1,2) 
## 
## Coefficients:
##           ma1      ma2
##       -0.9765  -0.0234
## s.e.   0.0765   0.0753
## 
## sigma^2 = 147615:  log likelihood = -1755.79
## AIC=3517.59   AICc=3517.69   BIC=3528
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE      MASE          ACF1
## Training set 2.113051 381.7885 236.1622 Inf  Inf 0.7080166 -0.0008223877
lmtest::coeftest(model1.da) #Terdapat parameter tidak signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ma1 -0.976514   0.076470 -12.7698   <2e-16 ***
## ma2 -0.023422   0.075327  -0.3109   0.7558    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(2,1,0)

model2.da<- forecast::Arima(train.diff, order=c(2,1,0),method="ML")
summary(model2.da) #AIC=3591.58
## Series: train.diff 
## ARIMA(2,1,0) 
## 
## Coefficients:
##           ar1      ar2
##       -0.5841  -0.3677
## s.e.   0.0604   0.0604
## 
## sigma^2 = 205658:  log likelihood = -1792.79
## AIC=3591.58   AICc=3591.69   BIC=3602
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set -5.432038 450.6401 264.1241 Inf  Inf 0.7918467 -0.1108034
lmtest::coeftest(model2.da) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.584145   0.060400 -9.6712 < 2.2e-16 ***
## ar2 -0.367718   0.060443 -6.0837 1.174e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(1,1,1)

model3.da<- forecast::Arima(train.diff, order=c(1,1,1),method="ML")
summary(model3.da) #AIC= 3517.61
## Series: train.diff 
## ARIMA(1,1,1) 
## 
## Coefficients:
##          ar1      ma1
##       0.0169  -1.0000
## s.e.  0.0649   0.0132
## 
## sigma^2 = 147614:  log likelihood = -1755.8
## AIC=3517.61   AICc=3517.71   BIC=3528.03
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
## Training set 2.153531 381.7868 236.2388 Inf  Inf 0.7082461 0.004893674
lmtest::coeftest(model3.da) #Terdapat parameter tidak signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ar1  0.016907   0.064946   0.2603   0.7946    
## ma1 -0.999999   0.013226 -75.6068   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(1,1,0)

model4.da<- forecast::Arima(train.diff, order=c(1,1,0),method="ML")
summary(model4.da) #AIC=3623.91 
## Series: train.diff 
## ARIMA(1,1,0) 
## 
## Coefficients:
##           ar1
##       -0.4264
## s.e.   0.0586
## 
## sigma^2 = 236855:  log likelihood = -1809.96
## AIC=3623.91   AICc=3623.96   BIC=3630.86
## 
## Training set error measures:
##                     ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set -3.466796 484.6369 290.1364 Inf  Inf 0.8698318 -0.1570939
lmtest::coeftest(model4.da) #seluruh parameter signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.426366   0.058594 -7.2766 3.422e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(2,1,2)

model5.da<- forecast::Arima(train.diff, order=c(2,1,2),method="ML")
summary(model5.da) #AIC=3515.81
## Series: train.diff 
## ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2
##       0.6127  -0.1363  -1.6122  0.6123
## s.e.  0.3447   0.0742   0.3491  0.3488
## 
## sigma^2 = 144886:  log likelihood = -1752.91
## AIC=3515.81   AICc=3516.07   BIC=3533.17
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
## Training set 12.86134 376.6363 238.1308 Inf  Inf 0.7139185 0.005676002
lmtest::coeftest(model5.da) #terdapat parameter tidak signifikan namun masih dapat dianggap signifikan pada tarafnyata 10%
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.612710   0.344686  1.7776   0.07547 .  
## ar2 -0.136293   0.074237 -1.8359   0.06637 .  
## ma1 -1.612224   0.349078 -4.6185 3.865e-06 ***
## ma2  0.612266   0.348846  1.7551   0.07924 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(3,1,3)

model6.da<- forecast::Arima(train.diff, order=c(3,1,3),method="ML")
summary(model6.da) #AIC=3517.67
## Series: train.diff 
## ARIMA(3,1,3) 
## 
## Coefficients:
##           ar1      ar2      ar3     ma1      ma2      ma3
##       -1.0868  -0.4561  -0.1762  0.1119  -0.7598  -0.3521
## s.e.   0.3099   0.2828   0.0723  0.3142   0.1808   0.2982
## 
## sigma^2 = 145006:  log likelihood = -1751.83
## AIC=3517.67   AICc=3518.15   BIC=3541.97
## 
## Training set error measures:
##                    ME     RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set 6.352709 375.1792 236.1572 -Inf  Inf 0.7080014 0.001226521
lmtest::coeftest(model6.da) #terdapat parameter tidak signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -1.086837   0.309939 -3.5066 0.0004538 ***
## ar2 -0.456087   0.282798 -1.6128 0.1067956    
## ar3 -0.176192   0.072292 -2.4372 0.0148005 *  
## ma1  0.111914   0.314174  0.3562 0.7216792    
## ma2 -0.759842   0.180831 -4.2019 2.646e-05 ***
## ma3 -0.352065   0.298226 -1.1805 0.2377894    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan pendugaan parameter di atas, nilai AIC terkecil dimana sleuruh parameternya masih bisa dikatakan signifikan dimiliki oleh model ARIMA(2,1,2) sehingga model yang dipilih adalah model ARIMA(2,1,2) atau model tentatif ke 5.

Overfirtting

Overfitting dilakukan dengan menambahkan nilai MA yang semula dari ARIMA(2,1,2) menjadi ARIMA(2,1,3). Dan juga mencoba menambahkan nilai AR yang semula dari ARIMA(2,1,2) menjadi ARIMA(3,1,2)

model7.da<- forecast::Arima(train.diff, order=c(2,1,3),method="ML")
summary(model7.da) #AIC= 3519.89
## Series: train.diff 
## ARIMA(2,1,3) 
## 
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): NaNs produced
##          ar1      ar2      ma1      ma2      ma3
##       -0.308  -0.4124  -0.6549  -0.0312  -0.3139
## s.e.     NaN      NaN      NaN   0.2174      NaN
## 
## sigma^2 = 147120:  log likelihood = -1753.95
## AIC=3519.89   AICc=3520.26   BIC=3540.73
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
## Training set 2.885187 378.7167 234.3079 Inf  Inf 0.7024575 -0.01517315
lmtest::coeftest(model7.da) #terdapat parameter tidak signifikan (model eror)
## Warning in sqrt(diag(se)): NaNs produced
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 -0.307997        NaN     NaN      NaN
## ar2 -0.412397        NaN     NaN      NaN
## ma1 -0.654852        NaN     NaN      NaN
## ma2 -0.031166   0.217422 -0.1433    0.886
## ma3 -0.313893        NaN     NaN      NaN
model8.da<- forecast::Arima(train.diff, order=c(3,1,2),method="ML")
summary(model8.da) #AIC= 3516.19
## Series: train.diff 
## ARIMA(3,1,2) 
## 
## Coefficients:
##           ar1      ar2      ar3      ma1      ma2
##       -0.7580  -0.1272  -0.1865  -0.2285  -0.7715
## s.e.   0.1698   0.0809   0.0680   0.1632   0.1628
## 
## sigma^2 = 144651:  log likelihood = -1752.1
## AIC=3516.19   AICc=3516.56   BIC=3537.03
## 
## Training set error measures:
##                    ME     RMSE      MAE  MPE MAPE    MASE        ACF1
## Training set 8.232221 375.5252 237.3174 -Inf  Inf 0.71148 0.006492241
lmtest::coeftest(model8.da) #Terdapat parameter tidak signifikan
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.757960   0.169849 -4.4626 8.099e-06 ***
## ar2 -0.127188   0.080947 -1.5713  0.116123    
## ar3 -0.186457   0.067985 -2.7426  0.006095 ** 
## ma1 -0.228484   0.163169 -1.4003  0.161427    
## ma2 -0.771507   0.162769 -4.7399 2.138e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diperoleh nilai aic overfitting lebih besar sehingga tetap dipakai model sebelum overfitting.

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 <- 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}\). Kemudian dapat dilihat juga lebar pita sisaan yang cenderung tidak sama menandakan bahwa sisaan memiliki ragam yang heterogen. Plot ACF dan PACF sisaan ARIMA(2,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 
ks.test(sisaan.da,"pnorm")  #tak tolak H0 > sisaan menyebar normal
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan.da
## D = 0.53975, p-value < 2.2e-16
## alternative hypothesis: two-sided

Selain dengan eksplorasi, asumsi tersebut dapat diuji menggunakan uji formal. Pada tahapan ini uji formal yang digunakan untuk normalitas adalah uji Kolmogorov-Smirnov (KS). Hipotesis pada uji KS adalah sebagai berikut.

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

Berdasarkan uji KS tersebut, didapat p-value sebesar 0.00 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.0077969, df = 1, p-value = 0.9296

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

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.0003517 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa ragam sisaan tidak 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.52712, df = 238, p-value = 0.5986
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -35.20516  60.92784
## sample estimates:
## mean of x 
##  12.86134

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

Peramalan

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

Plot Peramalan pada model

#---FORECAST---#
ramalan.da <- forecast::forecast(model5.da, h = 30) 
ramalan.da
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 241      105.82091 -382.9883 594.6302 -641.7483 853.3901
## 242       77.84336 -410.9667 566.6534 -669.7271 825.4138
## 243       43.14379 -449.8932 536.1808 -710.8912 797.1788
## 244       25.69616 -468.8778 520.2701 -730.6895 782.0818
## 245       19.73514 -475.0434 514.5137 -736.9633 776.4336
## 246       18.46076 -476.3225 513.2441 -738.2450 775.1665
## 247       18.49238 -476.2898 513.2746 -738.2117 775.1965
## 248       18.68545 -476.0979 513.4688 -738.0204 775.3913
## 249       18.79943 -475.9848 513.5836 -737.9077 775.5065
## 250       18.84296 -475.9416 513.6275 -737.8647 775.5506
## 251       18.85409 -475.9306 513.6387 -737.8537 775.5619
## 252       18.85498 -475.9297 513.6396 -737.8529 775.5628
## 253       18.85401 -475.9307 513.6387 -737.8538 775.5618
## 254       18.85329 -475.9314 513.6379 -737.8545 775.5611
## 255       18.85298 -475.9317 513.6376 -737.8548 775.5608
## 256       18.85289 -475.9318 513.6375 -737.8549 775.5607
## 257       18.85288 -475.9318 513.6375 -737.8550 775.5607
## 258       18.85288 -475.9318 513.6375 -737.8550 775.5607
## 259       18.85289 -475.9318 513.6375 -737.8550 775.5607
## 260       18.85289 -475.9318 513.6376 -737.8550 775.5607
## 261       18.85289 -475.9318 513.6376 -737.8550 775.5607
## 262       18.85289 -475.9318 513.6376 -737.8550 775.5607
## 263       18.85289 -475.9318 513.6376 -737.8550 775.5607
## 264       18.85289 -475.9318 513.6376 -737.8550 775.5607
## 265       18.85289 -475.9318 513.6376 -737.8550 775.5607
## 266       18.85289 -475.9318 513.6376 -737.8550 775.5607
## 267       18.85289 -475.9318 513.6376 -737.8550 775.5607
## 268       18.85289 -475.9318 513.6376 -737.8550 775.5608
## 269       18.85289 -475.9318 513.6376 -737.8550 775.5608
## 270       18.85289 -475.9318 513.6376 -737.8550 775.5608
data.ramalan.da <- ramalan.da$mean
plot(ramalan.da)

Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa ramalan ARIMA(2,1,2) cenderung stabil hingga akhir periode.

Plot Peramalan pada data awal

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

Dari hasil ini terlihat peramalan yang stabil pula yang mengikuti trend pada data dan cukup menggambarkan ramalan dari datanya.

Akurasi

Selanjutnya, dapat dicari nilai akurasi antara hasil ramalan dengan data uji sebagai berikut.

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,]  14084       14210.82
##  [2,]  14229       14288.66
##  [3,]  14572       14331.81
##  [4,]  14468       14357.50
##  [5,]  14310       14377.24
##  [6,]  14496       14395.70
##  [7,]  14491       14414.19
##  [8,]  14374       14432.88
##  [9,]  14307       14451.68
## [10,]  14199       14470.52
## [11,]  14340       14489.37
## [12,]  14269       14508.23
## [13,]  14381       14527.08
## [14,]  14371       14545.94
## [15,]  14349       14564.79
## [16,]  14418       14583.64
## [17,]  14544       14602.50
## [18,]  14848       14621.35
## [19,]  14958       14640.20
## [20,]  14875       14659.05
## [21,]  15247       14677.91
## [22,]  15542       14696.76
## [23,]  15737       14715.61
## [24,]  15731       14734.47
## [25,]  14979       14753.32
## [26,]  15274       14772.17
## [27,]  15062       14791.02
## [28,]  14751       14809.88
## [29,]  14969       14828.73
## [30,]  15026       14847.58
forecast::accuracy(ts(hasil[-1]), head(test.ts, n=30))
##                ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
## Test set 136.6796 370.6676 265.8283 0.8640157 1.765892 0.7839405  1.626279

Dari hasil akurasi tersebut dapat dilihat MAPE dengan nilai 1.765892 yang terbilang kecil artinya model memiliki akurasi yang baik untuk meramalkan data. Hal ini dapat dimaknai bahwa data 240 periode awal atau 20 tahun terakhir dapat meramalkan dengan baik data ujinya sebanyak 30 periode dengan modle ARIMA(2,1,2)

sekian untuk analisis kali ini dan terimakasih.