Metode Peramalan Deret Waktu

Tugas Individu

Library

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
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(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
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(TTR)
## Warning: package 'TTR' was built under R version 4.3.3
library(aTSA)
## Warning: package 'aTSA' was built under R version 4.3.2
## 
## 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)
library(readxl)
## Warning: package 'readxl' was built under R version 4.3.2

Data

data <- read.csv("C:/Users/Muhammad Rizqa Salas/Downloads/Semester 5/MPDW/Project/Data Air MPDW.csv")
dim(data)
## [1] 100   2
data.ts <- ts(data)
str(data.ts)
##  Time-Series [1:100, 1:2] from 1 to 100: 1 2 3 4 5 6 7 8 9 10 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr [1:2] "Periode" "Air"

Eksplorasi Data

Plot Data Penuh

data.ts <- ts(data$Air) 
ts.plot(data.ts, xlab="Time Period", ylab="Air", main="Time Series Plot")

mean(data.ts)
## [1] 185.85

Plot deret waktu di atas menunjukkan bahwa data cenderung tidak stasioner dalam rataan, ditandai dengan data yang tidak menyebar di sekitar nilai tengahnya (185.5). Namun, plot tersebut menunjukan indikasi data stasioner dalam ragam, ditandai dengan lebar pita yang cenderung sama.

lattice::densityplot(as.vector(data.ts))

Split Data

train <- data[1:80,]
test <- data[81:100,]
train.ts <- ts(train$Air)
test.ts <- ts(test$Air)

Plot Data Latih

plot(train.ts, col="blue",main="Plot data training")

mean(train.ts)
## [1] 187.6625

Plot deret waktu data train di atas menunjukkan bahwa data cenderung tidak stasioner dalam rataan, ditandai dengan data yang tidak menyebar di sekitar nilai tengahnya (187.6625). Plot tersebut juga menunjukan indikasi data stasioner dalam ragam, ditandai dengan lebar pita yang cenderung sama.

Plot Data Uji

plot(test.ts, col="blue",main="Plot data testing")

mean(test.ts)
## [1] 178.6

Uji Stasioneritas Data

Plot ACF

acf(train.ts)

Berdasarkan plot ACF, terlihat bahwa plot ACF pada data tersebut cenderung tails off slowly menandakan data tidak stasioner dalam rataan

Uji ADF

tseries::adf.test(train.ts)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.ts
## Dickey-Fuller = -4.0619, Lag order = 4, p-value = 0.01112
## 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.01112 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa data tidak stasioner dalam rataan. Hal ini sesuai dengan hasil plot ACF

Plot Box-Cox

index <- seq(1:80)
bc = boxcox(train.ts~index, lambda = seq(-5,20,by=0.1))

#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda
## [1] 2
#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]
##  [1] 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3
## [20] 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6

Gambar di atas menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 0.95 dan pada selang kepercayaan 95% nilai memiliki batas bawah 0.5 dan batas atas 3.6. Selang tersebut memuat nilai satu sehingga dapat dikatakan bahwa data keseluruhan stasioner dalam ragam.

Dikarenakan data train telah stasioner dalam rataan maupun ragam, maka tidak diperlukan penanganan ketidakstasioneran pada data sehingga dapat dilanjutkan dengan identifikasi model.

Penanganan Ketidakstasioneran dalam Rataan

Menggunakan differencing

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

acf(train.diff)

tseries::adf.test(train.diff)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.diff
## Dickey-Fuller = -3.7856, Lag order = 4, p-value = 0.02387
## alternative hypothesis: stationary

sudah stasioner p-value < 0.05

Identifikasi Model

Plot ACF

acf(train.diff, 20)

Berdasarkan plot tersebut, terlihat bahwa plot ACF cenderung tails off pada lag ke 3, sehingga jika plot PACF dianggap tails of, maka model tentatifnya adalah ARIMA(0,0,3). Namun plot tersebut juga menunjukan adanya pola tails off, sehingga jika plot PACF dianggap tails of, maka model tentatifnya adalah ARIMA(0,0,0). Dikarenakan perbedaan tersebut, maka dilakukan pembandingan keduanya untuk melihat model mana yang merupakan model terbaik.

Plot PACF

pacf(train.diff, 20)

Berdasarkan plot tersebut, terlihat bahwa plot PACF cenderung cuts off pada lag ke 1, sehingga model tentatifnya adalah ARIMA(1,1,2) dan 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 x x o o o x x x x x  x  x  x 
## 1 x x o o o o o o o o o  o  o  o 
## 2 x x o o o o x o o o o  o  o  o 
## 3 x o o o o o o o o o o  o  o  o 
## 4 x x o x o o o o o o o  o  o  o 
## 5 x o o o o x o o o o o  o  o  o 
## 6 o x o o o x o o o o o  o  o  o 
## 7 o x o o o x o o x 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 pada plot EACF adalah ARIMA(1,1,2), ARIMA(2,1,1), ARIMA(0,1,3), ARIMA(1,1,3), ARIMA(2,1,2)

Sehingga didapatkan 5 model tentatif yaitu, ARIMA(1,1,2), ARIMA(2,1,1), ARIMA(0,1,3), ARIMA(1,1,3), dan ARIMA(2,1,2).

Pendugaan Parameter Model Tentatif

ARIMA(1,1,2)

model1.da=Arima(train.diff, order=c(1,1,2),method="ML")
summary(model1.da) #AIC = 530.56
## Series: train.diff 
## ARIMA(1,1,2) 
## 
## Coefficients:
##          ar1      ma1     ma2
##       0.1550  -0.6732  0.5009
## s.e.  0.2379   0.2094  0.1258
## 
## sigma^2 = 49.05:  log likelihood = -261.28
## AIC=530.56   AICc=531.11   BIC=539.99
## 
## Training set error measures:
##                      ME     RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set 0.01421396 6.823736 5.575637 -Inf  Inf 0.9079325 0.004119044
lmtest::coeftest(model1.da)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.15495    0.23791  0.6513  0.514851    
## ma1 -0.67324    0.20942 -3.2148  0.001305 ** 
## ma2  0.50092    0.12584  3.9806 6.875e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(2,1,1)

model2.da=Arima(train.diff, order=c(2,1,1),method="ML")
summary(model2.da) #AIC = 532.3
## Series: train.diff 
## ARIMA(2,1,1) 
## 
## Coefficients:
##           ar1     ar2      ma1
##       -0.1050  0.3278  -0.3596
## s.e.   0.3264  0.1844   0.3263
## 
## sigma^2 = 50.31:  log likelihood = -262.15
## AIC=532.3   AICc=532.84   BIC=541.72
## 
## Training set error measures:
##                      ME    RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set 0.02474834 6.91107 5.494428 -Inf  Inf 0.8947085 -0.01043149
lmtest::coeftest(model2.da)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value Pr(>|z|)  
## ar1 -0.10505    0.32645 -0.3218  0.74761  
## ar2  0.32781    0.18441  1.7776  0.07546 .
## ma1 -0.35964    0.32635 -1.1020  0.27046  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(0,1,3)

model3.da=Arima(train.diff, order=c(0,1,3),method="ML")
summary(model3.da) #AIC = 530.64 
## Series: train.diff 
## ARIMA(0,1,3) 
## 
## Coefficients:
##           ma1     ma2     ma3
##       -0.5253  0.4216  0.0590
## s.e.   0.1084  0.1093  0.1039
## 
## sigma^2 = 49.1:  log likelihood = -261.32
## AIC=530.64   AICc=531.19   BIC=540.07
## 
## Training set error measures:
##                      ME     RMSE     MAE  MPE MAPE      MASE       ACF1
## Training set 0.01411674 6.827184 5.57133 -Inf  Inf 0.9072311 0.00988099
lmtest::coeftest(model3.da)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.525338   0.108415 -4.8456 1.262e-06 ***
## ma2  0.421569   0.109311  3.8566  0.000115 ***
## ma3  0.058997   0.103902  0.5678  0.570159    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ARIMA(1,1,3)

model4.da=Arima(train.diff, order=c(1,1,3),method="ML")
summary(model4.da) #AIC = 532.56
## Series: train.diff 
## ARIMA(1,1,3) 
## 
## Coefficients:
##          ar1      ma1     ma2      ma3
##       0.1809  -0.6992  0.5157  -0.0128
## s.e.  0.5920   0.5857  0.3362   0.2730
## 
## sigma^2 = 49.71:  log likelihood = -261.28
## AIC=532.56   AICc=533.4   BIC=544.35
## 
## Training set error measures:
##                      ME     RMSE      MAE  MPE MAPE     MASE        ACF1
## Training set 0.01424017 6.823663 5.576751 -Inf  Inf 0.908114 0.004707977
lmtest::coeftest(model4.da)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1  0.180852   0.592012  0.3055   0.7600
## ma1 -0.699245   0.585747 -1.1938   0.2326
## ma2  0.515677   0.336213  1.5338   0.1251
## ma3 -0.012814   0.273042 -0.0469   0.9626

ARIMA(2,1,2)

model5.da=Arima(train.diff, order=c(2,1,2),method="ML")
summary(model5.da) #AIC = 1227.01
## Series: train.diff 
## ARIMA(2,1,2) 
## 
## Coefficients:
##          ar1     ar2      ma1     ma2
##       0.1536  0.0270  -0.6730  0.4825
## s.e.  0.2383  0.2215   0.2129  0.2060
## 
## sigma^2 = 49.7:  log likelihood = -261.27
## AIC=532.55   AICc=533.38   BIC=544.33
## 
## Training set error measures:
##                      ME     RMSE     MAE  MPE MAPE      MASE        ACF1
## Training set 0.01445081 6.823303 5.58623 -Inf  Inf 0.9096575 0.008623022
lmtest::coeftest(model5.da)
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)   
## ar1  0.153608   0.238252  0.6447 0.519102   
## ar2  0.027028   0.221546  0.1220 0.902902   
## ma1 -0.673049   0.212886 -3.1615 0.001569 **
## ma2  0.482534   0.206041  2.3419 0.019184 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan pendugaan parameter di atas, nilai AIC terkecil dimiliki oleh model ARIMA(0,1,3) dan parameter model ARIMA(0,1,3) juga seluruhnya signifikan kecuali ma3 sehingga model yang dipilih adalah model ARIMA(0,1,3).

Analisis Sisaan

Eksplorasi Sisaan

sisaan.da <- model3.da$residuals 
qqnorm(sisaan.da) 
qqline(sisaan.da, col = "blue", lwd = 2) 

plot(c(1:length(sisaan.da)),sisaan.da) 

acf(sisaan.da) 

pacf(sisaan.da)

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

Uji Formal Sisaan

Asumsi Normalitas

ks.test(sisaan.da,"pnorm") 
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan.da
## D = 0.40945, p-value = 1.723e-12
## alternative hypothesis: two-sided

\(H_0\) : Sisaan menyebar normal

\(H_1\) : Sisaan tidak menyebar normal

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

Asumsi Sisaan Saling Bebas

Box.test(sisaan.da, type = "Ljung") 
## 
##  Box-Ljung test
## 
## data:  sisaan.da
## X-squared = 0.0080097, df = 1, p-value = 0.9287

\(H_0\) : Sisaan saling bebas

\(H_1\) : Sisaan tidak tidak saling bebas

Berdasarkan uji Ljung-Box tersebut, didapat p-value sebesar 0.9287 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan tidak cukup bukti untuk mengatakan bahwa sisaan saling bebas. Menandakan bahwa asumsi tidak terpenuhi.

Asumsi Heteroskedastisitas

Box.test((sisaan.da)^2, type = "Ljung") 
## 
##  Box-Ljung test
## 
## data:  (sisaan.da)^2
## X-squared = 0.084187, df = 1, p-value = 0.7717

\(H_0\) : Ragam sisaan homogen

\(H_1\) : Ragam sisaan tidak homogen

Berdasarkan uji Ljung-Box terhadap sisaan kuadrat tersebut, didapat p-value sebesar 0.7717 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan tidak cukup bukti untuk mengatakan bahwa ragam sisaan homogen. Menandakan bahwa asumsi terpenuhi. Hal ini sesuai dengan hasil eksplorasi.

Asumsi Nilai Tengah Sisaan = 0

t.test(sisaan.da, mu = 0, conf.level = 0.95) 
## 
##  One Sample t-test
## 
## data:  sisaan.da
## t = 0.018262, df = 78, p-value = 0.9855
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.524857  1.553091
## sample estimates:
##  mean of x 
## 0.01411674

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

Overfitting

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

model3a.da=Arima(train.ts, order=c(1,1,3),method="ML")
summary(model3a.da) #AIC = 1220.87
## Series: train.ts 
## ARIMA(1,1,3) 
## 
## Coefficients:
##          ar1      ma1     ma2     ma3
##       0.6897  -0.3849  0.4379  0.1437
## s.e.  0.1226   0.1431  0.1021  0.1089
## 
## sigma^2 = 43.29:  log likelihood = -259.57
## AIC=529.13   AICc=529.95   BIC=540.98
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.1900154 6.370927 5.005916 0.2143808 2.796443 0.6547473
##                    ACF1
## Training set 0.02305274
lmtest::coeftest(model3a.da)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1  0.68974    0.12260  5.6259 1.845e-08 ***
## ma1 -0.38493    0.14309 -2.6901  0.007143 ** 
## ma2  0.43790    0.10210  4.2889 1.796e-05 ***
## ma3  0.14367    0.10891  1.3192  0.187103    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model3b.da=Arima(train.ts, order=c(0,1,4),method="ML")
summary(model3b.da) #AIC = 1221.42
## Series: train.ts 
## ARIMA(0,1,4) 
## 
## Coefficients:
##          ma1     ma2     ma3     ma4
##       0.2756  0.7671  0.3252  0.4532
## s.e.  0.1012  0.1440  0.1768  0.1179
## 
## sigma^2 = 42.08:  log likelihood = -258.63
## AIC=527.27   AICc=528.09   BIC=539.12
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE     MAPE      MASE       ACF1
## Training set 0.234005 6.281192 4.927955 0.1826049 2.735506 0.6445504 0.07164885
lmtest::coeftest(model3b.da)
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ma1  0.27558    0.10120  2.7231  0.006467 ** 
## ma2  0.76706    0.14398  5.3276 9.953e-08 ***
## ma3  0.32521    0.17681  1.8393  0.065874 .  
## ma4  0.45325    0.11790  3.8442  0.000121 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan kedua model hasil overfitting di atas, model ARIMA(1,1,3) dan ARIMA(0,1,4) memiliki AIC yang lebih kecil dibandingkan dengan model ARIMA(0,1,3) parameter kedua model juga banyak yang signifikan. Oleh karena itu, model ARIMA(0,1,4) akan digunakan untuk melakukan peramalan karena memiliki AIC terkecil.

Peramalan

ramalan.da <- forecast::forecast(model3b.da, h = 20) 
ramalan.da
##     Point Forecast    Lo 80    Hi 80     Lo 95    Hi 95
##  81       232.0085 223.6949 240.3222 219.29389 244.7232
##  82       237.6277 224.1526 251.1027 217.01938 258.2360
##  83       240.4698 218.7913 262.1484 207.31537 273.6243
##  84       242.0127 212.7300 271.2954 197.22868 286.7968
##  85       242.0127 204.4954 279.5300 184.63494 299.3905
##  86       242.0127 197.7677 286.2577 174.34583 309.6796
##  87       242.0127 191.9359 292.0896 165.42677 318.5987
##  88       242.0127 186.7157 297.3098 157.44318 326.5823
##  89       242.0127 181.9475 302.0780 150.15084 333.8746
##  90       242.0127 177.5309 306.4946 143.39627 340.6292
##  91       242.0127 173.3980 310.6274 137.07558 346.9499
##  92       242.0127 169.5003 314.5251 131.11456 352.9109
##  93       242.0127 165.8017 318.2238 125.45801 358.5674
##  94       242.0127 162.2744 321.7510 120.06355 363.9619
##  95       242.0127 158.8968 325.1287 114.89782 369.1276
##  96       242.0127 155.6511 328.3744 109.93396 374.0915
##  97       242.0127 152.5230 331.5024 105.15003 378.8754
##  98       242.0127 149.5007 334.5248 100.52775 383.4977
##  99       242.0127 146.5740 337.4514  96.05178 387.9737
## 100       242.0127 143.7344 340.2910  91.70905 392.3164
data.ramalan.da <- ramalan.da$mean
plot(ramalan.da)

Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa ramalan ARIMA(0,1,4) cenderung tidak stabil hingga akhir periode ditandai dengan lebar pita yang semakin mengecil. Selanjutnya, dapat dicari nilai akurasi antara hasil ramalan dengan data uji sebagai berikut.

plot(train.ts, xlim = c(1, length(train.ts) + length(test.ts)), ylim = range(c(train.ts, test.ts, data.ramalan.da)), col = "black", main = "Data Latih, Uji, dan Ramalan", ylab = "Nilai", xlab = "Waktu")
lines(seq(length(train.ts) + 1, length(train.ts) + length(test.ts)), test.ts, col = "red")
lines(seq(length(train.ts) + 1, length(train.ts) + length(test.ts)), data.ramalan.da, col = "blue", lty = 2)
legend("topright", legend = c("Data Latih", "Data Uji", "Ramalan"), col = c("black", "red", "blue"), lty = c(1, 1, 2))

perbandingan.da <- matrix(data=c(head(test.ts, n=20), data.ramalan.da), 
                          nrow = 20, ncol = 2)
colnames(perbandingan.da) <- c("Aktual", "Hasil Forecast")
perbandingan.da
##       Aktual Hasil Forecast
##  [1,]    227       232.0085
##  [2,]    227       237.6277
##  [3,]    223       240.4698
##  [4,]    214       242.0127
##  [5,]    200       242.0127
##  [6,]    188       242.0127
##  [7,]    170       242.0127
##  [8,]    155       242.0127
##  [9,]    142       242.0127
## [10,]    135       242.0127
## [11,]    137       242.0127
## [12,]    135       242.0127
## [13,]    180       242.0127
## [14,]    158       242.0127
## [15,]    165       242.0127
## [16,]    175       242.0127
## [17,]    178       242.0127
## [18,]    185       242.0127
## [19,]    188       242.0127
## [20,]    190       242.0127
accuracy(ts(data.ramalan.da), head(test.ts, n=20))
##                 ME     RMSE      MAE       MPE     MAPE     ACF1 Theil's U
## Test set -62.61612 69.60504 62.61612 -38.85718 38.85718 0.798054  4.826294

Kesimpulan

Berdasarkan hasil identifikasi model, diperoleh model terbaik yaitu menggunakan model ARIMA(0,1,4). Namun, hasil peramalan menunjukkan perbedaan yang cukup signifikan, terlihat dari semakin menyempitnya lebar pita hasil peramalan seiring berjalannya waktu. Hal ini juga dibuktikan dengan nilai MAPE 3.886%, model yang dihasilkan cukup baik untuk digunakan dikarenakan nilai MAPE yang masih bisa dikatakan cukup baik. Namun tentu saja masih dapat ditingkatkan melalui optimasi lebih lanjut atau penyesuaian model lainnya sehingga akurasi model dapat lebih ditingkatkan.