Tugas Praktikum 8

Khusnia Nurul Khikmah (G1501211049)

4/4/2022

Soal

Lakukan pemodelan deret waktu untuk data musiman yang terdapat pada file latihan Praktikum 8.xlsx (Ada 3 gugus data) Catatan :

  1. Gunakan Tahapan Identifikasi Model ARIMA Musiman (plot, identifikasi, pendugaan dan seleksi model)

  2. Gunakan \(80\%\) data yang ada sebagai data training, sisanya \(20\%\) data untuk testing. Hitung nilai akurasinya menggunakan MAD,MSE, dan MAPE.

  3. Lakukan forecasting untuk 6 periode data kedepan.

Jawaban Gugus Data 1

Membaca data

setwd("D:/S2/IPB/STA542 Analisis Deret Waktu/8/Res")
data <- read.xlsx("latihan Praktikum 8.xlsx", "data1", header=TRUE)
head(data)
##   tahun  Yt
## 1  1903 221
## 2  1904 279
## 3  1905 276
## 4  1906 111
## 5  1907 503
## 6  1908 360

Forming data time series

Di sini dilakukan perubahan data Yt dari data.frame menjadi data ts atau time series dengan fungsi ts.

data.ts <- ts(data, start=1903)

Visualisasi data time series

Data yang digunakan terdiri dari \(114\) observasi atau yang diukur dalam \(114\) tahun. Hasil dari visualisasi data tersaji pada plot berikut:

ts.plot(data.ts[,"Yt"], type="l", ylab="Yt", col="blue")
title(main = "Time Series Plot of Yt", cex.sub = 0.8)
points(data.ts[,"Yt"], pch = 20, col = "blue")

Mengecek seasonality data time series

seasonplot(data.ts[,"Yt"], 5, main="Yt", ylab="YT", year.labels = TRUE, col=rainbow(18))

seasonplot(data.ts[,"Yt"], 8, main="Yt", ylab="YT", year.labels = TRUE, col=rainbow(18))

seasonplot(data.ts[,"Yt"], 10, main="Yt", ylab="YT", year.labels = TRUE, col=rainbow(18))

seasonplot(data.ts[,"Yt"], 12, main="Yt", ylab="YT", year.labels = TRUE, col=rainbow(18))

Berdasarkan percobaan plot seasonal di atas dapat diketahui bahwa tidak menunjukkan pola data yang seasonal maka akan dilakukan analisis model ARIMA reguler.

Splitting data

Data yang akan digunakan untuk menentukan model terbaik berdasarkan kandidat model akan dilakukan splitting data terlebih dahulu yaitu dengan membagi data menjadi dua bagian bagian pertama yaitu sebagai data training sebesar \(80\%\) atau \(91\) 0bservasi dan bagian kedua sebagai data testing sebesar \(20\%\) atau \(23\) observasi. Splitting data ini menggunakan fungsi subset.

# Training set
train.Yt.ts = subset(data.ts[,"Yt"], start=1,end=91)

# Testing set
test.Yt.ts = subset(data.ts[,"Yt"], start=92,end=114)

Visualisasi data training yang akan digunakan untuk mencari model terbaik

ts.plot(train.Yt.ts, col="blue", ylab = "Yt", xlab = "Year")
title(main = "Train Time Series Plot of Yt", cex.sub = 0.8)
points(train.Yt.ts, pch = 20, col = "blue")

Visualisasi data testing yang akan digunakan untuk mencari model terbaik

ts.plot(test.Yt.ts, col="blue", ylab = "Yt", xlab = "Year")
title(main = "Test Time Series Plot of Yt", cex.sub = 0.8)
points(test.Yt.ts, pch = 20, col = "blue")

Model ARIMA

Model ARIMA terbaik dari data Yt sebagai berikut:

1. Cek kestasioneran data

Mengidentifikasi kestasioneran data dilakukan dengan uji ADF.

Uji Augmented Dickey-Fuller (ADF) adalah salah satu alat identifikasi untuk menguji kestasioneran data dan memiliki hipotesis sebagai berikut:

\(H_{0}\) : Data tidak stasioner (unit roots mendekati 1)

\(H_{1}\) : Data stasioner (unit roots tidak mendekati 1)

adf.test(train.Yt.ts, alternative = "stationary", k=trunc((length(train.Yt.ts)-1)^(1/3)))
## 
##  Augmented Dickey-Fuller Test
## 
## data:  train.Yt.ts
## Dickey-Fuller = -3.2448, Lag order = 4, p-value = 0.08539
## alternative hypothesis: stationary

Berdasarkan hasil uji ADF dengan taraf signifikansi \(\alpha = 5% = 0.05\) di atas dapat diketahui bersama bahwa, \(p-value = 0.08539\) \(>\) \(\alpha = 0.05\) sehingga keputusan yang diambil adalah Tak Tolak \(H_{0}\) artinya data yang digunakan merupakan data yang tidak stasioner, sehingga dilakukan differencing.

Differencing

  • Differencing ordo 1
yt.diff1 <- diff(train.Yt.ts, difference=1)
adf.test(yt.diff1, alternative = "stationary", k=trunc((length(yt.diff1)-1)^(1/3)))
## Warning in adf.test(yt.diff1, alternative = "stationary", k =
## trunc((length(yt.diff1) - : p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  yt.diff1
## Dickey-Fuller = -5.9477, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Berdasarkan hasil uji ADF dengan taraf signifikansi \(\alpha = 5% = 0.05\) di atas dapat diketahui bersama bahwa, \(p-value = 0.01\) \(<\) \(\alpha = 0.05\) sehingga keputusan yang diambil adalah Tolak \(H_{0}\) artinya data yang telah didifferencing 1 kali merupakan data yang stasioner, sehingga bisa untuk dilakukan analisis lebih lanjut.

2. Identifikasi kandidat model

Identifikasi kandidat model diperoleh berdasarkan nilai \(p\) dan \(q\) dimana nilai \(d=1\)

acf(yt.diff1, lag.max = 30, col = "blue")

Berdasarkan plot ACF terlihat bahwa plot cut-off di lag 1, \(ARIMA(0,1,1)\).

pacf(yt.diff1, lag.max = 20, col = "blue")

Berdasarkan plot PACF di atas kandidat model yang diperoleh adalah \(ARIMA(2,1,0)\)

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

Berdasarkan hasil EACF di atas kandidat model yang diperoleh adalah \(ARIMA(2,1,1)\), \(ARIMA(1,1,2)\) dan \(ARIMA(3,1,2)\)

Karena data yang digunakan adalah data yang hasil differencing maka \(d=1\) dan dapat ditentukan beberapa kandidat model yaitu: \(ARIMA(0,1,1)\), \(ARIMA(2,1,0)\), \(ARIMA(2,1,1)\), \(ARIMA(1,1,2)\) dan \(ARIMA(3,1,2)\).

# ARIMA(0,1,1)
arima011 <- arima(yt.diff1, order = c(0,0,1), include.mean = TRUE, method = "ML")
arima011
## 
## Call:
## arima(x = yt.diff1, order = c(0, 0, 1), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##           ma1  intercept
##       -0.6816    18.3087
## s.e.   0.0979     6.3849
## 
## sigma^2 estimated as 34450:  log likelihood = -598.14,  aic = 1200.29
# ARIMA(2,1,0)
arima210 <- arima(yt.diff1, order = c(2,0,0), include.mean = TRUE, method = "ML")
arima210
## 
## Call:
## arima(x = yt.diff1, order = c(2, 0, 0), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##           ar1      ar2  intercept
##       -0.5357  -0.2165    18.4607
## s.e.   0.1027   0.1022    11.6234
## 
## sigma^2 estimated as 36866:  log likelihood = -601.04,  aic = 1208.07
# ARIMA(2,1,1)
arima211 <- arima(yt.diff1, order = c(2,0,1), include.mean = TRUE, method = "ML")
arima211
## 
## Call:
## arima(x = yt.diff1, order = c(2, 0, 1), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1     ar2      ma1  intercept
##       0.2805  0.1646  -1.0000    19.2583
## s.e.  0.1041  0.1039   0.0325     1.2234
## 
## sigma^2 estimated as 31052:  log likelihood = -594.92,  aic = 1197.85
# ARIMA(1,1,2)
arima112 <- arima(yt.diff1, order = c(1,0,2), include.mean = TRUE, method = "ML")
arima112
## 
## Call:
## arima(x = yt.diff1, order = c(1, 0, 2), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ma1     ma2  intercept
##       0.7037  -1.4174  0.4174    19.2082
## s.e.  0.1792   0.2217  0.2199     1.3159
## 
## sigma^2 estimated as 31112:  log likelihood = -594.93,  aic = 1197.86
# ARIMA(3,1,2)
arima312 <- arima(yt.diff1, order = c(3,0,2), include.mean = TRUE, method = "ML")
arima312
## 
## Call:
## arima(x = yt.diff1, order = c(3, 0, 2), include.mean = TRUE, method = "ML")
## 
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): NaNs produced
##          ar1     ar2     ar3      ma1     ma2  intercept
##       0.2804  0.1544  0.0336  -1.0051  0.0051    19.2416
## s.e.     NaN     NaN     NaN      NaN     NaN     1.2490
## 
## sigma^2 estimated as 31038:  log likelihood = -594.87,  aic = 1201.74

Model Terbaik

AICKandidatModel <- c(arima011$aic, arima210$aic, arima211$aic, arima112$aic, arima312$aic)
KandidatModelARIMA <- c("ARIMA(0,1,1)", "ARIMA(2,1,0)", "ARIMA(2,1,1)", "ARIMA(1,1,1)", "ARIMA(3,1,2)")
compmodelARIMA <- cbind(KandidatModelARIMA, AICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA
##   Kandidat Model        Nilai AIC
## 1   ARIMA(0,1,1)  1200.2873778929
## 2   ARIMA(2,1,0) 1208.07443005334
## 3   ARIMA(2,1,1) 1197.84904150017
## 4   ARIMA(1,1,1) 1197.85679424979
## 5   ARIMA(3,1,2) 1201.74203668152

Model terbaik diperoleh Berdasarkan nilai AIC terkecil dari kandidat model. Oleh karena itu, Model terbaik yang diperoleh yaitu \(ARIMA(2,1,1)\)

Model terbaik ini selanjutnya dilakukan uji hipotesis untuk mengetahui signifikansi parameternya, dimana:

\(H_{0}\) : parameter \(= 0\)

\(H_{1}\) : parameter \(\neq\) \(0\)

Pengambilan keputusan yaitu Tolak \(H_{0}\) jika \(|t_{hitung}| > t_{\frac{\alpha}{2},(n-1)}\)

T-hitung parameter dugaan model \(ARIMA(1,0,0)\) \[t_{hitung}=\frac{parameter estimasi}{SE parameter estimasi}\]

intercept_mu <- 19.2583/1.2234
intercept_mu
## [1] 15.74162
ar_phi1 <- 0.2805/0.1041
ar_phi1
## [1] 2.694524
ar_phi2 <- 0.1646/0.1039
ar_phi2
## [1] 1.584216
ma_theta1 <- -1.0000/0.0325
ma_theta1
## [1] -30.76923
prmtrdugaanarima211 <- c("mu", "phi(1)", "phi(2)", "theta(1)")
thitungprmtrarima211 <- c(abs(intercept_mu), abs(ar_phi1), abs(ar_phi2), abs(ma_theta1))
ttabel <- c("1.661961084", "1.661961084", "1.661961084")
keputusan <- c("Signifikan", "Signifikan", "Tidak Signifikan", "Signifikan")
tksignimodelarima211 <- cbind(prmtrdugaanarima211, thitungprmtrarima211, ttabel, keputusan)
## Warning in cbind(prmtrdugaanarima211, thitungprmtrarima211, ttabel, keputusan):
## number of rows of result is not a multiple of vector length (arg 3)
colnames(tksignimodelarima211) <- c("Parameter dugaan", "T-Hitung", "T-Tabel", "Keputusan")
tksignimodelarima211 <- as.data.frame(tksignimodelarima211)
tksignimodelarima211
##   Parameter dugaan         T-Hitung     T-Tabel        Keputusan
## 1               mu 15.7416217099886 1.661961084       Signifikan
## 2           phi(1) 2.69452449567723 1.661961084       Signifikan
## 3           phi(2)  1.5842155919153 1.661961084 Tidak Signifikan
## 4         theta(1) 30.7692307692308 1.661961084       Signifikan
lmtest::coeftest(arima211)
## 
## z test of coefficients:
## 
##            Estimate Std. Error  z value  Pr(>|z|)    
## ar1        0.280537   0.104133   2.6940  0.007059 ** 
## ar2        0.164620   0.103907   1.5843  0.113123    
## ma1       -1.000000   0.032499 -30.7704 < 2.2e-16 ***
## intercept 19.258323   1.223358  15.7422 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Berdasarkan tabel di atas maka pengambilan keputusan dari paramater dugaan dari model \(ARIMA(2,1,1)\) semua parameternya signifikan, kecuali \(\phi_{2}\).

Selanjutnya, model terbaik tersebut digunakan untuk peramalan.

Dibandingkan auto.arima

auto.arima(train.Yt.ts)
## Series: train.Yt.ts 
## ARIMA(0,1,1) with drift 
## 
## Coefficients:
##           ma1    drift
##       -0.6816  18.3086
## s.e.   0.0979   6.3849
## 
## sigma^2 estimated as 35233:  log likelihood=-598.14
## AIC=1202.29   AICc=1202.57   BIC=1209.79

Dibandingkan Arima

model1 <- Arima(yt.diff1,order=c(0,0,1))
summary(model1)
## Series: yt.diff1 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##           ma1     mean
##       -0.6816  18.3086
## s.e.   0.0979   6.3849
## 
## sigma^2 estimated as 35233:  log likelihood=-598.14
## AIC=1202.29   AICc=1202.57   BIC=1209.79
## 
## Training set error measures:
##                    ME    RMSE      MAE MPE MAPE      MASE       ACF1
## Training set 0.307445 185.607 131.3133 Inf  Inf 0.4972083 0.04781962
model2 <- Arima(yt.diff1,order=c(2,0,0))
summary(model2)
## Series: yt.diff1 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##           ar1      ar2     mean
##       -0.5357  -0.2164  18.4590
## s.e.   0.1027   0.1022  11.6236
## 
## sigma^2 estimated as 38137:  log likelihood=-601.04
## AIC=1210.07   AICc=1210.55   BIC=1220.07
## 
## Training set error measures:
##                    ME     RMSE     MAE  MPE MAPE      MASE        ACF1
## Training set 0.185087 192.0051 131.859 -Inf  Inf 0.4992745 -0.04017621
model3 <- Arima(yt.diff1,order=c(2,0,1))
summary(model3)
## Series: yt.diff1 
## ARIMA(2,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1     ar2      ma1     mean
##       0.2805  0.1646  -1.0000  19.2583
## s.e.  0.1041  0.1039   0.0325   1.2234
## 
## sigma^2 estimated as 32497:  log likelihood=-594.92
## AIC=1199.85   AICc=1200.56   BIC=1212.35
## 
## Training set error measures:
##                     ME    RMSE      MAE  MPE MAPE     MASE        ACF1
## Training set -4.670714 176.217 127.1697 -Inf  Inf 0.481519 -0.01401929
model4 <- Arima(yt.diff1,order=c(1,0,1))
summary(model4)
## Series: yt.diff1 
## ARIMA(1,0,1) with non-zero mean 
## 
## Coefficients:
##          ar1     ma1     mean
##       0.3315  -1.000  19.2873
## s.e.  0.1002   0.036   1.0480
## 
## sigma^2 estimated as 32912:  log likelihood=-596.17
## AIC=1200.33   AICc=1200.8   BIC=1210.33
## 
## Training set error measures:
##                     ME     RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set -5.700151 178.3676 130.3274 -Inf  Inf 0.4934752 -0.06305828
model5 <- Arima(yt.diff1,order=c(3,0,2))
summary(model5)
## Series: yt.diff1 
## ARIMA(3,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ar2     ar3      ma1      ma2     mean
##       -0.5829  0.4072  0.1328  -0.1349  -0.8647  19.2648
## s.e.   0.4180  0.1594  0.1291   0.4037   0.4033   1.2150
## 
## sigma^2 estimated as 33264:  log likelihood=-594.91
## AIC=1203.83   AICc=1205.19   BIC=1221.33
## 
## Training set error measures:
##                     ME     RMSE     MAE  MPE MAPE      MASE        ACF1
## Training set -4.727184 176.2008 127.227 -Inf  Inf 0.4817358 -0.01492458
AICKandidatModel <- c(model1$aic, model2$aic, model3$aic, model4$aic, model5$aic)
AICcKandidatModel <- c(model1$aicc, model2$aicc, model3$aicc, model4$aicc, model5$aicc)
BICKandidatModel <- c(model1$bic, model2$bic, model3$bic, model4$bic, model5$bic)
KandidatModelARIMA <- c("ARIMA(0,1,1)", "ARIMA(2,1,0)", "ARIMA(2,1,1)", "ARIMA(1,1,1)", "ARIMA(3,1,2)")
compmodelARIMA <- cbind(KandidatModelARIMA, AICKandidatModel, AICcKandidatModel, BICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC", "Nilai AICc", "Nilai BIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA
##   Kandidat Model        Nilai AIC       Nilai AICc        Nilai BIC
## 1   ARIMA(0,1,1) 1202.28737789285 1202.56644766029 1209.78680690384
## 2   ARIMA(2,1,0) 1210.07442886563 1210.54501710093 1220.07366754695
## 3   ARIMA(2,1,1) 1199.84904150017 1200.56332721446 1212.34808985183
## 4   ARIMA(1,1,1) 1200.33390441973 1200.80449265503 1210.33314310105
## 5   ARIMA(3,1,2) 1203.82740849962 1205.19326215815 1221.32607619193

Model terbaik yang didapatkan adalah:

  • \(ARIMA(2,1,1)\) \(1197.85\) dengan fungsi arima

  • \(ARIMA(2,1,1)\) \(1199.85\) dengan fungsi Arima

  • \(ARIMA(0,1,1)\) \(1202.29\) dengan fungsi auto.arima

Berdasarkan model terbaik dari fungsi auto.arima yaitu \(ARIMA(0,1,1)\) nilai AIC sebesar \(1202.29\) dan model \(ARIMA(2,1,1)\) dengan fungsi arima dan Arima nilai AICnya sebesar \(1197.85\) dan \(1199.85\) diperoleh bahwa nilai AIC terkecil yaitu model \(ARIMA(2,1,1)\) dari fungsi arima dan Arima.

Oleh karena itu, model yang disarankan untuk dianalisis lebih lanjut yaitu dengan menggunakan model \(ARIMA(2,1,1)\) yang mana ini adalah model terbaik (dengan nilai AIC terkecil) dari kandidat model berdasarkan ACF, PACF, dan EACF.

Model \(ARIMA(2,1,1)\), maka \(Y_{t}\) diperoleh dari penjabaran operator backshift sehingga untuk model \(ARIMA(2,1,1)\): \[\phi_{p} (1-B)^{d} Y_{t}=\mu+\theta_{q}(B)e_{t}\] \[\phi_{2}(B)(1-B)^{1} Y_{t}=\mu+\theta_{1}(B)e_{t}\] \[(1-\phi_{1}B-\phi_{2}B^{2})(1-B)Y_{t}=\mu+e_{t}+\theta_{1}e_{t-1}\] \[(1-B-\phi_{1}B+\phi_{1}B^{2}-\phi_{2}B^{2}+\phi_{2}B^{3})Y_{t}=\mu+e_{t}+\theta_{1}e_{t-1}\] \[(1-(1+\phi_{1})B-(\phi_{2}-\phi_{1})B^{2}+\phi_{2}B^{3})Y_{t}=\mu+e_{t}+\theta_{1}e_{t-1}\] \[Y_{t}-(1+\phi_{1})Y_{t-1}-(\phi_{2}-\phi_{1})Y_{t-2}+\phi_{2}Y_{t-3}=\mu+e_{t}+\theta_{1}e_{t-1}\] \[Y_{t}=\mu+(1+\phi_{1})Y_{t-1}+(\phi_{2}-\phi_{1})Y_{t-2}-\phi_{2}Y_{t-3}+e_{t}+\theta_{1}e_{t-1}\] Penduga parameter \(\widehat{\mu}=19.2583\), \(\widehat{\phi}_{1}=0.2805\), \(\widehat{\phi}_{2}=0.1646\), \(\widehat{\theta}_{1}=1.0000\), dengan sigma^2 estimated as \(31052\) adalah nilai dugaan bagi \(\sigma_{e}^2\). Model terbaik yang diperoleh yaitu model \(ARIMA(2,1,1)\).

\[Y_{t}=19.2583+(1+0.2805)Y_{t-1}+(0.1646-0.2805)Y_{t-2}-0.1646Y_{t-3}+e_{t}+1.0000e_{t-1}\] \[Y_{t}=19.2583+(1.2805)Y_{t-1}-0.1159Y_{t-2}-0.1646Y_{t-3}+e_{t}+e_{t-1}\]

Diagnostic model

ARIMA211diag <- stats::arima(yt.diff1, order = c(2,0,1), method = "ML")
checkresiduals(ARIMA211diag)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,1) with non-zero mean
## Q* = 2.3947, df = 6, p-value = 0.8801
## 
## Model df: 4.   Total lags used: 10

Berdasarkan pl0t di atas terlihat bahwa sisaan tidak mengikuti sebaran normal. Selanjutnya, ditinjau dari plot ACF dan PACF terlihat bahwa tidak ada lag yang signifikan. Hal tersebut menunjukkan bahwa tidak ada gejala autokorelasi pada sisaan. Selanjutnya, untuk memastikan kembali akan dilakukan uji asumsi secara formal:

sisaan <- arima211$residuals

Uji formal normalitas data

\(H_{0}\) : sisaan menyebar normal

\(H_{1}\) : sisaan tidak mengikuti sebaran normal

# Uji formal normalitas data
ks.test(sisaan,"pnorm")
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.49999, p-value = 6.661e-16
## alternative hypothesis: two-sided

Hasil pengujian, Normalitas Data: Hasil : \(p-value=6.661e-16 < \alpha=0.05\) yang berarti TOLAK \(H_{0}\). Sisaan tidak menyebar normal

Uji nilai tengah sisaan

\(H_{0}:\mu = 0\)

\(H_{1}:\mu \neq 0\)

# Uji nilai tengah sisaan
t.test(sisaan, mu = 0, alternative = "two.sided")
## 
##  One Sample t-test
## 
## data:  sisaan
## t = -0.25014, df = 89, p-value = 0.8031
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -41.77238  32.43096
## sample estimates:
## mean of x 
## -4.670714

Nilai Tengah Sisaan Hasil : \(p-value=0.8031 > \alpha=0.05\) yang berarti TAK TOLAK \(H_{0}\). Nilai tengah sisaan sama dengan \(0\)

Uji autokorelasi

\(H_{0}\) : tidak ada autokorelasi

\(H_{1}\) : terdapat autokorelasi

# Uji autokorelasi
Box.test(sisaan, lag = 23 ,type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 18.307, df = 23, p-value = 0.7406

Autokorelasi Hasil : \(p-value=0.7406 > \alpha=0.05\) yang berarti TAK TOLAK \(H_{0}\). Tidak terdapat gejala autokorelasi

Kesimpulan: Asumsi terpenuhi, kecuali sisaan tidak menyebar normal.

Overfitting

Overfitting model ARIMA akan diaplikasikan pada \(ARIMA(2,1,2)\) dan \(ARIMA(3,1,1)\)

# Overfitting model
ovefitarima212<- arima(train.Yt.ts, order = c(2,0,2), method = "ML")
ovefitarima311<- arima(train.Yt.ts, order = c(3,0,1), method = "ML")
modelakursii <- c("AIC ARIMA(2,0,2)", "AIC ARIMA(3,0,1)")
ModelARIMAaku <- c(ovefitarima212$aic, ovefitarima311$aic)
Akurasiall <- cbind(modelakursii, ModelARIMAaku)
colnames(Akurasiall) <- c("Akurasi", "Nilai")
Akurasiall <- as.data.frame(Akurasiall)
Akurasiall
##            Akurasi            Nilai
## 1 AIC ARIMA(2,0,2) 1227.51089074401
## 2 AIC ARIMA(3,0,1) 1227.58226890449

Berdasarkan hasil overfitting model \(ARIMA(2,1,2)\) dan \(ARIMA(3,1,1)\) nilai AICnya sebesar \(1227.5108\) dan \(1227.5822\) diperoleh bahwa nilai AIC terkecil yaitu model \(ARIMA(2,1,1)\) dari fungsi arima dan Arima.

Model \(ARIMA(2,1,1)\) ini selanjutnya digunakan untuk peramalan.

Peramalan

Selanjutnya yaitu melakukan peramalan hasil dari pemodelan dengan menggunakan data training, yaitu sebagai berikut:

Fitted model \(ARIMA(2,1,1)\) dengan data training

predictARIMA211 <- fitted(arima211)
fittedARIMA211 <- cbind(yt.diff1, predictARIMA211)

Visualisasi hasil fitted model \(ARIMA(2,1,1)\) dengan data training

ts.plot(yt.diff1, col="blue", ylab = "Wind speed", xlab = "Year")
title(main = "Fitted Time Series Plot of Yt", cex.sub = 0.8)
points(yt.diff1, pch = 20, col = "blue")
par(col="red")
lines(predictARIMA211)

forecasting <- predict(arima211, n.ahead = 6)
hasilforcasting <- as.data.frame(forecasting)
hasilforcasting
##        pred       se
## 1 -8.670237 177.1599
## 2 27.513614 217.1457
## 3 16.976623 217.2297
## 4 19.977212 218.3785
## 5 19.084383 218.4996
## 6 19.327870 218.5745

Visualisasi hasil peramalan model \(ARIMA(2,1,1)\) \(6\) waktu ke depan

plot(arima211, n.ahead=6, col="blue", ylab = "Yt", xlab = "Year")
title(main = "Fitted Time Series Plot of Yt", cex.sub = 0.8)
points(yt.diff1, pch = 20, col = "blue")
par(col="red")
lines(predictARIMA211)

Fitting \(ARIMA(2,1,1)\) dengan data testing

arima211test <- arima(test.Yt.ts, order = c(2,1,1), include.mean = TRUE, method = "ML")
arima211test
## 
## Call:
## arima(x = test.Yt.ts, order = c(2, 1, 1), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ar2     ma1
##       0.0949  -0.3062  -0.413
## s.e.  0.6435   0.3120   0.704
## 
## sigma^2 estimated as 46593:  log likelihood = -149.67,  aic = 305.33
predictARIMA211test <- fitted(arima211test)
fittedARIMA211test <- cbind(test.Yt.ts, predictARIMA211test)

Visualisasi hasil fitted model \(ARIMA(2,1,1)\) dengan data testing

ts.plot(test.Yt.ts, col="blue", ylab = "Yt", xlab = "Day")
title(main = "Fitted Time Series Plot of Yt", cex.sub = 0.8)
points(test.Yt.ts, pch = 20, col = "blue")
par(col="red")
lines(predictARIMA211test)

forecastingtest <- predict(arima211test, n.ahead = 6)
hasilforcastingtest <- as.data.frame(forecastingtest)

Akurasi hasil peramalan

Akurasi hasil peramalan model \(ARIMA(2,1,1)\). Akurasi hasil peramalan model \(ARIMA(2,1,1)\) dibandingin dengan data aktual adalah Data testing.

akurasi.arima211training <- accuracy(yt.diff1, predictARIMA211)
akurasi.arima211testing <- accuracy(test.Yt.ts, predictARIMA211test)

compmodelforecast <- cbind(akurasi.arima211training, akurasi.arima211testing)
Akurasi <- c("ME", "RMSE", "MAE", "MPE", "MAPE")
ModelARIMA211train <- c(4.670714, 176.217, 127.1697, -0.4472064, 421.1571)
ModelARIMA211test <- c(47.70344, 211.111, 155.7373, 2.151168, 8.428226)
Akurasiall <- rbind(Akurasi, ModelARIMA211train, ModelARIMA211test)
Akurasiall <- as.data.frame(Akurasiall)
Akurasiall
##                          V1      V2       V3         V4       V5
## Akurasi                  ME    RMSE      MAE        MPE     MAPE
## ModelARIMA211train 4.670714 176.217 127.1697 -0.4472064 421.1571
## ModelARIMA211test  47.70344 211.111 155.7373   2.151168 8.428226

Berdasarkan hasil fitting model dengan menggunakan metode ARIMA yang diperoleh model terbaik yaitu \(ARIMA(2,1,1)\). Pemilihan model terbaik ini didapatkan berdasarkan nilai akurasi AIC yang terkecil dari model tentatif dan diperoleh akurasi dari model berdasarkan nilai MAPE sebesar \(421.1571\) untuk data train dan MAPE sebesar \(8.428226\) untuk data testing.

Berdasarkan range nilai MAPE bahwa kemampuan model dalam meramal:

rangeaku <- c("<10%", "10-20%", "20-50%", ">50%")
arti <- c("Kemampuan peramalan sangat baik", "Kemampuan peramalan baik", "Kemampuan peramalan layak", "Kemampuan peramalan buruk")
Akurasiall <- cbind(rangeaku, arti)
colnames(Akurasiall) <- c("Range MAPE", "Arti")
Akurasiall <- as.data.frame(Akurasiall)
Akurasiall
##   Range MAPE                            Arti
## 1       <10% Kemampuan peramalan sangat baik
## 2     10-20%        Kemampuan peramalan baik
## 3     20-50%       Kemampuan peramalan layak
## 4       >50%       Kemampuan peramalan buruk

Sehingga berdasarkan range tersebut model \(ARIMA(2,1,1)\) dengan nilai MAPE sebesar \(421.1571\) untuk data train dan MAPE sebesar \(8.428226\) untuk data testing memiliki tingkat peramalan yang sangat baik untuk data testing dan buruk untuk data train.

Jawaban Gugus Data 2

Pendahuluan

Seasonal Autoregressive Integrated Moving Average (SARIMA) merupakan pengembangan dari model Autoregressive Integrated Moving Average (ARIMA) pada data deret waktu yang memiliki pola musiman.

Model ARIMA musiman menggabungkan faktor-faktor non-musiman (regular) dan musiman dalam model multiplikasi, dengan notasi sebagai berikut : \(ARIMA(p,d,q)\times(P,D,Q)_{s}\) dengan :

  • p = non-seasonal AR order,

  • d = non-seasonal differencing,

  • q = non-seasonal MA order,

  • P = seasonal AR order,

  • D = seasonal differencing,

  • Q = seasonal MA order,

  • s = time span of repeating seasonal pattern.

dimana tahapan identifikasi model SARIMA sama halnya seperti yang dilakukan pada model ARIMA regular atau model ARIMA non-seasonal, yaitu :

  • Plot time series

  • Identifikasi model

  • Pendugaan parameter model

  • Seleksi Model

  • Melakukan peramalan menggunakan model terbaik

Membaca Data

setwd("D:/S2/IPB/STA542 Analisis Deret Waktu/8/Res")
data2 <- read.xlsx("latihan Praktikum 8.xlsx", "data2", header=TRUE)
head(data2)
##   tahun bulan.     Xt
## 1  1979      1 17.083
## 2  1979      2 19.146
## 3  1979      3 20.341
## 4  1979      4 15.431
## 5  1979      5 22.361
## 6  1979      6 21.511
length(data2$Xt)
## [1] 150

Forming data time series

Di sini dilakukan perubahan data Xt dari data.frame menjadi data ts atau time series dengan fungsi ts

data.ts <- ts(data2, start=c(1979,1),frequency=12)

Visualisasi data time series

Data yang digunakan terdiri dari \(150\) observasi atau yang diukur dalam \(150\) bulan. Hasil dari visualisasi data tersaji pada plot berikut:

ts.plot(data.ts[,"Xt"], type="l", ylab="Xt", col="blue")
title(main = "Time Series Plot of Xt", cex.sub = 0.8)
points(data.ts, pch = 20, col = "blue")

Plot data asal memperlihatkan pola musiman dengan \(s = 12\), serta terdapat perilaku nonstasioner baik dalam rataan maupun ragam.

Visualisasi per Musim

library(forecast)
seasonplot(data.ts[,"Xt"],12,main="Seasonal Plot of Xt", ylab="times",year.labels = TRUE, 
col=rainbow(18))

Visualisasi deskriptif

monthplot(data.ts[,"Xt"],ylab="Xt", col="blue")

boxplot(data.ts[,"Xt"]~data.ts[,"bulan."],ylab="Xt",xlab="Month", col="blue")

Boxplot menunjukkan ukuran keragaman (panjang box) dan pemusatan (median) yang sama pada masing-masing bulan.

Identifikasi model

Berdasarkan plot time series di atas bahwa data kemungkinan stasioner terhadap ragam maupun rataan, maka akan dilakukan uji asumsi formal.

Uji kehomogenan ragam

Uji asumsi formal terhadap kehomogenan ragam yang digunakan yaitu Fligner-Killen test, dimana:

\(H_{0}\) : Ragam homogen

\(H_{1}\) : Ragam tidak homogen

library(car)
fligner.test(Xt ~ tahun, data=data.ts)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Xt by tahun
## Fligner-Killeen:med chi-squared = 6.0862, df = 12, p-value = 0.9117

Berdasarkan hasil uji Fligner-Killeen dengan menggunakan taraf signifikansi \(\alpha = 5% = 0.05\), nilai \(p-value\) yang diperoleh yaitu \(0.91168\), maka \(p-value=0.91168 > \alpha = 5% = 0.05\) sehingga tolak \(H_{0}\) atau dengan kata lain ragam data sudah stasioner.

Splitting data

Data yang akan digunakan untuk menentukan model terbaik berdasarkan kandidat model akan dilakukan splitting data terlebih dahulu yaitu dengan membagi data menjadi dua bagian bagian pertama yaitu sebagai data training sebanyak \(80\%\) atau \(120\) observasi dan bagian kedua sebagai data testing sebanyak \(20\%\) atau \(30\) observasi. Splitting data ini menggunakan fungsi subset. Data train inilah yang nantinya digunakan untuk membangun model.

train.ts <- subset(data.ts[,"Xt"],start=1,end=120)
test.ts <- subset(data.ts[,"Xt"],start=121,end=150)

Visualisasi data training yang akan digunakan untuk mencari model terbaik

ts.plot(train.ts, col="blue", ylab = "Xt", xlab = "Year")
title(main = "Time Series Train Plot of Xt", cex.sub = 0.8)
points(train.ts, pch = 20, col = "blue")

Visualisasi data testing yang akan digunakan untuk melihat akurasi dari model terbaik

ts.plot(test.ts, col="blue", ylab = "Xt", xlab = "Day")
title(main = "Time Series Testing Plot of Xt", cex.sub = 0.8)
points(test.ts, pch = 20, col = "blue")

ARIMA NON-SEASONAL

Identifikasi kestasioneran data

acf0 <- acf(train.ts,main="ACF",lag.max=48,xaxt="n", col="blue")
axis(1, at=0:48/12, labels=0:48)

acf0$lag <- acf0$lag * 12
acf0.1 <- as.data.frame(cbind(acf0$acf,acf0$lag))
acf0.2 <- acf0.1[which(acf0.1$V2%%12==0),]
barplot(height = acf0.2$V1, 
names.arg=acf0.2$V2, ylab="ACF", xlab="Lag")

Plot ACF menunjukkan data Xt tidak stasioner sehingga perlu dilakukan proses differencing.

Differencing data

d1 <- diff(train.ts)
ts.plot(d1, type="l", ylab="d1 Xt", col="blue")

Differencing non-seasonal \(d = 1\) jika dilihat berdasarkan plot di atas berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen non-seasonal.

Identifikasi kestasioneran data hasil differencing

acf1 <- acf(d1,lag.max=48,xaxt="n", main="ACF d1", col="blue")
axis(1, at=0:48/12, labels=0:48)

acf2 <- acf1$lag <- acf1$lag * 12
acf1.1 <- as.data.frame(cbind(acf1$acf,acf1$lag))
acf1.2 <- acf1.1[which(acf1.1$V2%%12==0),]
barplot(height = acf1.2$V1, names.arg=acf1.2$V2, ylab="ACF", xlab="Lag")

Plot ACF data non-seasonal differencing \(d = 1\) dan untuk mengkonfirmasi kestasioneran komponen non-seasonal (namun perhatikan lag 12,24, dst), pada series seasonal belum stasioner.

ARIMA SEASONAL

Differencing data seasonal

Melakukan proses differencing pada series seasonal

D1 <- diff(train.ts,12)
ts.plot(D1, type="l", ylab="D1 Xt", col="blue")

acf2<-acf(D1,lag.max=48,xaxt="n", main="ACF D1", col="blue")

acf2$lag <- acf2$lag * 12
acf2.1 <- as.data.frame(cbind(acf2$acf,acf2$lag))
acf2.2 <- acf2.1[which(acf2.1$V2%%12==0),]
barplot(height = acf2.2$V1, names.arg=acf2.2$V2, ylab="ACF", xlab="Lag")

Non-seasonal differencing D = 12 berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen seasonalnya (namun tidak untuk komponen n0n-seasonalnya).

Differencing

Melakukan differencing pertama pada data yang telah dilakukan differencing pada series seasonalya.

d1D1 <- diff(D1)
ts.plot(d1D1, type="l", ylab="d1 D1 Xt", col="blue")

Jika dilihat berdasarkan plot di atas maka data sudah stasioner terhadap rataan.

Identifikasi Model

acf3 <- acf(d1D1,lag.max=48,xaxt="n", main="ACF d1D1", col="blue")
axis(1, at=0:48/12, labels=0:48)

pacf3 <- pacf(d1D1,lag.max=48,xaxt="n", main="PACF d1D1", col="blue")
axis(1, at=0:48/12, labels=0:48)

Berdasarkan plot acf dan pacf di atas menujukkan pola cut-off di lag 2 dan pacfnya tails off di lag 2.

acf3$lag <- acf3$lag * 12
acf3.1 <- as.data.frame(cbind(acf3$acf,acf3$lag))
acf3.2 <- acf3.1[which(acf3.1$V2%%12==0),]
barplot(height = acf3.2$V1, names.arg=acf3.2$V2, ylab="ACF", 
xlab="Lag")

pacf3$lag <- pacf3$lag * 12
pacf3.1 <- as.data.frame(cbind(pacf3$acf,pacf3$lag))
pacf3.2 <- pacf3.1[which(pacf3.1$V2%%12==0),]
barplot(height = pacf3.2$V1, names.arg=pacf3.2$V2, ylab="PACF", xlab="Lag")

Untuk seasonalnya diperoleh \(ARIMA(0,1,1)_{12}\)

EACF

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

Karena, kedua komponen telah stasioner. Identifikasi komponen non-seasonal adalah \(ARIMA(0,1,2)\), \(ARIMA(2,1,2)\), \(ARIMA(2,1,0)\), \(ARIMA(1,1,1)\), dan \(ARIMA(0,1,1)\). Identifikasi komponen seasonal adalah \(ARIMA(0,1,1)_{12}\), sehingga model tentatif yang diperoleh adalah:

  • \(ARIMA(0,1,2)\times(0,1,1)12\)

  • \(ARIMA(2,1,2)\times(0,1,1)12\)

  • \(ARIMA(2,1,0)\times(0,1,1)12\)

  • \(ARIMA(1,1,1)\times(0,1,1)12\)

  • \(ARIMA(0,1,1)\times(0,1,1)12\)

tmodel1 <- Arima(train.ts,order=c(0,1,2),seasonal=c(0,1,1))
summary(tmodel1)
## Series: train.ts 
## ARIMA(0,1,2)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     ma2     sma1
##       -0.8570  0.2182  -0.4201
## s.e.   0.0939  0.0851   0.0963
## 
## sigma^2 estimated as 1.335:  log likelihood=-167.3
## AIC=342.6   AICc=342.99   BIC=353.29
## 
## Training set error measures:
##                      ME     RMSE       MAE        MPE     MAPE       MASE
## Training set 0.08552043 1.075466 0.7926163 0.09748527 1.197007 0.06069222
##                     ACF1
## Training set -0.02992964
#checkresiduals(tmodel1,lag=c(10,30))

tmodel2 <- Arima(train.ts,order=c(2,1,2),seasonal=c(0,1,1))
summary(tmodel2)
## Series: train.ts 
## ARIMA(2,1,2)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ar2     ma1      ma2     sma1
##       -1.1697  -0.1758  0.3640  -0.5464  -0.5058
## s.e.   0.1380   0.1363  0.1097   0.0929   0.0943
## 
## sigma^2 estimated as 1.26:  log likelihood=-163.85
## AIC=339.7   AICc=340.54   BIC=355.74
## 
## Training set error measures:
##                     ME     RMSE       MAE       MPE     MAPE       MASE
## Training set 0.1055238 1.034696 0.7605394 0.1129866 1.140696 0.05823603
##                     ACF1
## Training set -0.02840871
#checkresiduals(tmodel,lag=c(10,30))

tmodel3 <- Arima(train.ts,order=c(2,1,0),seasonal=c(0,1,1))
summary(tmodel3)
## Series: train.ts 
## ARIMA(2,1,0)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ar2     sma1
##       -0.7852  -0.2438  -0.3905
## s.e.   0.0943   0.0937   0.0983
## 
## sigma^2 estimated as 1.392:  log likelihood=-169.29
## AIC=346.59   AICc=346.98   BIC=357.28
## 
## Training set error measures:
##                      ME     RMSE       MAE        MPE     MAPE       MASE
## Training set 0.05500423 1.098231 0.8028162 0.05257322 1.222689 0.06147325
##                     ACF1
## Training set -0.06216435
#checkresiduals(tmodel3,lag=c(10,30))

tmodel4 <- Arima(train.ts,order=c(1,1,1),seasonal=c(0,1,1))
summary(tmodel4)
## Series: train.ts 
## ARIMA(1,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ma1     sma1
##       -0.3490  -0.5248  -0.4213
## s.e.   0.1233   0.1101   0.0956
## 
## sigma^2 estimated as 1.322:  log likelihood=-166.79
## AIC=341.58   AICc=341.97   BIC=352.27
## 
## Training set error measures:
##                      ME     RMSE       MAE        MPE     MAPE       MASE
## Training set 0.08785643 1.070214 0.7869391 0.09801171 1.190308 0.06025751
##                      ACF1
## Training set -0.004395688
#checkresiduals(tmodel4,lag=c(10,30))

tmodel5 <- Arima(train.ts,order=c(0,1,1),seasonal=c(0,1,1))
summary(tmodel5)
## Series: train.ts 
## ARIMA(0,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ma1     sma1
##       -0.6986  -0.3596
## s.e.   0.0570   0.0955
## 
## sigma^2 estimated as 1.408:  log likelihood=-170.29
## AIC=346.58   AICc=346.82   BIC=354.6
## 
## Training set error measures:
##                      ME     RMSE       MAE       MPE     MAPE      MASE
## Training set 0.09156746 1.109937 0.8307775 0.1083214 1.264212 0.0636143
##                    ACF1
## Training set -0.1859463
#checkresiduals(tmodel5,lag=c(10,30))

Model terbaik

AICKandidatModel <- c(tmodel1$aic, tmodel2$aic, tmodel3$aic, tmodel4$aic, tmodel5$aic)
AICcKandidatModel <- c(tmodel1$aicc, tmodel2$aicc, tmodel3$aicc, tmodel4$aicc, tmodel5$aicc)
BICKandidatModel <- c(tmodel1$bic, tmodel2$bic, tmodel3$bic, tmodel4$bic, tmodel5$bic)
KandidatModelARIMA <- c("ARIMA(0,1,2)(0,1,1)12", "ARIMA(2,1,2)(0,1,1)12", "ARIMA(2,1,0)(0,1,1)12", "ARIMA(1,1,1)(0,1,1)12", "ARIMA(2,1,1)(0,1,1)12")
compmodelARIMA <- cbind(KandidatModelARIMA, AICKandidatModel, AICcKandidatModel, BICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC", "Nilai AICc", "Nilai BIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA
##          Kandidat Model        Nilai AIC       Nilai AICc        Nilai BIC
## 1 ARIMA(0,1,2)(0,1,1)12 342.601426515043 342.993583377788  353.29274185289
## 2 ARIMA(2,1,2)(0,1,1)12 339.704960279606 340.544960279606 355.741933286378
## 3 ARIMA(2,1,0)(0,1,1)12 346.589242048885  346.98139891163 357.280557386732
## 4 ARIMA(1,1,1)(0,1,1)12 341.576304570268 341.968461433013 352.267619908116
## 5 ARIMA(2,1,1)(0,1,1)12 346.582020679753 346.815030388491 354.600507183139

Model terbaik diperoleh Berdasarkan nilai AIC dan AICc terkecil dari kandidat model. Oleh karena itu, Model terbaik yang diperoleh yaitu \(ARIMA(2,1,2)\times(0,1,1)12\)

Dibandingkan auto.arima

Pemilihan model terbaik juga akan dilakukan dengan menggunakan fungsi auto.arima

auto.arima(train.ts)
## Series: train.ts 
## ARIMA(1,1,1)(0,1,1)[12] 
## 
## Coefficients:
##           ar1      ma1     sma1
##       -0.3490  -0.5248  -0.4213
## s.e.   0.1233   0.1101   0.0956
## 
## sigma^2 estimated as 1.322:  log likelihood=-166.79
## AIC=341.58   AICc=341.97   BIC=352.27

Jika berdasarkan fungsi auto.arima model terbaik yang diperoleh yaitu \(ARIMA(1,1,1)\times(0,1,1)12\).

Model terbaik yang didapatkan adalah:

  • \(ARIMA(2,1,2)\times(0,1,1)12\) \(339.70\) dengan fungsi Arima

  • \(ARIMA(1,1,1)\times(0,1,1)12\) \(341.97\) dengan fungsi auto.arima

Berdasarkan model terbaik dari fungsi auto.arima yaitu \(ARIMA(1,1,1)\times(0,1,1)12\) nilai AIC sebesar \(341.97\) dan model \(ARIMA(2,1,2)\times(0,1,1)12\) dengan fungsi Arima nilai AICnya sebesar \(339.70\) diperoleh bahwa nilai AIC terkecil yaitu model \(ARIMA(2,1,2)\times(0,1,1)12\) \(339.70\). Oleh karena itu, model yang disarankan untuk dianalisis lebih lanjut yaitu dengan menggunakan model \(ARIMA(2,1,2)\times(0,1,1)12\) yang mana ini adalah model terbaik (dengan nilai AIC terkecil) dari kandidat model berdasarkan ACF, PACF, dan EACF.

Pengujian parameter model

Di sini digunakan User Defined Function printstatarima

printstatarima <- function (x, digits = 4,se=TRUE){
       if (length(x$coef) > 0) {
         cat("\nCoefficients:\n")
         coef <- round(x$coef, digits = digits)
         if (se && nrow(x$var.coef)) {
           ses <- rep(0, length(coef))
           ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits = digits)
           coef <- matrix(coef, 1, dimnames = list(NULL, names(coef)))
           coef <- rbind(coef, s.e. = ses)
           statt <- coef[1,]/ses
           pval  <- 2*pt(abs(statt), df=length(x$residuals)-1, lower.tail = FALSE)
           coef <- rbind(coef, t=round(statt,digits=digits),sign.=round(pval,digits=digits))
           coef <- t(coef)
         }
         print.default(coef, print.gap = 2)
       }
     }
printstatarima(tmodel2)
## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -1.1697  0.1380  -8.4761  0.0000
## ar2   -0.1758  0.1363  -1.2898  0.1996
## ma1    0.3640  0.1097   3.3181  0.0012
## ma2   -0.5464  0.0929  -5.8816  0.0000
## sma1  -0.5058  0.0943  -5.3637  0.0000

Model terbaik adalah model \(ARIMA(2,1,2)\times(0,1,1)12\) karena semua dugaan parameter model berpengaruh nyata, keculai \(\widehat{\phi_{2}}\).

Penduga parameter \(\widehat{\phi}_{1}=-1.1697\), \(\widehat{\phi}_{2}=-0.1758\), \(\widehat{\theta}_{1}=-0.3640\), \(\widehat{\theta}_{2}=0.5464\), \(\widehat{\Theta}_{1}=0.5058\) dan dengan sigma^2 estimated as \(1.26\) adalah nilai dugaan bagi \(\sigma_{e}^2\). Model terbaik yang diperoleh yaitu model \(ARIMA(2,1,2)\times(0,1,1)12\).

Model \(ARIMA(2,1,2)\times(0,1,1)12\), maka Xt diperoleh dari penjabaran operator backshift sehingga untuk model \(ARIMA(2,1,2)\times(0,1,1)12\):

  • \(p=2\), \(d=1\), \(q=2\)

  • \(P=0\), \(D=1\), \(Q=1\), \(s=12\)

\[\phi_{p}(B)\Phi_{P}(B^{s})(1-B)^{d}(1-B^{s})^{D}X_{t}=\theta_{q}(B)\Theta_{Q}(B^{s})e_{t}\]

\[\phi_{2}(B)\Phi_{0}(B^{12})(1-B)^{1}(1-B^{12})^{1}X_{t}=\theta_{2}(B)\Theta_{1}(B^{12})e_{t}\]

\[\phi_{2}(B)(1-B)(1-B^{12})X_{t}=\theta_{2}(B)\Theta_{1}(B^{12})e_{t}\]

\[(1-\phi_{1}B-\phi_{2}B^{2})(1-B)(1-B^{12})X_{t}=(1+\theta_{1}B+\theta_{2}B^{2})(1+\Theta_{1}B^{12})e_{t}\]

\((1-B+\phi_{1}B+\phi_{1}B^{2}-\phi_{2}B^{2}+\phi_{2}B^{3}-B^{12}+B^{13}+\phi_{1}B^{13}-\phi_{1}B^{14}+\phi_{2}B^{14}-\phi_{2}B^{15})X_{t}=\) \((1+\theta_{1}B+\theta_{2}B^{2}+\Theta_{1}B^{12}+\theta_{1}\Theta_{1}B^{13}+\theta_{2}\Theta_{1}B^{14})e_{t}\)

\((1-(1+\phi_{1})B+(\phi_{1}-\phi_{2})B^{2}+\phi_{2}B^{3}-B^{12}+(1+\phi_{1})B^{13}-(\phi_{1}-\phi_{2})B^{14}-\phi_{2}B^{15})X_{t}=\) \((1+\theta_{1}B+\theta_{2}B^{2}+\Theta_{1}B^{12}+\theta_{1}\Theta_{1}B^{13}+\theta_{2}\Theta_{1}B^{14})e_{t}\)

\(X_{t}-(1+\phi_{1})X_{t-1}+(\phi_{1}-\phi_{2})X_{t-2}+\phi_{2}X_{t-3}-B^{12}+(1+\phi_{1})X_{t-13}-(\phi_{1}-\phi_{2})X_{t-14}-\phi_{2}X_{t-15}=\) \(e_{t}+\theta_{1}e_{t-1}+\theta_{2}e_{t-2}+\Theta_{1}e_{t-12}+\theta_{1}\Theta_{1}e_{t-13}+\theta_{2}\Theta_{1}e_{t-14}\)

\(X_{t}=\) \((1+\phi_{1})X_{t-1}-(\phi_{1}-\phi_{2})X_{t-2}-\phi_{2}X_{t-3}+X_{t-12}-(1+\phi_{1})X_{t-13}+(\phi_{1}-\phi_{2})X_{t-14}+\phi_{2}X_{t-15}e_{t}\) \(+\theta_{1}e_{t-1}+\theta_{2}e_{t-2}+\Theta_{1}e_{t-12}+\theta_{1}\Theta_{1}e_{t-13}+\theta_{2}\Theta_{1}e_{t-14}\)

\(X_{t}=\) \((1+(-1.1697))X_{t-1}-((-1.1697)-(-0.1758))X_{t-2}-(-0.1758)X_{t-3}\)+X_{t-12}$ \(-(1+(-1.1697))X_{t-13}+((-1.1697)-(-0.1758))X_{t-14}+(-0.1758)X_{t-15}e_{t}\) \(+(-0.3640)e_{t-1}+0.5464e_{t-2}+0.5058e_{t-12}+(-0.3640)(0.5058)e_{t-13}+0.54640.5058e_{t-14}\)

\(X_{t}=\) \((-0.1697)X_{t-1}+0.9939X_{t-2}+0.1758X_{t-3}+X_{t-12}+0.1697X_{t-13}-0.9939X_{t-14}-0.1758X_{t-15}e_{t}\) \(-0.3640e_{t-1}+0.5464e_{t-2}+0.5058e_{t-12}-0.1841112e_{t-13}+0.2763691e_{t-14}\)

Diagnostic model

tsdisplay(residuals(tmodel2), lag.max=45, main='ARIMA(2,1,2)(0,1,1)12 Model Residuals', col="blue")

Berdasarkan pl0t di atas terlihat bahwa sisaan tidak mengikuti sebaran normal. Selanjutnya, ditinjau dari plot ACF dan PACF terlihat bahwa ada lag yang signifikan. Hal tersebut menunjukkan bahwa kemungkinan ada gejala autokorelasi pada sisaan. Selanjutnya, untuk memastikan kembali akan dilakukan uji asumsi secara formal:

library(portes)
ljbtest2 <- LjungBox(residuals(tmodel2),lags=seq(5,30,5),order=0,season=1,squared.residuals=FALSE)
ljbtest2
##  lags statistic df   p-value
##     5  1.020759  5 0.9608741
##    10 14.483269 10 0.1520671
##    15 21.787550 15 0.1134968
##    20 25.619672 20 0.1787435
##    25 27.505092 25 0.3311334
##    30 30.897944 30 0.4204368

Berdasarkan hasil uji LjungBox di atas terdapat autokorelasi pada sisaan, karena nilai \(p-value\) ada lag yang tidak signifikan atau \(p-value > \alpha=5%=0.05\).

Selanjutnya dilakukan uji asumsi formal terhadap kenormalan sisan dengan menggunakan uji Jarque Bera.

library(tseries)
jarque.bera.test(residuals(tmodel2))
## 
##  Jarque Bera Test
## 
## data:  residuals(tmodel2)
## X-squared = 2.3179, df = 2, p-value = 0.3138

Berdasarkan hasil uji kenormalan dengan uji Jarque Bera sisaan menyebar normal, karena nilai \(p-value=0.3138 > \alpha=5%=0.05\).

Sehingga, berdasarkan uji asumsi secara formal model \(ARIMA(2,1,2)\times(0,1,1)12\) akan digunakan analisis lebih lanjut / forecasting.

Forecasting

Validasi model

ramalan_arima2 = forecast(tmodel2, 6)
ramalan_arima2
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 1989       144.9433 143.5050 146.3816 142.7436 147.1429
## Feb 1989       152.0529 150.5877 153.5180 149.8121 154.2936
## Mar 1989       149.7392 148.1574 151.3209 147.3200 152.1583
## Apr 1989       146.3791 144.7401 148.0181 143.8724 148.8858
## May 1989       149.3032 147.5683 151.0381 146.6499 151.9565
## Jun 1989       156.5131 154.7245 158.3017 153.7776 159.2485
accuracy(ramalan_arima2,test.ts)
##                      ME      RMSE       MAE        MPE      MAPE       MASE
## Training set  0.1055238 1.0346962 0.7605394  0.1129866 1.1406959 0.05823603
## Test set     -0.4924512 0.9040121 0.7366223 -0.3313611 0.4983915 0.05640465
##                     ACF1 Theil's U
## Training set -0.02840871        NA
## Test set      0.12013822 0.1879534
plot(ramalan_arima2, col="blue")

forecast_arima2 <- cbind(ramalan_arima2$mean,ramalan_arima2$lower,ramalan_arima2$upper)
forecast_arima2
##          ramalan_arima2$mean ramalan_arima2$lower.80% ramalan_arima2$lower.95%
## Jan 1989            144.9433                 143.5050                 142.7436
## Feb 1989            152.0529                 150.5877                 149.8121
## Mar 1989            149.7392                 148.1574                 147.3200
## Apr 1989            146.3791                 144.7401                 143.8724
## May 1989            149.3032                 147.5683                 146.6499
## Jun 1989            156.5131                 154.7245                 153.7776
##          ramalan_arima2$upper.80% ramalan_arima2$upper.95%
## Jan 1989                 146.3816                 147.1429
## Feb 1989                 153.5180                 154.2936
## Mar 1989                 151.3209                 152.1583
## Apr 1989                 148.0181                 148.8858
## May 1989                 151.0381                 151.9565
## Jun 1989                 158.3017                 159.2485

Transformasi balik

forecast_arima2 <- cbind(data2[145:150,1:2],forecast_arima2)
     forecast_arima2 <- as.data.frame(forecast_arima2)

Plot ramalan perbulan

plot(forecast_arima2[,2],forecast_arima2[,1], xlab='Month', ylab='Xt)', type='l', col='blue', lwd=1, ylim=c(140, 160))
     lines(forecast_arima2[,2],forecast_arima2[,3], col='red', lwd=1)
     lines(forecast_arima2[,2],forecast_arima2[,5], col='blue', lwd=1)
     lines(forecast_arima2[,2],forecast_arima2[,7], col='green',lwd=1)
     legend("topleft", c("Forecast","Lower 95%","Upper 95%"), 
     col=c("red","blue","green"), lty=1,lwd=1,cex=0.6, inset=0.02, box.lty=0)

Jawaban Gugus Data 3

Membaca data

setwd("D:/S2/IPB/STA542 Analisis Deret Waktu/8/Res")
data3 <- read.xlsx("latihan Praktikum 8.xlsx", "data3", header=TRUE)
head(data3)
##   tahun bulan.  Zt
## 1  1984      1 308
## 2  1984      2 317
## 3  1984      3 316
## 4  1984      4 259
## 5  1984      5 277
## 6  1984      6 308
length(data3$Zt)
## [1] 300

Forming data time series

Di sini dilakukan perubahan data Zt dari data.frame menjadi data ts atau time series dengan fungsi ts

data.ts <- ts(data3, start=c(1984,1),frequency=12)

Visualisasi data time series

Data yang digunakan terdiri dari \(300\) observasi atau yang diukur dalam \(300\) bulan. Hasil dari visualisasi data tersaji pada plot berikut:

ts.plot(data.ts[,"Zt"], type="l", ylab="Zt", col="blue")
title(main = "Time Series Plot of Zt", cex.sub = 0.8)
points(data.ts, pch = 20, col = "blue")

Plot data asal memperlihatkan pola musiman dengan \(s = 12\), serta terdapat perilaku nonstasioner baik dalam rataan maupun ragam.

Visualisasi per Musim

library(forecast)
seasonplot(data.ts[,"Zt"],12,main="Seasonal Plot of Zt", ylab="times",year.labels = TRUE, 
col=rainbow(18))

Visualisasi deskriptif

monthplot(data.ts[,"Zt"],ylab="Zt", col="blue")

boxplot(data.ts[,"Zt"]~data.ts[,"bulan."],ylab="Zt",xlab="Month", col="blue")

Boxplot menunjukkan ukuran keragaman (panjang box) dan pemusatan (median) yang berbeda pada masing-masing bulan.

Identifikasi model

Berdasarkan plot time series di atas bahwa data kemungkinan tidak stasioner terhadap ragam maupun rataan, maka akan dilakukan uji asumsi formal.

Uji kehomogenan ragam

Uji asumsi formal terhadap kehomogenan ragam yang digunakan yaitu Fligner-Killen test, dimana:

\(H_{0}\) : Ragam homogen

\(H_{1}\) : Ragam tidak homogen

library(car)
fligner.test(Zt ~ tahun, data=data.ts)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  Zt by tahun
## Fligner-Killeen:med chi-squared = 43.369, df = 24, p-value = 0.009025

Berdasarkan hasil uji Fligner-Killeen dengan menggunakan taraf signifikansi \(\alpha = 5% = 0.05\), nilai \(p-value\) yang diperoleh yaitu \(0.009025\), maka \(p-value=0.009025\) \(<\) \(\alpha = 5% = 0.05\) sehingga tak tolak \(H_{0}\) atau dengan kata lain ragam data tidak stasioner.

Sehingga perlu dilakukan proses transformasi data

Transformasi data

logdat3 <- log10(data.ts[,"Zt"])
ts.plot(logdat3, type="l", ylab="Log Zt", main = "Monthly logged Zt Data", xlab = "Year", col="blue")

Berdasarkan plot di atas transformasi logaritma tidak berhasil mengatasi ketidakstasioneran dalam ragam dan ketidakstasioneran dalam rataan masih nampak.

Deskriptif Data

monthplot(logdat3,ylab="Zt", col="blue")

boxplot(logdat3 ~ cycle(logdat3), col="blue", xlab = "Month", ylab = "log Zt", main = "Monthly Zt Data - Boxplot")

Boxplot menunjukkan ukuran keragaman (panjang box) hampir sama akan tetapi ukuran pemusatan (median) yang berbeda pada masing-masing bulan

Uji Kehomogenan Ragam

\(H_{0}\) : Ragam homogen

\(H_{1}\) : Ragam tidak homogen

fligner.test(logdat3~ tahun, data=data.ts)
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  logdat3 by tahun
## Fligner-Killeen:med chi-squared = 68.761, df = 24, p-value = 3.363e-06

Hasil uji Fligner-Killeen menunjukkan bahwa ragam semakin tidak stasioner setelah dilakukan transformasi log, karena nilai kurang dari \(5\%\) sehinga tolak \(H_{0}\) atau ragam tidak homogen.

Berdasarkan hasil Uji Fligner-Killeen di atas transformasi logaritma tidak berhasil mengatasi ketidakstasioneran dalam ragam sehingga data yang digunakan tetap menggunakan data Zt.

Splitting data

Data yang akan digunakan untuk menentukan model terbaik berdasarkan kandidat model akan dilakukan splitting data terlebih dahulu yaitu dengan membagi data menjadi dua bagian bagian pertama yaitu sebagai data training sebanyak \(80\%\) atau \(240\) observasi dan bagian kedua sebagai data testing sebesar \(20\%\) sebanyak \(60\) observasi. Splitting data ini menggunakan fungsi subset. Data train inilah yang nantinya digunakan untuk membangun model.

train.ts <- subset(data.ts[,"Zt"],start=1,end=240)
test.ts <- subset(data.ts[,"Zt"],start=241,end=300)

Visualisasi data training yang akan digunakan untuk mencari model terbaik

ts.plot(train.ts, col="blue", ylab = "Zt", xlab = "Year")
title(main = "Time Series Train Plot of Zt", cex.sub = 0.8)
points(train.ts, pch = 20, col = "blue")

Visualisasi data testing yang akan digunakan untuk melihat akurasi dari model terbaik

ts.plot(test.ts, col="blue", ylab = "Zt", xlab = "Day")
title(main = "Time Series Testing Plot of Zt", cex.sub = 0.8)
points(test.ts, pch = 20, col = "blue")

Identifikasi kestasioneran data

acf0 <- acf(train.ts,main="ACF",lag.max=48,xaxt="n", col="blue")
axis(1, at=0:48/12, labels=0:48)

acf0$lag <- acf0$lag * 12
acf0.1 <- as.data.frame(cbind(acf0$acf,acf0$lag))
acf0.2 <- acf0.1[which(acf0.1$V2%%12==0),]
barplot(height = acf0.2$V1, 
names.arg=acf0.2$V2, ylab="ACF", xlab="Lag")

Plot ACF menunjukkan data Zt tidak stasioner sehingga perlu dilakukan proses differencing.

Differencing data

d1 <- diff(train.ts)
ts.plot(d1, type="l", ylab="d1 Zt", col="blue")

Differencing non-seasonal \(d = 1\) jika dilihat berdasarkan plot di atas berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen non-seasonal.

Series: train.ts ARIMA(0,1,1)(1,0,2)[12] with drift

Identifikasi kestasioneran data hasil differencing

acf1 <- acf(d1,lag.max=48,xaxt="n", main="ACF d1", col="blue")
axis(1, at=0:48/12, labels=0:48)

acf2 <- acf1$lag <- acf1$lag * 12
acf1.1 <- as.data.frame(cbind(acf1$acf,acf1$lag))
acf1.2 <- acf1.1[which(acf1.1$V2%%12==0),]
barplot(height = acf1.2$V1, names.arg=acf1.2$V2, ylab="ACF", xlab="Lag")

Plot ACF data non-seasonal differencing \(d = 1\) \(cut-off\) di lag 1, mengkonfirmasi kestasioneran komponen non musiman (namun perhatikan lag 12,24, dst), pada series musiman sudah stasioner.

Identifikasi model

acf2 <- acf(d1,lag.max=48,xaxt="n", main="ACF d1D0", col = "blue")
axis(1, at=0:48/12, labels=0:48)

pacf2 <- pacf(d1,lag.max=48,xaxt="n", main="PACF d1D0", col = "blue")
axis(1, at=0:48/12, labels=0:48)

Berdasarkan plot acf dan pacf di atas menujukkan pola tails-off di lag 1 dan pacfnya tails off di lag 3.

acf2$lag <- acf2$lag * 12
acf2.1 <- as.data.frame(cbind(acf2$acf,acf2$lag))
acf2.2 <- acf2.1[which(acf2.1$V2%%12==0),]
barplot(height = acf2.2$V1, names.arg=acf2.2$V2, ylab="ACF", xlab="Lag")

pacf2$lag <- pacf2$lag * 12
pacf2.1 <- as.data.frame(cbind(pacf2$acf,pacf2$lag))
pacf2.2 <- pacf2.1[which(pacf2.1$V2%%12==0),]
barplot(height = pacf2.2$V1, names.arg=pacf2.2$V2, ylab="PACF", xlab="Lag")

dari plot di atas diperoleh nilai \(P=1\) dan \(Q=0\), sehingga \(ARIMA(P,D,Q)_{s}=ARIMA(1,0,0)_{12}\)

EACF

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

Karena, kedua komponen telah stasioner. Identifikasi komponen non-seasonal adalah \(ARIMA(0,1,1)\), \(ARIMA(1,1,2)\), \(ARIMA(2,1,2)\), dan \(ARIMA(3,1,3)\). Identifikasi komponen seasonal adalah \(ARIMA(1,0,0)_{12}\), sehingga model tentatif yang diperoleh adalah:

  • \(ARIMA(0,1,1)\times (1,0,0)_{12}\)

  • \(ARIMA(1,1,2)\times (1,0,0)_{12}\)

  • \(ARIMA(2,1,2)\times (1,0,0)_{12}\)

  • \(ARIMA(3,1,3)\times (1,0,0)_{12}\)

tmodel1 <- Arima(train.ts,order=c(0,1,1),seasonal=c(1,0,0))
summary(tmodel1)
## Series: train.ts 
## ARIMA(0,1,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1     sar1
##       -0.5434  -0.0029
## s.e.   0.0586   0.0667
## 
## sigma^2 estimated as 1364:  log likelihood=-1200.91
## AIC=2407.81   AICc=2407.91   BIC=2418.24
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 3.149686 36.70631 28.15673 0.08575194 6.492074 0.5398364
##                     ACF1
## Training set 0.005287994
#checkresiduals(tmodel1,lag=c(10,30))

tmodel2 <- Arima(train.ts,order=c(1,1,2),seasonal=c(1,0,0))
summary(tmodel2)
## Series: train.ts 
## ARIMA(1,1,2)(1,0,0)[12] 
## 
## Coefficients:
##           ar1     ma1      ma2    sar1
##       -0.8182  0.2587  -0.3910  0.0026
## s.e.   0.1278  0.1434   0.0976  0.0671
## 
## sigma^2 estimated as 1367:  log likelihood=-1200.16
## AIC=2410.31   AICc=2410.57   BIC=2427.7
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 2.988257 36.59124 28.09873 0.06044021 6.508181 0.5387245
##                    ACF1
## Training set 0.03216061
#checkresiduals(tmodel,lag=c(10,30))

tmodel3 <- Arima(train.ts,order=c(2,1,2),seasonal=c(1,0,0))
summary(tmodel3)
## Series: train.ts 
## ARIMA(2,1,2)(1,0,0)[12] 
## 
## Coefficients:
##           ar1     ar2     ma1      ma2    sar1
##       -0.5364  0.2019  0.0137  -0.4562  0.0073
## s.e.   0.2226  0.1127  0.2139   0.1182  0.0671
## 
## sigma^2 estimated as 1357:  log likelihood=-1198.79
## AIC=2409.58   AICc=2409.94   BIC=2430.44
## 
## Training set error measures:
##                   ME     RMSE      MAE       MPE     MAPE      MASE        ACF1
## Training set 3.44405 36.37927 27.88885 0.1520029 6.405166 0.5347005 -0.00925425
#checkresiduals(tmodel3,lag=c(10,30))

tmodel4 <- Arima(train.ts,order=c(3,1,3),seasonal=c(1,0,0))
summary(tmodel4)
## Series: train.ts 
## ARIMA(3,1,3)(1,0,0)[12] 
## 
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): NaNs produced
##           ar1     ar2      ar3     ma1      ma2      ma3    sar1
##       -0.6577  0.0876  -0.0029  0.1374  -0.4131  -0.0617  0.0088
## s.e.      NaN     NaN      NaN     NaN   0.1840      NaN  0.0671
## 
## sigma^2 estimated as 1369:  log likelihood=-1198.75
## AIC=2413.49   AICc=2414.12   BIC=2441.3
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
## Training set 3.407212 36.37233 27.92903 0.145034 6.414213 0.5354709 -0.01139828
#checkresiduals(tmodel4,lag=c(10,30))

Model Terbaik

AICKandidatModel <- c(tmodel1$aic, tmodel2$aic, tmodel3$aic, tmodel4$aic)
AICcKandidatModel <- c(tmodel1$aicc, tmodel2$aicc, tmodel3$aicc, tmodel4$aicc)
BICKandidatModel <- c(tmodel1$bic, tmodel2$bic, tmodel3$bic, tmodel4$bic)
KandidatModelARIMA <- c("ARIMA(0,1,1)(1,0,0)12", "ARIMA(1,1,2)(1,0,0)12", "ARIMA(2,1,2)(1,0,0)12", "ARIMA(3,1,3)(1,0,0)12")
compmodelARIMA <- cbind(KandidatModelARIMA, AICKandidatModel, AICcKandidatModel, BICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC", "Nilai AICc", "Nilai BIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA
##          Kandidat Model        Nilai AIC       Nilai AICc        Nilai BIC
## 1 ARIMA(0,1,1)(1,0,0)12 2407.81003670379 2407.91216436337 2418.23942735959
## 2 ARIMA(1,1,2)(1,0,0)12 2410.31447495331 2410.57198568293 2427.69679271297
## 3 ARIMA(2,1,2)(1,0,0)12 2409.57795294597 2409.94002191149 2430.43673425756
## 4 ARIMA(3,1,3)(1,0,0)12 2413.49296573106 2414.11905268758 2441.30467414651

Model terbaik diperoleh Berdasarkan nilai AIC dan AICc terkecil dari kandidat model. Oleh karena itu, Model terbaik yang diperoleh yaitu \(ARIMA(0,1,1)\times(1,0,0)12\)

Dibandingkan auto.arima

Pemilihan model terbaik juga akan dilakukan dengan menggunakan fungsi auto.arima

auto.arima(data.ts[,"Zt"])
## Series: data.ts[, "Zt"] 
## ARIMA(0,1,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1     sar1
##       -0.5081  -0.0208
## s.e.   0.0517   0.0602
## 
## sigma^2 estimated as 1402:  log likelihood=-1506.63
## AIC=3019.27   AICc=3019.35   BIC=3030.37

Jika berdasarkan fungsi auto.arima model terbaik yang diperoleh yaitu \(ARIMA(0,1,1)\times(1,0,0)12\).

Dipilih model terbaik berdasarkan nilai AIC yang terkecil yaitu model \(ARIMA(0,1,1)\times(1,0,0)12\) yang akan dilanjutkan analisisnya.

Pengujian parameter model

Di sini digunakan User Defined Function printstatarima

printstatarima <- function (x, digits = 4,se=TRUE){
       if (length(x$coef) > 0) {
         cat("\nCoefficients:\n")
         coef <- round(x$coef, digits = digits)
         if (se && nrow(x$var.coef)) {
           ses <- rep(0, length(coef))
           ses[x$mask] <- round(sqrt(diag(x$var.coef)), digits = digits)
           coef <- matrix(coef, 1, dimnames = list(NULL, names(coef)))
           coef <- rbind(coef, s.e. = ses)
           statt <- coef[1,]/ses
           pval  <- 2*pt(abs(statt), df=length(x$residuals)-1, lower.tail = FALSE)
           coef <- rbind(coef, t=round(statt,digits=digits),sign.=round(pval,digits=digits))
           coef <- t(coef)
         }
         print.default(coef, print.gap = 2)
       }
     }
printstatarima(tmodel1)
## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -0.5434  0.0586  -9.2730  0.0000
## sar1  -0.0029  0.0667  -0.0435  0.9654

Model terbaik \(ARIMA(0,1,1)\times(1,0,0)12\) dimana dugaan parameter model \(ARIMA(0,1,1)\times(1,0,0)12\) \(\Phi_{1}\) tidak berpengaruh nyata.

Model \(ARIMA(0,1,1)\times(1,0,0)12\), maka \(Z_{t}\) diperoleh dari penjabaran operator backshift sehingga untuk model \(ARIMA(0,1,1)\times(1,0,0)12\):

  • \(p=0\), \(d=1\), \(q=1\)

  • \(P=1\), \(D=0\), \(Q=0\), \(s=12\)

\[\phi_{p}(B)\Phi_{P}(B^{s})(1-B)^{d}(1-B^{s})^{D}Z_{t}=\theta_{q}(B)\Theta_{Q}(B^{s})e_{t}\] \[\phi_{0}(B)\Phi_{1}(B^{12})(1-B)^{1}(1-B^{s})^{0}Z_{t}=\theta_{1}(B)\Theta_{0}(B^{12})e_{t}\] \[\Phi_{1}(B^{12})(1-B)Z_{t}=\theta_{1}(B)e_{t}\] \[(1-\Phi_{1}B^{12})(1-B)Z_{t}=(1+\theta_{1}B)e_{t}\] \[(1-B-\Phi_{1}B^{12}+\Phi_{1}B^{13})Z_{t}=(1+\theta_{1}B)e_{t}\] \[Z_{t}-Z_{t-1}-\Phi_{1}Z_{t-12}+\Phi_{1}Z_{t-13}=e_{t}+\theta_{1}e_{t-1}\] \[Z_{t}=Z_{t-1}+\Phi_{1}Z_{t-12}-\Phi_{1}Z_{t-13}e_{t}+\theta_{1}e_{t-1}\] Penduga parameter \(\widehat{\theta}_{1}=0.5434\), \(\widehat{\Phi}_{1}=-0.0029\), dengan sigma^2 estimated as \(1364\) adalah nilai dugaan bagi \(\sigma_{e}^2\). Model terbaik yang diperoleh yaitu model \(ARIMA(0,1,1)\times(1,0,0)12\). \[Z_{t}=Z_{t-1}-0.0029Z_{t-12}+0.0029Z_{t-13}e_{t}+0.5434e_{t-1}\]

Diagnostic model

tsdisplay(residuals(tmodel1), col="blue", lag.max=45, main='ARIMA(0,1,1)(1,0,0)12 Model Residuals')

Berdasarkan pl0t di atas terlihat bahwa plot ACF dan PACF terlihat bahwa tidak ada lag yang signifikan. Hal tersebut menunjukkan bahwa tidak ada gejala autokorelasi pada sisaan. Selanjutnya, untuk memastikan kembali akan dilakukan uji asumsi secara formal:

library(portes)
ljbtest1 <- LjungBox(residuals(tmodel1),lags=seq(5,30,5),order=0,season=1,squared.residuals=FALSE)
ljbtest1
##  lags statistic df    p-value
##     5  5.679413  5 0.33867328
##    10  9.145621 10 0.51834169
##    15 16.133395 15 0.37324198
##    20 26.753845 20 0.14233251
##    25 35.808522 25 0.07455586
##    30 47.104749 30 0.02430405

Berdasarkan hasil uji LjungBox di atas terdapat autokorelasi pada sisaan, karena nilai \(p-value\) semua lag tidak signifikan atau \(p-value > \alpha=5\%=0.05\).

Selanjutnya dilakukan uji asumsi formal terhadap kenormalan sisan dengan menggunakan uji Jarque Bera.

library(tseries)
jarque.bera.test(residuals(tmodel1))
## 
##  Jarque Bera Test
## 
## data:  residuals(tmodel1)
## X-squared = 1.7241, df = 2, p-value = 0.4223

Berdasarkan hasil uji kenormalan dengan uji Jarque Bera sisaan menyebar normal, karena nilai \(p-value=0.4223 > \alpha=5\%=0.05\).

Sehingga, berdasarkan uji asumsi secara formal model \(ARIMA(1,1,1)\times(2,0,1)12\) akan digunakan analisis lebih lanjut / forecasting.

Forecasting

Validasi model

ramalan_arima1 = forecast(tmodel1, 6)
ramalan_arima1
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2004       652.9313 605.5935 700.2692 580.5344 725.3283
## Feb 2004       652.8529 600.8138 704.8919 573.2660 732.4397
## Mar 2004       652.9023 596.5529 709.2516 566.7233 739.0812
## Apr 2004       653.0214 592.6688 713.3741 560.7200 745.3229
## May 2004       652.7279 588.6214 716.8343 554.6855 750.7702
## Jun 2004       652.9081 585.2558 720.5603 549.4429 756.3733
accuracy(ramalan_arima1,test.ts)
##                     ME     RMSE      MAE        MPE     MAPE      MASE
## Training set  3.149686 36.70631 28.15673 0.08575194 6.492074 0.5398364
## Test set     65.942693 68.20289 65.94269 9.11995328 9.119953 1.2642898
##                     ACF1 Theil's U
## Training set 0.005287994        NA
## Test set     0.073548472  3.358747
plot(ramalan_arima1, col="blue")

forecast_arima1 <- cbind(ramalan_arima1$mean,ramalan_arima1$lower,ramalan_arima1$upper)
forecast_arima1
##          ramalan_arima1$mean ramalan_arima1$lower.80% ramalan_arima1$lower.95%
## Jan 2004            652.9313                 605.5935                 580.5344
## Feb 2004            652.8529                 600.8138                 573.2660
## Mar 2004            652.9023                 596.5529                 566.7233
## Apr 2004            653.0214                 592.6688                 560.7200
## May 2004            652.7279                 588.6214                 554.6855
## Jun 2004            652.9081                 585.2558                 549.4429
##          ramalan_arima1$upper.80% ramalan_arima1$upper.95%
## Jan 2004                 700.2692                 725.3283
## Feb 2004                 704.8919                 732.4397
## Mar 2004                 709.2516                 739.0812
## Apr 2004                 713.3741                 745.3229
## May 2004                 716.8343                 750.7702
## Jun 2004                 720.5603                 756.3733

Plot ramalan perbulan

forecast_arima1 <- cbind(data3[295:300,1:2],forecast_arima1)
     forecast_arima1 <- as.data.frame(forecast_arima1)
plot(forecast_arima1[,2],forecast_arima1[,1], xlab='Month', ylab='Xt)', type='l', col='black', lwd=1, ylim=c(550, 770))
     lines(forecast_arima1[,2],forecast_arima1[,3], col='red', lwd=1)
     lines(forecast_arima1[,2],forecast_arima1[,5], col='blue', lwd=1)
     lines(forecast_arima1[,2],forecast_arima1[,7], col='green',lwd=1)
     legend("topleft", c("Forecast","Lower 95%","Upper 95%"), 
     col=c("red","blue","green"), lty=1,lwd=1,cex=0.6, inset=0.02, box.lty=0)