G1401211002 DIVA NISFU MUSTIKA
Packages
library(ggplot2)
library(tsibble)##
## Attaching package: 'tsibble'## The following objects are masked from 'package:base':
##
## intersect, setdiff, unionlibrary(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(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':
##
## tarlibrary(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':
##
## identifylibrary(graphics)
library(readxl)Data
Data yang digunakan adalah harga jual gas bumi dan setiap periode data ditambahkan dengan angka 162.
data <- read_excel("data.xlsx")
head(data)## # A tibble: 6 × 2
## Periode Data
## <dbl> <dbl>
## 1 1 3980
## 2 2 4039
## 3 3 3878
## 4 4 3809
## 5 5 3829
## 6 6 3798dim(data)## [1] 151 2#Menghapus data periode
data <- data$DataMengubah data menjadi data time series
data.ts <- ts(data)Eksplorasi Data
Plot Data Penuh
plot.ts(data.ts, lty=1, xlab="waktu", ylab="harga jual gas bumi", main="Plot Data Harga Jual Gas Bumi")
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%.
Pembagian Data
Data kemudian dibagi menjadi data latih dan data uji. Pembagian kali ini dilakukan dengan proporsi / perbandingan, yaitu 80:20.
data.train <- data[1:120]
train.ts <- ts(data.train)
data.test <- data[121:151]
test.ts <- ts(data.test)Plot Data Latih
train.ts<-ts(data.train)
plot.ts(train.ts, lty=1, xlab="waktu", ylab="harga jual gas bumi", main="Plot Harga Jual Gas Bumi Train")
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
test.ts<-ts(data.test)
plot.ts(test.ts, lty=1, xlab="waktu", ylab="harga jual gas bumi", main="Plot Harga Jual Gas Bumi 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.5599, Lag order = 4, p-value = 0.344
## 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.344 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,10,by=1))
#Nilai Rounded Lambda
lambda <- bc$x[which.max(bc$y)]
lambda## [1] 2.333333#SK
bc$x[bc$y > max(bc$y) - 1/2 * qchisq(.95,1)]## [1] 1.555556 1.666667 1.777778 1.888889 2.000000 2.111111 2.222222 2.333333
## [9] 2.444444 2.555556 2.666667 2.777778 2.888889 3.000000 3.111111Plot Boxcox menunjukkan nilai rounded value (\(\lambda\)) optimum sebesar 2,333 dan pada selang kepercayaan 95% nilai memiliki batas bawah 1.555556 dan batas atas 3.1111117. 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 Harga Jual Gas Bumi", main="Plot Difference Harga Jual Gas Bumi")
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 = -4.5242, 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 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).
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 o o o o x o o o o o o o x
## 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 o o o o o o o o o o o o
## 4 x o o o o o o o o o o o o o
## 5 x o x x x o o o o o o o o o
## 6 o x x x x o o o o o o o o o
## 7 x x x x o o x o o o o o o oIdentifikasi model menggunakan plot EACF dilakukan dengan melihat ujung segitiga pada pola segitiga nol. Dalam hal ini model tentatif yang terbentuk adalah ARIMA(0,1,1), ARIMA(1,1,1), 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=1418.07## Series: train.diff
## ARIMA(0,1,1)
##
## Coefficients:
## ma1
## -1.0000
## s.e. 0.0229
##
## sigma^2 = 9080: log likelihood = -707.04
## AIC=1418.07 AICc=1418.18 BIC=1423.61
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 0.5237226 94.48326 72.30748 Inf Inf 0.6249841 -0.1848378lmtest::coeftest(model1.da) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -1.000000 0.022854 -43.756 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(1,1,0)
model2.da=Arima(train.diff, order=c(1,1,0),method="ML")
summary(model2.da) #AIC=1470.5## Series: train.diff
## ARIMA(1,1,0)
##
## Coefficients:
## ar1
## -0.5662
## s.e. 0.0766
##
## sigma^2 = 14696: log likelihood = -733.25
## AIC=1470.5 AICc=1470.61 BIC=1476.04
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -1.435216 120.2047 89.1368 Inf Inf 0.770447 -0.249072lmtest::coeftest(model2.da) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.566243 0.076567 -7.3954 1.41e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(1,1,1)
model3.da=Arima(train.diff, order=c(1,1,1),method="ML")
summary(model3.da) #AIC=1415.86## Series: train.diff
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## -0.1886 -1.000
## s.e. 0.0910 0.024
##
## sigma^2 = 8809: log likelihood = -704.93
## AIC=1415.86 AICc=1416.07 BIC=1424.17
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.53219 92.66424 71.10624 Inf Inf 0.6146013 -0.01374354lmtest::coeftest(model3.da) #seluruh parameter signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.188595 0.090952 -2.0736 0.03812 *
## ma1 -0.999994 0.024011 -41.6468 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1ARIMA(2,1,1)
model4.da=Arima(train.diff, order=c(2,1,1),method="ML")
summary(model4.da) #AIC=1417.25## Series: train.diff
## ARIMA(2,1,1)
##
## Coefficients:
## ar1 ar2 ma1
## -0.2026 -0.0727 -1.0000
## s.e. 0.0924 0.0930 0.0247
##
## sigma^2 = 8828: log likelihood = -704.63
## AIC=1417.25 AICc=1417.61 BIC=1428.34
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.636641 92.36665 71.00341 Inf Inf 0.6137125 0.004252496lmtest::coeftest(model4.da) #ada parameter yang tidak signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.202601 0.092398 -2.1927 0.02833 *
## ar2 -0.072689 0.092966 -0.7819 0.43428
## ma1 -1.000000 0.024664 -40.5442 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Berdasarkan pendugaan parameter di atas, nilai AIC terkecil dimiliki oleh model ARIMA(1,1,1) dan juga seluruhnya signifikan sehingga model yang dipilih adalah model ARIMA(1,1,1).
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 <- model3.da$residuals
par(mfrow=c(1,1))
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 tidak sama menandakan bahwa sisaan memiliki ragam yang heterogen. Plot ACF dan PACF sisaan ARIMA(1,1,1) juga 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
ks.test(sisaan.da,"pnorm") #tak tolak H0 > sisaan menyebar normal##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: sisaan.da
## D = 0.50419, p-value < 2.2e-16
## alternative hypothesis: two-sidedSelain 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 2.2e-16 yang kurang dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak menyebar normal. Hal ini tidak 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.023049, df = 1, p-value = 0.8793Selanjutnya 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.04117 yang lebih kecil dari taraf nyata 5% sehingga tolak \(H_0\) dan menandakan bahwa sisaan tidak saling bebas.
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.8471 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa sisaan saling bebas.Hal ini berbeda 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 = 0.0013599, df = 1, p-value = 0.9706Hipotesis 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.6754 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.17964, df = 118, p-value = 0.8577
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -15.35808 18.42246
## sample estimates:
## mean of x
## 1.53219Terakhir, 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.9521 yang lebih besar dari taraf nyata 5% sehingga tak tolak \(H_0\) dan menandakan bahwa nilai tengah sisaan sama dengan nol.
Overfitting
Tahapan selanjutnya adalah overfitting dilakukan dengan menaikkan orde AR(p) dan MA(q) dari model ARIMA(1,1,1) untuk melihat apakah terdapat model lain yang lebih baik dari model saat ini. Kandidat model overfitting adalah ARIMA(1,1,2) dan ARIMA(2,1,1).
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.0077 -1.2076 0.2076
## s.e. 0.3606 0.3494 0.3485
##
## sigma^2 = 8853: log likelihood = -704.77
## AIC=1417.54 AICc=1417.89 BIC=1428.62
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.609121 92.49418 71.0488 Inf Inf 0.6141047 -0.0009336968lmtest::coeftest(model.overfit1) #ar1 tidak signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.0076959 0.3606071 0.0213 0.9829733
## ma1 -1.2075707 0.3493580 -3.4565 0.0005472 ***
## ma2 0.2075831 0.3485135 0.5956 0.5514262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model ARIMA(2,1,1)
model.overfit2=Arima(train.diff, order=c(2,1,1),method="ML")
summary(model.overfit2) #798.74## Series: train.diff
## ARIMA(2,1,1)
##
## Coefficients:
## ar1 ar2 ma1
## -0.2026 -0.0727 -1.0000
## s.e. 0.0924 0.0930 0.0247
##
## sigma^2 = 8828: log likelihood = -704.63
## AIC=1417.25 AICc=1417.61 BIC=1428.34
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.636641 92.36665 71.00341 Inf Inf 0.6137125 0.004252496lmtest::coeftest(model.overfit2) #ma3 tidak signifikan##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.202601 0.092398 -2.1927 0.02833 *
## ar2 -0.072689 0.092966 -0.7819 0.43428
## ma1 -1.000000 0.024664 -40.5442 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1model yang dipilih adalah model awal, yaitu ARIMA(1,1,1) # Peramalan
Peramalan dilakukan menggunakan fungsi forecast() .
Contoh peramalan berikut ini dilakukan untuk 31 hari ke depan.
#---FORECAST---#
ramalan.da <- forecast::forecast(model4.da, h = 31)
ramalan.da## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 121 7.362690 -113.5582 128.2836 -177.5698 192.2952
## 122 -9.361101 -132.5376 113.8154 -197.7433 179.0211
## 123 -15.884931 -139.0958 107.3259 -204.3196 172.5498
## 124 -13.347562 -136.6010 109.9058 -201.8474 175.1523
## 125 -13.387423 -136.6395 109.8647 -201.8852 175.1104
## 126 -13.563786 -136.8151 109.6875 -202.0604 174.9328
## 127 -13.525158 -136.7767 109.7264 -202.0222 174.9719
## 128 -13.520164 -136.7718 109.7314 -202.0172 174.9769
## 129 -13.523984 -136.7756 109.7276 -202.0210 174.9730
## 130 -13.523573 -136.7751 109.7280 -202.0206 174.9734
## 131 -13.523379 -136.7750 109.7282 -202.0204 174.9736
## 132 -13.523448 -136.7750 109.7281 -202.0205 174.9736
## 133 -13.523448 -136.7750 109.7281 -202.0205 174.9736
## 134 -13.523443 -136.7750 109.7281 -202.0205 174.9736
## 135 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 136 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 137 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 138 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 139 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 140 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 141 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 142 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 143 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 144 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 145 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 146 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 147 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 148 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 149 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 150 -13.523444 -136.7750 109.7281 -202.0205 174.9736
## 151 -13.523444 -136.7750 109.7281 -202.0205 174.9736data.ramalan.da <- ramalan.da$mean
plot(ramalan.da)
Berdasarkan hasil plot ramalan di atas, dapat dilihat bahwa ramalan ARIMA(1,1,1) 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
hasil## Time Series:
## Start = 120
## End = 151
## Frequency = 1
## [1] 2355.000 2362.363 2353.002 2337.117 2323.769 2310.382 2296.818 2283.293
## [9] 2269.773 2256.249 2242.725 2229.202 2215.678 2202.155 2188.631 2175.108
## [17] 2161.584 2148.061 2134.538 2121.014 2107.491 2093.967 2080.444 2066.920
## [25] 2053.397 2039.873 2026.350 2012.827 1999.303 1985.780 1972.256 1958.733#has.1 sama hasilnya dengan: cumsum(c(pt_1,hasil.forc.Diff))
ts.plot(train.ts,hasil)
perbandingan.da<-matrix(data=c(head(test.ts, n=31), hasil[-1]),
nrow = 31, ncol = 2)
colnames(perbandingan.da)<-c("Aktual","Hasil Forecast")
perbandingan.da## Aktual Hasil Forecast
## [1,] 2356 2362.363
## [2,] 2302 2353.002
## [3,] 2272 2337.117
## [4,] 2324 2323.769
## [5,] 2269 2310.382
## [6,] 2299 2296.818
## [7,] 2284 2283.293
## [8,] 2279 2269.773
## [9,] 2326 2256.249
## [10,] 2351 2242.725
## [11,] 2335 2229.202
## [12,] 2360 2215.678
## [13,] 2269 2202.155
## [14,] 2275 2188.631
## [15,] 2226 2175.108
## [16,] 2208 2161.584
## [17,] 2191 2148.061
## [18,] 2149 2134.538
## [19,] 2126 2121.014
## [20,] 2152 2107.491
## [21,] 2187 2093.967
## [22,] 2141 2080.444
## [23,] 2089 2066.920
## [24,] 2107 2053.397
## [25,] 2031 2039.873
## [26,] 1984 2026.350
## [27,] 1983 2012.827
## [28,] 1981 1999.303
## [29,] 2016 1985.780
## [30,] 1951 1972.256
## [31,] 1951 1958.733accuracy(ts(hasil[-1]), head(test.ts, n=31))## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 24.68389 56.2917 43.5358 1.07675 1.958092 0.714279 1.384581Dari hasil diperoleh nilai MAPE 1.958092 sangat baik untuk menggambarkan model deret waktu