Jawaban ini tersedia di : rpubs
Soal 2
Melalui program R, bangkitkan data \(Y_{t}\), (n=200), berupa model \(ARIMA(1,2,2)\) dengan \(\mu=0.35\), \(\phi=0.75\), \(\theta_{1}=-0.65\), dan \(\theta_{2}=0.40\) serta \(e_{t} \sim Normal(0, 1)\). Gunakan 175 data terakhir dan lakukan proses berikut:
Membangkitkan Data
Membangkitkan data \(Y_{t}\), (n=200)
set.seed(101)
e <- rnorm(200, 0, 1)
n <- length(e)Model \(ARIMA(1,2,2)\), maka \(Y_{t}\) diperoleh dari penjabaran operator backshift sehingga untuk model \(ARIMA(1,2,2)\): \[\phi_{p} (1-B)^{d} Y_{t}=\mu+\theta_{q}(B)e_{t}\] \[(1-\phi_{1}(B))(1-B)^{2} Y_{t}=\mu+\theta_{2}(B)e_{t}\] \[(1-\phi_{1}(B))(1-2B+B^{2}) Y_{t}=\mu+(1-\theta_{1}(B)-\theta_{2}B^{2})e_{t}\] \[(1-2B+B^{2}-\phi_{1}B+2\phi_{1}B^{2}-\phi_{1}B^{3}) Y_{t}=\mu+e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}\] \[(1-(2+\phi_{1})B+(1+2\phi_{1})B^{2}-\phi_{1}B^{3}) Y_{t}=\mu+e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}\] \[Y_{t}-(2+\phi_{1})Y_{t-1}+(1+2\phi_{1})Y_{t-2}-\phi_{1}Y_{t-3} =\mu+e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}\] \[Y_{t} =\mu+(2+\phi_{1})Y_{t-1}-(1+2\phi_{1})Y_{t-2}+\phi_{1}Y_{t-3}+e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}\] Konstanta dalam bangkitkan data \(Y_{t}\), (n=200) dengan model \(ARIMA(1,2,2)\)
mu <- 0.35
phi <- 0.75
tetha1 <- -0.65
tetha2 <- 0.40yt <- c(1:n)
for(i in 4:n){yt[i] <-
mu + ((2+phi)*yt[i-1]) - ((1+2*phi)*yt[i-2]) + (phi*yt[i-3]) + (e[i]) - (tetha1*e[i-1]) - (tetha2*e[i-3])}Data yang digunakan adalah 175 data terakhir, maka:
yt.arima122 <- yt[-c(1:25)] #membuang 25 data pertamaVisualisasi hasil data bangkitan dengan mengambil 175 data terakhir
ts.plot(yt.arima122, col="blue", xlab = "Time Period")
title(main = "Plot time series dari data bangkitan", cex.sub = 0.8)
points(yt.arima122, pch = 20, col = "blue")
- Identifikasilah kestasioneran data melalui correlogram dan uji ADF, serta lakukan proses differencing jika data tidak stasioner.
Jawaban poin a.
Mengidentifikasi kestasioneran data pada latihan ini akan dilakukan dengan 2 cara yaitu dengan correlogram dan dengan uji ADF
- Melalui Correlogram
acf(yt.arima122, lag.max = 20, col = "blue")
acf(yt.arima122, lag.max = 20, plot = FALSE)##
## Autocorrelations of series 'yt.arima122', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0.983 0.966 0.949 0.932 0.915 0.898 0.881 0.864 0.847 0.830 0.813 0.796 0.780
## 14 15 16 17 18 19 20
## 0.763 0.746 0.729 0.713 0.696 0.679 0.663Berdasarkan plot Correlogram di atas data yang dibangkitkan merupakan data yang tidak stasioner. Hal ini diindikasikan berdasarkan pola yang ditunjukkan oleh plot Correlogram yaitu pola yang menurun secara perlahan.
- Melalui 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(yt.arima122, alternative = "stationary", k=trunc((length(yt.arima122)-1)^(1/3)))##
## Augmented Dickey-Fuller Test
##
## data: yt.arima122
## Dickey-Fuller = -1.0449, Lag order = 5, p-value = 0.9286
## alternative hypothesis: stationaryBerdasarkan hasil uji ADF dengan taraf signifikansi \(\alpha = 5% = 0.05\) di atas dapat diketahui bersama bahwa, \(p-value = 0.9286 > \alpha = 0.05\) sehingga keputusan yang diambil adalah Terima \(H_{0}\) artinya data yang dibangkitkan merupakan data yang tidak stasioner
Maka, untuk data yang dibangkitkan agar bisa dilakukan analisis lebih lanjut maka perlunya dilakukan proses differencing
- Differencing ordo 1
yt.arima122.diff1 <- diff(yt.arima122, difference=1)Visualisasi hasil data bangkitan dengan mengambil 175 data terakhir dan dilakukan differencing 1 kali
ts.plot(yt.arima122.diff1, col="blue", xlab = "Time Period")
title(main = "Plot time series data bangkitan dg differencing 1 kali", cex.sub = 0.8)
points(yt.arima122.diff1, pch = 20, col = "blue")
Mengidentifikasi kestasioneran data pada data yang telah dibangkitkan dan dilakukan differencing 1 kali ini akan dilakukan dengan 2 cara yaitu dengan correlogram dan dengan uji ADF
- Melalui Correlogram
acf(yt.arima122.diff1, lag.max = 20, col = "blue")
acf(yt.arima122.diff1, lag.max = 20, plot = FALSE)##
## Autocorrelations of series 'yt.arima122.diff1', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0.981 0.962 0.942 0.923 0.904 0.886 0.868 0.850 0.834 0.817 0.801 0.785 0.770
## 14 15 16 17 18 19 20
## 0.754 0.738 0.722 0.707 0.691 0.675 0.659Berdasarkan plot Correlogram di atas data yang telah dibangkitkan dan dilakukan differencing 1 kali merupakan data yang tidak stasioner. Hal ini diindikasikan berdasarkan pola yang ditunjukkan oleh plot Correlogram yaitu pola yang menurun secara perlahan.
- Melalui 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(yt.arima122.diff1, alternative = "stationary", k=trunc((length(yt.arima122)-1)^(1/3)))##
## Augmented Dickey-Fuller Test
##
## data: yt.arima122.diff1
## Dickey-Fuller = -1.6651, Lag order = 5, p-value = 0.7161
## alternative hypothesis: stationaryBerdasarkan hasil uji ADF dengan taraf signifikansi \(\alpha = 5% = 0.05\) di atas dapat diketahui bersama bahwa, \(p-value = 0.71613 > \alpha = 0.05\) sehingga keputusan yang diambil adalah Terima \(H_{0}\) artinya data yang telah dibangkitkan dan dilakukan differencing 1 kali merupakan data yang tidak stasioner.
Maka, untuk data perlu dilakukan proses differencing ordo 2
- Differencing ordo 2
yt.arima122.diff2 <- diff(yt.arima122.diff1, difference=1)Visualisasi hasil data bangkitan dengan mengambil 175 data terakhir dan dilakukan differencing 2 kali
ts.plot(yt.arima122.diff2, col="blue", xlab = "Time Period")
title(main = "Plot time series data bangkitan dg differencing 2 kali", cex.sub = 0.8)
points(yt.arima122.diff2, pch = 20, col = "blue")
Mengidentifikasi kestasioneran data pada data yang telah dibangkitkan dan dilakukan differencing 2 kali ini akan dilakukan dengan 2 cara yaitu dengan correlogram dan dengan uji ADF
- Melalui Correlogram
acf(yt.arima122.diff2, lag.max = 20, col = "blue")
acf(yt.arima122.diff2, lag.max = 20, plot = FALSE)##
## Autocorrelations of series 'yt.arima122.diff2', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.769 0.359 0.064 -0.018 -0.011 0.042 0.102 0.137 0.124 0.067 0.020
## 12 13 14 15 16 17 18 19 20
## -0.001 0.013 0.034 0.061 0.056 0.034 -0.024 -0.086 -0.124Berdasarkan plot Correlogram di atas data yang telah dibangkitkan dan dilakukan differencing 2 kali merupakan data yang kemungkinan stasioner. Hal ini diindikasikan berdasarkan pola yang ditunjukkan oleh plot Correlogram yaitu pola yang terdapat lag. Maka, untuk memastikan hal ini akan dilakukan uji ADF
- Melalui 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(yt.arima122.diff2, alternative = "stationary", k=trunc((length(yt.arima122)-1)^(1/3)))##
## Augmented Dickey-Fuller Test
##
## data: yt.arima122.diff2
## Dickey-Fuller = -3.692, Lag order = 5, p-value = 0.02674
## alternative hypothesis: stationaryBerdasarkan hasil uji ADF dengan taraf signifikansi \(\alpha = 5% = 0.05\) di atas dapat diketahui bersama bahwa, \(p-value = 0.026735 < \alpha = 0.05\) sehingga keputusan yang diambil adalah Tolak \(H_{0}\) artinya data yang telah dibangkitkan dan dilakukan differencing 2 kali merupakan data yang stasioner
Maka, data yang telah dibangkitkan dan dilakukan differencing 2 kali dapat dilanjutkan untuk analisis.
- Selanjutnya, berdasarkan ACF, PACF, dan EACF, identifikasilah kandidat model yang sesuai.
Jawaban poin b.
Identifikasi kandidat model diperoleh berdasarkan nilai \(p\) dan \(q\) dimana nilai \(d=2\) (yaitu telah dilakukan differencing sebanyak 2 kali)
acf(yt.arima122.diff2, lag.max = 20, col = "blue")
acf(yt.arima122.diff2, lag.max = 20, plot = FALSE)##
## Autocorrelations of series 'yt.arima122.diff2', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.769 0.359 0.064 -0.018 -0.011 0.042 0.102 0.137 0.124 0.067 0.020
## 12 13 14 15 16 17 18 19 20
## -0.001 0.013 0.034 0.061 0.056 0.034 -0.024 -0.086 -0.124Berdasarkan plot ACF terlihat bahwa plot cuts off setelah lag ke-2, model ARIMA yang terbentuk adalah \(ARIMA(0,2,2)\)
pacf(yt.arima122.diff2, lag.max = 20, col = "blue")
pacf(yt.arima122.diff2, lag.max = 20, plot = FALSE)##
## Partial autocorrelations of series 'yt.arima122.diff2', by lag
##
## 1 2 3 4 5 6 7 8 9 10 11
## 0.769 -0.570 0.254 0.054 -0.151 0.249 -0.048 -0.004 0.058 -0.105 0.111
## 12 13 14 15 16 17 18 19 20
## -0.063 0.058 0.003 0.021 -0.057 0.036 -0.137 0.005 -0.008Berdasarkan plot PACF di atas model yang bisa diidentifikasi adalah, \(ARIMA(3,2,0)\)
eacf(yt.arima122.diff2)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x o o o o o o o o o o o o
## 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 o x o x o o o o o o o o o o
## 4 x o x x x o o o o o o o o o
## 5 x o x x o o o o o o o o o o
## 6 o o x x o o o x o o o o o o
## 7 x o o o o x o x o o o o o o
Untuk plot EACF tanda silang mulai berhenti pada pertemuan angka \(0\) dengan \(2\) dan \(1\) dengan \(2\). Karena telah dilakukan differencing sebanyak dua kali maka \(d=2\), maka dapat ditentukan beberapa kandidat model yaitu: \(ARIMA(0,2,2), ARIMA(3,2,0), ARIMA(3,2,2), ARIMA(1,2,2)\)
- Berdasarkan kandidat model tersebut, tentukan model terbaik berdasarkan nilai AIC-nya
Jawaban poin c.
# ARIMA(0,2,2)
arima022 <- arima(yt.arima122.diff2, order = c(0,0,2), include.mean = TRUE, method = "ML")
arima022##
## Call:
## arima(x = yt.arima122.diff2, order = c(0, 0, 2), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ma1 ma2 intercept
## 1.2989 0.7834 1.2348
## s.e. 0.0518 0.0459 0.2306
##
## sigma^2 estimated as 0.9789: log likelihood = -244.96, aic = 495.92# ARIMA(3,2,0)
arima320 <- arima(yt.arima122.diff2, order = c(3,0,0), include.mean = TRUE, method = "ML")
arima320##
## Call:
## arima(x = yt.arima122.diff2, order = c(3, 0, 0), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 intercept
## 1.3641 -0.8929 0.2615 1.2291
## s.e. 0.0732 0.1071 0.0730 0.2865
##
## sigma^2 estimated as 1.031: log likelihood = -249.06, aic = 506.13# ARIMA(3,2,2)
arima322 <- arima(yt.arima122.diff2, order = c(3,0,2), include.mean = TRUE, method = "ML")
arima322##
## Call:
## arima(x = yt.arima122.diff2, order = c(3, 0, 2), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 intercept
## 0.2371 0.0962 -0.0944 1.1900 0.6655 1.2294
## s.e. 0.1322 0.1328 0.1023 0.1127 0.0840 0.2738
##
## sigma^2 estimated as 0.9329: log likelihood = -240.77, aic = 493.54# ARIMA(1,2,2)
arima122 <- arima(yt.arima122.diff2, order = c(1,0,2), include.mean = TRUE, method = "ML")
arima122##
## Call:
## arima(x = yt.arima122.diff2, order = c(1, 0, 2), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ma1 ma2 intercept
## 0.2908 1.1362 0.6688 1.230
## s.e. 0.0979 0.0834 0.0691 0.289
##
## sigma^2 estimated as 0.9377: log likelihood = -241.19, aic = 490.38AIC <- c(495.92, 506.13, 493.54, 490.38)
KMARIMA <- c("ARIMA(0,2,2)", "ARIMA(3,2,0)","ARIMA(3,2,2)", "ARIMA(1,2,2)")
compmodel <- cbind(KMARIMA, AIC)
colnames(compmodel) <- c("Kandidat Model", "Nilai AIC")
compmodel <- as.data.frame(compmodel)
compmodel## Kandidat Model Nilai AIC
## 1 ARIMA(0,2,2) 495.92
## 2 ARIMA(3,2,0) 506.13
## 3 ARIMA(3,2,2) 493.54
## 4 ARIMA(1,2,2) 490.38Model terbaik diperoleh Berdasarkan nilai AIC terkecil dari kandidat model. Oleh karena itu, Model terbaik yang diperoleh yaitu \(ARIMA(1,2,2)\)
- Bandingkan penduga parameter yang diperoleh untuk model terbaik pada poin (c) tersebut dengan nilai parameter yang sesungguhnya. Apa kesimpulan Anda?
Jawaban poin d.
mu <- 0.35
phi <- 0.75
tetha1 <- -0.65
tetha2 <- 0.40
arima122##
## Call:
## arima(x = yt.arima122.diff2, order = c(1, 0, 2), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ma1 ma2 intercept
## 0.2908 1.1362 0.6688 1.230
## s.e. 0.0979 0.0834 0.0691 0.289
##
## sigma^2 estimated as 0.9377: log likelihood = -241.19, aic = 490.38Penduga parameter \(\widehat{\mu}=1.23001\), \(\widehat{\phi}=0.29078\), \(\widehat{\theta}_{1}=-1.13615\), \(\widehat{\theta}_{2}=-0.66875\). Jika dibandingkan dengan data \(Y\) yang dibangkitkan yaitu ARIMA(1,2,2) dengan parameter \(\mu=0.35\), \(\phi=0.75\), \(\theta_{1}=-0.65\), \(\theta_{2}=0.40\). sigma^2 estimated as 0.9377 adalah nilai dugaan bagi \(\sigma_{e}^2\). Model terbaik yang diperoleh yaitu model \(ARIMA(1,2,2)\) adalah model yang sama dengan model awal yang digunakan untuk membangkitkan data. Oleh karena itu, model terbaik tersebut selanjutnya bisa digunakan untuk peramalan.
Model \(ARIMA(1,2,2)\), maka \(Y_{t}\) diperoleh dari penjabaran operator backshift sehingga untuk model \(ARIMA(1,2,2)\):
\[\phi_{p} (1-B)^{d} Y_{t}=\mu+\theta_{q}(B)e_{t}\] \[(1-\phi_{1}(B))(1-B)^{2} Y_{t}=\mu+\theta_{2}(B)e_{t}\] \[(1-\phi_{1}(B))(1-2B+B^{2}) Y_{t}=\mu+(1-\theta_{1}(B)-\theta_{2}B^{2})e_{t}\] \[(1-2B+B^{2}-\phi_{1}B+2\phi_{1}B^{2}-\phi_{1}B^{3}) Y_{t}=\mu+e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}\] \[(1-(2+\phi_{1})B+(1+2\phi_{1})B^{2}-\phi_{1}B^{3}) Y_{t}=\mu+e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}\] \[Y_{t}-(2+\phi_{1})Y_{t-1}+(1+2\phi_{1})Y_{t-2}-\phi_{1}Y_{t-3} =\mu+e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}\] \[Y_{t} =\mu+(2+\phi_{1})Y_{t-1}-(1+2\phi_{1})Y_{t-2}+\phi_{1}Y_{t-3}+e_{t}-\theta_{1}e_{t-1}-\theta_{2}e_{t-2}\] \[Y_{t} =1.23001+(2+0.29078)Y_{t-1}-(1+2(0.29078))Y_{t-2}+0.29078Y_{t-3}+e_{t}-(-1.13615)e_{t-1}-(0.66875)e_{t-2}\] \[Y_{t} =1.23001+2.29078Y_{t-1}-1.58156Y_{t-2}+0.29078Y_{t-3}+e_{t}+1.13615e_{t-1}+0.66875e_{t-2}\]
Model terbaik ini selanjutnya dilakukan uji hipotesis untuk mengetahui signifikansi sebuah parameter.
\(H_{0}\) : parameter \(= 0\)
\(H_{1}\) : parameter \(\neq\) \(0\)
Pengambilan keputusan yaitu Tolak \(H_{0}\) jika \(|t_{hitung}| > t_{\frac{\alpha}{2},(n-1)}\)
\(\alpha\) yang digunakan sebesar \(10\%\)
T-hitung parameter dugaan model \(ARIMA(1,2,2)\) \[t_{hitung}=\frac{parameter estimasi}{SE parameter estimasi}\]
intercept_mu <- 1.23001/0.28903
intercept_mu## [1] 4.255648ar_phi <- 0.29078/0.09791
ar_phi## [1] 2.96987ma_theta1 <- -1.13615/0.08342
ma_theta1## [1] -13.61964ma_theta2 <- -0.66875/0.06910
ma_theta2## [1] -9.678003prmtrdugaan <- c("mu", "phi", "theta(1)", "theta(2)")
thitung <- c(abs(intercept_mu), abs(ar_phi), abs(ma_theta1), abs(ma_theta2))
ttabel <- c("1.97369144", "1.97369144","1.97369144", "1.97369144")
keputusan <- c("Signifikan", "Signifikan", "Signifikan", "Signifikan")
tksignimodel <- cbind(prmtrdugaan, thitung, ttabel, keputusan)
colnames(tksignimodel) <- c("Parameter dugaan", "T-Hitung", "T-Tabel", "Keputusan")
tksignimodel <- as.data.frame(tksignimodel)
tksignimodel## Parameter dugaan T-Hitung T-Tabel Keputusan
## 1 mu 4.25564820260873 1.97369144 Signifikan
## 2 phi 2.96987028904096 1.97369144 Signifikan
## 3 theta(1) 13.6196355789978 1.97369144 Signifikan
## 4 theta(2) 9.67800289435601 1.97369144 SignifikanBerdasarkan tabel di atas maka pengambilan keputusan dari semua paramater dugaan dari model \(ARIMA(1,2,2)\) adalah tolak \(H_{0}\) artinya semua paremeter dugaannya nyata atau signifikan terhadap model.
lmtest::coeftest(arima122)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.290781 0.097905 2.9700 0.002978 **
## ma1 1.136152 0.083415 13.6205 < 2.2e-16 ***
## ma2 0.668750 0.069099 9.6781 < 2.2e-16 ***
## intercept 1.230007 0.289034 4.2556 2.085e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Berdasarkan tabel di atas maka pengambilan keputusan dari semua paramater dugaan dari model \(ARIMA(1,0,0)\) semua parameternya signifikan.
Diagnostic model
ARIMA122diag <- stats::arima(yt.arima122.diff2, order = c(1,2,2), method = "ML")
checkresiduals(ARIMA122diag)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,2,2)
## Q* = 36.867, df = 7, p-value = 4.97e-06
##
## Model df: 3. Total lags used: 10Berdasarkan pl0t di atas terlihat bahwa sisaan 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 <- arima122$residuals# Uji formal normalitas data
ks.test(sisaan,"pnorm")##
## One-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.032685, p-value = 0.9926
## alternative hypothesis: two-sidedHasil pengujian, Normalitas Data:
\(H_{0}\) : sisaan menyebar normal
\(H_{1}\) : sisaan tidak mengikuti sebaran normal
Hasil : \(p-value=0.99264 > \alpha=0.05\) yang berarti TAK TOLAK \(H_{0}\). Sisaan menyebar normal
# Uji nilai tengah sisaan
t.test(sisaan, mu = 0, alternative = "two.sided")##
## One Sample t-test
##
## data: sisaan
## t = 0.043935, df = 172, p-value = 0.965
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.1424965 0.1489845
## sample estimates:
## mean of x
## 0.003243983Nilai Tengah Sisaan
\(H_{0}:\mu = 0\)
\(H_{1}:\mu \neq 0\)
Hasil : \(p-value=0.96501 > \alpha=0.05\) yang berarti TAK TOLAK \(H_{0}\). Nilai tengah sisaan sama dengan \(0\)
# Uji autokorelasi
Box.test(sisaan, lag = 23 ,type = "Ljung")##
## Box-Ljung test
##
## data: sisaan
## X-squared = 13.053, df = 23, p-value = 0.9508Autokorelasi
\(H_{0}\) : tidak ada autokorelasi
\(H_{1}\) : terdapat autokorelasi
Hasil : \(p-value=0.95083 > \alpha=0.05\) yang berarti TAK TOLAK \(H_{0}\). Tidak terdapat gejala autokorelasi
Kesimpulan : Semua asumsi terpenuhi.
Soal 3.
Penerapan pada Data AKTUAL
Melalui program R, kerjakan: Exercise 5.11 (Montgomery, hlm. 415):
5.11 Table B.2 contains data on pharmaceutical product sales.
Membaca Data
Data yang digunakan adalah pada buku Montgomery, halaman 415. Data ini berisi total penjualan produk farmasi perminggu dalam ratusan ribu.
data<-read_xlsx("D:/S2/IPB/STA542 Analisis Deret Waktu/6/data.xlsx")
data <- data$`Sales (In Thousands)`Visualisasi Data
Visualisasi data dari total penjualan produk farmasi perminggu dalam ratusan ribu tersaji pada plot berikut:
data.ts <- ts(data) #form a time-series object
data.ts## Time Series:
## Start = 1
## End = 120
## Frequency = 1
## [1] 10618.1 10537.9 10209.3 10553.0 9934.9 10534.5 10196.5 10511.8 10089.6
## [10] 10371.2 10239.4 10472.4 10827.2 10640.8 10517.8 10154.2 9969.2 10260.4
## [19] 10737.0 10430.0 10689.0 10430.4 10002.4 10135.7 10096.2 10288.7 10289.1
## [28] 10589.9 10551.9 10208.3 10334.5 10480.1 10387.6 10202.6 10219.3 10382.7
## [37] 10820.5 10358.7 10494.6 10497.6 10431.5 10447.8 10684.4 10176.5 10616.0
## [46] 10627.7 10684.0 10246.7 10265.0 10090.4 9881.1 10449.7 10276.3 10175.2
## [55] 10212.5 10395.5 10545.9 10635.7 10265.2 10551.6 10538.2 10286.2 10171.3
## [64] 10393.1 10162.3 10164.5 10327.0 10365.1 10755.9 10463.6 10080.5 10479.6
## [73] 9980.9 10039.2 10246.1 10368.0 10446.3 10535.3 10786.9 9975.8 10160.9
## [82] 10422.1 10757.2 10463.8 10307.0 10134.7 10207.7 10488.0 10262.3 10785.9
## [91] 10375.4 10123.4 10462.7 10205.5 10522.7 10253.2 10428.7 10615.8 10417.3
## [100] 10445.4 10690.6 10271.8 10524.8 9815.0 10398.5 10553.1 10655.8 10199.1
## [109] 10416.6 10391.3 10210.1 10352.5 10423.8 10519.3 10596.7 10650.0 10741.6
## [118] 10246.0 10354.4 10155.4ts.plot(data.ts, col="blue", ylab = "Sales (In Thousands)", xlab = "Time Period (Week)")
title(main = "Time Series Plot of Data Sales", cex.sub = 0.8)
points(data.ts, pch = 20, col = "blue")
- Fit an ARIMA model to this time series, excluding the last 10 observations. Investigate model adequacy. Explain how this model would be used for forecasting.
Jawaban poin a.
Data yang akan digunakan untuk menentukan model terbaik berdasarkan kandidat model akan dilakukan splitting data terlebih dahulu yaitu dengan membagi data menjadi dua bagian yang pertama data ke-1 hingga data ke-110 sebagai data training dan data ke-111 hingga data ke-120 menjadi data testing.
# Training set
train.sales.ts = subset(data.ts, start=1,end=110)
# Testing set
test.sales.ts = subset(data.ts, start=111,end=120)Visualisasi data training yang akan digunakan untuk mencari model terbaik
ts.plot(train.sales.ts, col="blue", ylab = "Sales (In Thousands)", xlab = "Time Period (Week)")
title(main = "Time Series Plot of Train Data Sales", cex.sub = 0.8)
points(data.ts, pch = 20, col = "blue")
Model ARIMA terbaik dari data total penjualan produk farmasi perminggu dalam ratusan ribu akan doperoleh dengan melalui beberapa tahapan 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.sales.ts, alternative = "stationary", k=trunc((length(train.sales.ts)-1)^(1/3)))## Warning in adf.test(train.sales.ts, alternative = "stationary", k =
## trunc((length(train.sales.ts) - : p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: train.sales.ts
## Dickey-Fuller = -5.9291, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryBerdasarkan 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 digunakan merupakan data yang stasioner
2. Identifikasi kandidat model Identifikasi kandidat model diperoleh berdasarkan nilai \(p\) dan \(q\) dimana nilai \(d=0\)
acf(train.sales.ts, lag.max = 20, col = "blue")
Berdasarkan plot ACF terlihat bahwa plot cuts off setelah lag ke-3, \(ARIMA(0,0,3)\)
pacf(train.sales.ts, lag.max = 20, col = "blue")
Berdasarkan plot PACF di atas terlihat bahwa plot PACF menunjukkan model \(ARIMA(3,0,0)\)
eacf(train.sales.ts)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o x o o o o o o o o o o o
## 1 o o o o o o o o o o o o o o
## 2 o x o o o o o o o o o o o 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 o o
## 5 o o x o o o o o o o o o o o
## 6 x o o o o o o o o o o o o o
## 7 x x o o x o o o o o o o o o
Berdasarkan plot ACF, PACF, dan EACF dapat ditentukan beberapa kandidat model yaitu: \(ARIMA(0,0,3)\), \(ARIMA(3,0,0)\), \(ARIMA(3,0,3)\), \(ARIMA(1,0,2)\), \(ARIMA(2,0,2)\), dan \(ARIMA(3,0,1)\)
3. Pemilihan model terbaik berdasarkan nilai AIC terkecil dari kandidat model
# ARIMA(0,0,3)
arima003 <- arima(train.sales.ts, order = c(0,0,3), include.mean = TRUE, method = "ML")
arima003##
## Call:
## arima(x = train.sales.ts, order = c(0, 0, 3), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ma1 ma2 ma3 intercept
## 0.0141 -0.0518 -0.2326 10374.4514
## s.e. 0.1002 0.1171 0.0967 14.9337
##
## sigma^2 estimated as 45079: log likelihood = -745.56, aic = 1499.12# ARIMA(3,0,0)
arima300 <- arima(train.sales.ts, order = c(3,0,0), include.mean = TRUE, method = "ML")
arima300##
## Call:
## arima(x = train.sales.ts, order = c(3, 0, 0), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 intercept
## 0.0904 0.0086 -0.2312 10374.9540
## s.e. 0.0927 0.0929 0.0926 17.9023
##
## sigma^2 estimated as 44745: log likelihood = -745.15, aic = 1498.3# ARIMA(3,0,3)
arima303 <- arima(train.sales.ts, order = c(3,0,3), include.mean = TRUE, method = "ML")
arima303##
## Call:
## arima(x = train.sales.ts, order = c(3, 0, 3), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 intercept
## 1.7166 -1.4818 0.6168 -1.7132 1.3812 -0.6680 10370.9490
## s.e. 0.3399 0.4689 0.2187 0.3450 0.4829 0.2248 3.4577
##
## sigma^2 estimated as 41349: log likelihood = -742.04, aic = 1498.07# ARIMA(1,0,2)
arima102 <- arima(train.sales.ts, order = c(1,0,2), include.mean = TRUE, method = "ML")
arima102##
## Call:
## arima(x = train.sales.ts, order = c(1, 0, 2), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ma1 ma2 intercept
## -0.5789 0.6900 0.1485 10375.7967
## s.e. 0.5456 0.5358 0.1023 24.0048
##
## sigma^2 estimated as 46866: log likelihood = -747.63, aic = 1503.26# ARIMA(2,0,2)
arima202 <- arima(train.sales.ts, order = c(2,0,2), include.mean = TRUE, method = "ML")
arima202##
## Call:
## arima(x = train.sales.ts, order = c(2, 0, 2), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ar2 ma1 ma2 intercept
## 1.1590 -0.3063 -1.1583 0.1583 10370.7951
## s.e. 0.2954 0.2562 0.2999 0.2985 3.2127
##
## sigma^2 estimated as 42722: log likelihood = -743.85, aic = 1497.7# ARIMA(3,0,1)
arima301 <- arima(train.sales.ts, order = c(3,0,1), include.mean = TRUE, method = "ML")
arima301##
## Call:
## arima(x = train.sales.ts, order = c(3, 0, 1), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ma1 intercept
## 0.6498 -0.0437 -0.2197 -0.6094 10374.0433
## s.e. 0.2868 0.1149 0.1105 0.2950 12.8663
##
## sigma^2 estimated as 43484: log likelihood = -743.64, aic = 1497.28AICsoal3 <- c(1499.12, 1498.3, 1498.07, 1503.26, 1497.7, 1497.28)
KMARIMAsoal3 <- c("ARIMA(0,0,3)", "ARIMA(3,0,0)","ARIMA(3,0,3)", "ARIMA(1,0,2)", "ARIMA(2,0,2)", "ARIMA(3,0,1)")
compmodelsoal3 <- cbind(KMARIMAsoal3, AICsoal3)
colnames(compmodelsoal3) <- c("Kandidat Model", "Nilai AIC")
compmodelsoal3 <- as.data.frame(compmodelsoal3)
compmodelsoal3## Kandidat Model Nilai AIC
## 1 ARIMA(0,0,3) 1499.12
## 2 ARIMA(3,0,0) 1498.3
## 3 ARIMA(3,0,3) 1498.07
## 4 ARIMA(1,0,2) 1503.26
## 5 ARIMA(2,0,2) 1497.7
## 6 ARIMA(3,0,1) 1497.28Model terbaik diperoleh Berdasarkan nilai AIC terkecil dari kandidat model. Oleh karena itu, Model terbaik yang diperoleh yaitu ARIMA(3,0,1)
Model terbaik ini selanjutnya dilakukan uji hipotesis untuk mengetahui signifikansi sebuah parameter.
\(H_{0}\) : parameter = \(0\)
\(H_{1}\) : parameter $$ \(0\)
Pengambilan keputusan yaitu Tolak \(H_{0}\) jika \(|t_{hitung}| > t_{\frac{\alpha}{2},(n-1)}\)
\(\alpha\) yang digunakan sebesar \(10\%\)
T-hitung parameter dugaan model \(ARIMA(1,2,2)\) \[t_{hitung}=\frac{parameter estimasi}{SE parameter estimasi}\]
arima301 #model terbaik##
## Call:
## arima(x = train.sales.ts, order = c(3, 0, 1), include.mean = TRUE, method = "ML")
##
## Coefficients:
## ar1 ar2 ar3 ma1 intercept
## 0.6498 -0.0437 -0.2197 -0.6094 10374.0433
## s.e. 0.2868 0.1149 0.1105 0.2950 12.8663
##
## sigma^2 estimated as 43484: log likelihood = -743.64, aic = 1497.28intercept_mu_ <- 10374.04332/12.86633
intercept_mu_## [1] 806.2939ar_phi_1 <- 0.64980/0.28679
ar_phi_1## [1] 2.265769ar_phi_2 <- -0.04367/0.11492
ar_phi_2## [1] -0.3800035ar_phi_3 <- -0.21975/0.11049
ar_phi_3## [1] -1.988868ma_theta1 <- -0.60943/0.29502
ma_theta1## [1] -2.065724prmtrdugaann <- c("mu", "phi(1)", "phi(2)", "phi(3)", "theta(1)")
thitungg <- c(abs(intercept_mu_), abs(ar_phi_1), abs(ar_phi_2), abs(ar_phi_3), abs(ma_theta1))
ttabell <- c("1.98196749", "1.98196749","1.98196749", "1.98196749", "1.98196749")
keputusann <- c("Signifikan", "Signifikan", "Tidak Signifikan", "Signifikan", "Signifikan")
tksignimodel301 <- cbind(prmtrdugaann, thitungg, ttabell, keputusann)
colnames(tksignimodel301) <- c("Parameter dugaan", "T-Hitung", "T-Tabel", "Keputusan")
tksignimodel301 <- as.data.frame(tksignimodel301)
tksignimodel301## Parameter dugaan T-Hitung T-Tabel Keputusan
## 1 mu 806.293894218476 1.98196749 Signifikan
## 2 phi(1) 2.26576937829074 1.98196749 Signifikan
## 3 phi(2) 0.380003480682214 1.98196749 Tidak Signifikan
## 4 phi(3) 1.98886777083899 1.98196749 Signifikan
## 5 theta(1) 2.06572435767067 1.98196749 SignifikanBerdasarkan tabel di atas maka pengambilan keputusan dari paramater dugaan dari model \(ARIMA(3,0,1)\) adalah tolak \(H_{0}\) artinya semua paremeter dugaannya nyata atau signifikan terhadap model, kecuali untuk parameter ketiga yaitu \(\phi_{2}\) yang menunjukkan AR(2) tidak signifikan terhadap model.
4. Diagnostic model
ARIMA301diag <- stats::arima(train.sales.ts, order = c(3,0,1), method = "ML")
checkresiduals(ARIMA301diag)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(3,0,1) with non-zero mean
## Q* = 1.8854, df = 5, p-value = 0.8648
##
## Model df: 5. Total lags used: 10Berdasarkan pl0t di atas terlihat bahwa sisaan 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 <- arima301$residuals# Uji formal normalitas data
ks.test(sisaan,"pnorm")##
## One-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.5, p-value < 2.2e-16
## alternative hypothesis: two-sidedHasil pengujian, Normalitas Data:
\(H_{0}\) : sisaan menyebar normal
\(H_{1}\) : sisaan tidak mengikuti sebaran normal
Hasil : \(p-value=0.99264 > \alpha=0.05\) yang berarti TAK TOLAK \(H_{0}\). Sisaan menyebar normal
# Uji nilai tengah sisaan
t.test(sisaan, mu = 0, alternative = "two.sided")##
## One Sample t-test
##
## data: sisaan
## t = 0.037872, df = 109, p-value = 0.9699
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -38.83003 40.34288
## sample estimates:
## mean of x
## 0.7564248Nilai Tengah Sisaan
\(H_{0}:\mu = 0\)
\(H_{1}:\mu \neq 0\)
Hasil : \(p-value=0.96986 > \alpha=0.05\) yang berarti TAK TOLAK $H_{0}. Nilai tengah sisaan sama dengan \(0\)
# Uji autokorelasi
Box.test(sisaan, lag = 23 ,type = "Ljung")##
## Box-Ljung test
##
## data: sisaan
## X-squared = 8.8403, df = 23, p-value = 0.9965Autokorelasi
\(H_{0}\) : tidak ada autokorelasi
\(H_{1}\) : terdapat autokorelasi
Hasil : \(p-value=0.99647 > \alpha=0.05\) yang berarti TAK TOLAK \(H_{0}\). Tidak terdapat gejala autokorelasi
Kesimpulan : Asumsi terpenuhi, kecuali sisaan tidak menyebar normal.
- Forecast the last 10 observations.
Jawaban poin b.
ARIMAdata <- Arima(train.sales.ts, order=c(3,0,1), method="ML")
forecastARIMA301 <- forecast(ARIMAdata)
hasil_forecastARIMA301 <- as.data.frame(forecastARIMA301)
hasil_forecastARIMA301## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 111 10371.02 10097.49 10644.55 9952.692 10789.35
## 112 10361.97 10088.22 10635.72 9943.305 10780.64
## 113 10362.54 10088.75 10636.33 9943.808 10781.27
## 114 10367.76 10086.66 10648.86 9937.849 10797.67
## 115 10373.12 10088.65 10657.58 9938.063 10808.17
## 116 10376.24 10090.72 10661.76 9939.579 10812.91
## 117 10376.89 10091.37 10662.41 9940.230 10813.56
## 118 10376.00 10090.29 10661.72 9939.037 10812.97
## 119 10374.71 10088.73 10660.69 9937.338 10812.08
## 120 10373.76 10087.68 10659.84 9936.242 10811.29Visualisasi hasil forecast
plot(forecastARIMA301)
Akurasi hasil peramalan
data_forecastARIMA301 <- hasil_forecastARIMA301$`Point Forecast`
data_forecastARIMA301 <- as.data.frame(data_forecastARIMA301)
head(data_forecastARIMA301)## data_forecastARIMA301
## 1 10371.02
## 2 10361.97
## 3 10362.54
## 4 10367.76
## 5 10373.12
## 6 10376.24Data Aktual : Data hasil peramalan
dataaktual <- as.data.frame(test.sales.ts)
compmodelforecast <- cbind(dataaktual, data_forecastARIMA301)
colnames(compmodelforecast) <- c("Data Aktual", "Data hasil peramalan ARIMA(3,0,1")
compmodelforecast <- as.data.frame(compmodelforecast)
compmodelforecast## Data Aktual Data hasil peramalan ARIMA(3,0,1
## 1 10210.1 10371.02
## 2 10352.5 10361.97
## 3 10423.8 10362.54
## 4 10519.3 10367.76
## 5 10596.7 10373.12
## 6 10650.0 10376.24
## 7 10741.6 10376.89
## 8 10246.0 10376.00
## 9 10354.4 10374.71
## 10 10155.4 10373.76# Akurasi
akurasi.arima301 <- accuracy(data_forecastARIMA301$data_forecastARIMA301,dataaktual$x)
akurasi.arima301## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 53.57804 193.8149 161.3917 0.4821229 1.538676 0.3254855 1.007241- In Exercise 4.12, you were asked to use simple exponential smoothing with \(\lambda = 0.1\) to smooth the data, and to forecast the last 10 observations. Compare the ARIMA and exponential smoothing forecasts. Which forecasting method do you prefer?
Jawaban poin c.
Membaca Data
Data yang digunakan adalah pada buku Montgomery, halaman 415. Data ini berisi total penjualan produk farmasi perminggu dalam ratusan ribu.
data<-read_xlsx("D:/S2/IPB/STA542 Analisis Deret Waktu/6/data.xlsx")
data <- data$`Sales (In Thousands)`Visualisasi Data
Visualisasi data dari total penjualan produk farmasi perminggu dalam ratusan ribu tersaji pada plot berikut:
data.ts <- ts(data) #form a time-series object
data.ts## Time Series:
## Start = 1
## End = 120
## Frequency = 1
## [1] 10618.1 10537.9 10209.3 10553.0 9934.9 10534.5 10196.5 10511.8 10089.6
## [10] 10371.2 10239.4 10472.4 10827.2 10640.8 10517.8 10154.2 9969.2 10260.4
## [19] 10737.0 10430.0 10689.0 10430.4 10002.4 10135.7 10096.2 10288.7 10289.1
## [28] 10589.9 10551.9 10208.3 10334.5 10480.1 10387.6 10202.6 10219.3 10382.7
## [37] 10820.5 10358.7 10494.6 10497.6 10431.5 10447.8 10684.4 10176.5 10616.0
## [46] 10627.7 10684.0 10246.7 10265.0 10090.4 9881.1 10449.7 10276.3 10175.2
## [55] 10212.5 10395.5 10545.9 10635.7 10265.2 10551.6 10538.2 10286.2 10171.3
## [64] 10393.1 10162.3 10164.5 10327.0 10365.1 10755.9 10463.6 10080.5 10479.6
## [73] 9980.9 10039.2 10246.1 10368.0 10446.3 10535.3 10786.9 9975.8 10160.9
## [82] 10422.1 10757.2 10463.8 10307.0 10134.7 10207.7 10488.0 10262.3 10785.9
## [91] 10375.4 10123.4 10462.7 10205.5 10522.7 10253.2 10428.7 10615.8 10417.3
## [100] 10445.4 10690.6 10271.8 10524.8 9815.0 10398.5 10553.1 10655.8 10199.1
## [109] 10416.6 10391.3 10210.1 10352.5 10423.8 10519.3 10596.7 10650.0 10741.6
## [118] 10246.0 10354.4 10155.4ts.plot(data.ts, col="blue", ylab = "Sales (In Thousands)", xlab = "Time Period (Week)")
title(main = "Time Series Plot of Data Sales", cex.sub = 0.8)
points(data.ts, pch = 20, col = "blue")
Data yang akan digunakan untuk dimodelkan dengan metode Simple Exponential Smoothing kemudian akan dilakukan terlebih dahulu splitting data yaitu dengan membagi data menjadi dua bagian yang pertama data ke-1 hingga data ke-110 sebagai data training dan data ke-111 hingga data ke-120 menjadi data testing. Splitting data ini dilakukan sama seperti dalam memodelkan data pada metode ARIMA.
# Training set
train.sales.ts = subset(data.ts, start=1,end=110)
# Testing set
test.sales.ts = subset(data.ts, start=111,end=120)Visualisasi data training yang akan digunakan untuk metode Simple Exponential Smoothing dengan \(\lambda = 0.1\).
ts.plot(train.sales.ts, col="blue", ylab = "Sales (In Thousands)", xlab = "Time Period (Week)")
title(main = "Time Series Plot of Train Data Sales", cex.sub = 0.8)
points(data.ts, pch = 20, col = "blue")
data.ses<-HoltWinters(data.ts,alpha=0.1, beta=FALSE, gamma=FALSE)
data.ses## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = data.ts, alpha = 0.1, beta = FALSE, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.1
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 10407.6Akurasi
SSE <- data.ses$SSE
SSE## [1] 6338093MSE <- SSE/110
MSE## [1] 57619.03Hasil peramalan dengan menggunakan metode Simple Exponential Smoothing dengan \(\lambda = 0.1\)
forecastresultsess <- forecast(data.ses)
forecastresultses <- as.data.frame(forecastresultsess)
forecastresultses## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 121 10407.6 10111.46 10703.74 9954.697 10860.51
## 122 10407.6 10109.99 10705.22 9952.438 10862.77
## 123 10407.6 10108.52 10706.69 9950.191 10865.01
## 124 10407.6 10107.05 10708.15 9947.954 10867.25
## 125 10407.6 10105.60 10709.61 9945.728 10869.48
## 126 10407.6 10104.15 10711.05 9943.513 10871.69
## 127 10407.6 10102.71 10712.50 9941.308 10873.90
## 128 10407.6 10101.27 10713.93 9939.114 10876.09
## 129 10407.6 10099.85 10715.36 9936.930 10878.27
## 130 10407.6 10098.42 10716.78 9934.756 10880.45Visualisasi hasil forecast
plot(forecastresultsess)
Data Aktual : Data hasil peramalan
dataaktualses <- as.data.frame(test.sales.ts)
compmodelforecastses <- cbind(dataaktual, forecastresultses$`Point Forecast`)
colnames(compmodelforecastses) <- c("Data Aktual", "Data hasil peramalan SES dengan Lambda=1")
compmodelforecastses <- as.data.frame(compmodelforecastses)
compmodelforecastses## Data Aktual Data hasil peramalan SES dengan Lambda=1
## 1 10210.1 10407.6
## 2 10352.5 10407.6
## 3 10423.8 10407.6
## 4 10519.3 10407.6
## 5 10596.7 10407.6
## 6 10650.0 10407.6
## 7 10741.6 10407.6
## 8 10246.0 10407.6
## 9 10354.4 10407.6
## 10 10155.4 10407.6# Akurasi
akurasi.ses <- accuracy(forecastresultses$`Point Forecast`,dataaktual$x)
akurasi.ses## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 17.3778 187.9902 161.3 0.1346043 1.542825 0.3328423 0.9559034akurasi.arima301## ME RMSE MAE MPE MAPE ACF1 Theil's U
## Test set 53.57804 193.8149 161.3917 0.4821229 1.538676 0.3254855 1.007241Perbandingan Peramalan Model ARIMA dan Exponential Smoothing
modelakursi <- c("ME", "RMSE", "MAE", "MPE", "MAPE")
ModelARIMA301aku <- c(53.57804, 193.8149, 161.3917, 0.4821229, 1.538676)
ModelSES <- c(17.3778, 187.9902, 161.3, 0.1346043, 1.542825)
Akurasiall <- cbind(modelakursi, ModelARIMA301aku, ModelSES)
colnames(Akurasiall) <- c("Akurasi", "Model ARIMA(3,0,1)", "Model SES Lambda=0.1")
Akurasiall <- as.data.frame(Akurasiall)
Akurasiall## Akurasi Model ARIMA(3,0,1) Model SES Lambda=0.1
## 1 ME 53.57804 17.3778
## 2 RMSE 193.8149 187.9902
## 3 MAE 161.3917 161.3
## 4 MPE 0.4821229 0.1346043
## 5 MAPE 1.538676 1.542825Berdasarkan hasil peramalan dengan menggunakan metode ARIMA yang diperoleh model terbaik yaitu \(ARIMA(3,0,1)\) dan metode simple exponential smoothing dengan \(\lambda=0.1\) diperoleh model yang terbaik yaitu ARIMA(3,0,1). Penilihan model terbaik ini didapatkan berdasarkan nilai akurasi MAPE yang terkecil anatara kedua model yaitu sebesar 1.538676.
Jadi, motode terbaik untuk meramalkan data total penjualan produk farmasi perminggu dalam ratusan ribu adalah dengan metode ARIMA dan diperoleh model terbaiknya yaitu \(ARIMA(3,0,1)\) dimana nilai \(p=3\), \(d=0\), dan \(q=1\).
- How would prediction intervals be obtained for the ARIMA forecasts?
Jawaban poin d.
Prediksi interval untuk peramalan model ARIMA yaitu:
fitarima301 <- Arima(train.sales.ts, order=c(3,0,1), method="ML")
fitarima301## Series: train.sales.ts
## ARIMA(3,0,1) with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 ma1 mean
## 0.6498 -0.0437 -0.2197 -0.6094 10374.0433
## s.e. 0.2868 0.1149 0.1105 0.2950 12.8663
##
## sigma^2 estimated as 45555: log likelihood=-743.64
## AIC=1499.28 AICc=1500.1 BIC=1515.49hasil_forecastARIMA301## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 111 10371.02 10097.49 10644.55 9952.692 10789.35
## 112 10361.97 10088.22 10635.72 9943.305 10780.64
## 113 10362.54 10088.75 10636.33 9943.808 10781.27
## 114 10367.76 10086.66 10648.86 9937.849 10797.67
## 115 10373.12 10088.65 10657.58 9938.063 10808.17
## 116 10376.24 10090.72 10661.76 9939.579 10812.91
## 117 10376.89 10091.37 10662.41 9940.230 10813.56
## 118 10376.00 10090.29 10661.72 9939.037 10812.97
## 119 10374.71 10088.73 10660.69 9937.338 10812.08
## 120 10373.76 10087.68 10659.84 9936.242 10811.29interval.prediction <- hasil_forecastARIMA301[,-1]
interval.prediction## Lo 80 Hi 80 Lo 95 Hi 95
## 111 10097.49 10644.55 9952.692 10789.35
## 112 10088.22 10635.72 9943.305 10780.64
## 113 10088.75 10636.33 9943.808 10781.27
## 114 10086.66 10648.86 9937.849 10797.67
## 115 10088.65 10657.58 9938.063 10808.17
## 116 10090.72 10661.76 9939.579 10812.91
## 117 10091.37 10662.41 9940.230 10813.56
## 118 10090.29 10661.72 9939.037 10812.97
## 119 10088.73 10660.69 9937.338 10812.08
## 120 10087.68 10659.84 9936.242 10811.29Visualisasi
plot(hasil_forecastARIMA301$`Point Forecast`, type="n", ylim=range(hasil_forecastARIMA301$`Lo 80`, hasil_forecastARIMA301$`Hi 80`))
polygon(c(time(hasil_forecastARIMA301$`Point Forecast`),rev(time(hasil_forecastARIMA301$`Point Forecast`))), c(hasil_forecastARIMA301$`Hi 80`,rev(hasil_forecastARIMA301$`Lo 80`)),
col=rgb(0,0,0.6,0.2), border=FALSE)
lines(hasil_forecastARIMA301$`Point Forecast`)
lines(fitted(fitarima301),col='red')
out <- (hasil_forecastARIMA301$`Point Forecast`< hasil_forecastARIMA301$`Lo 80` | hasil_forecastARIMA301$`Point Forecast` > hasil_forecastARIMA301$`Hi 80`)
plot(hasil_forecastARIMA301$`Point Forecast`, type="n", ylim=range(hasil_forecastARIMA301$`Lo 95`, hasil_forecastARIMA301$`Hi 95`))
polygon(c(time(hasil_forecastARIMA301$`Point Forecast`),rev(time(hasil_forecastARIMA301$`Point Forecast`))), c(hasil_forecastARIMA301$`Hi 95`,rev(hasil_forecastARIMA301$`Lo 95`)),
col=rgb(0,0,0.6,0.2), border=FALSE)
lines(hasil_forecastARIMA301$`Point Forecast`)
lines(fitted(fitarima301),col='red')
out <- (hasil_forecastARIMA301$`Point Forecast`< hasil_forecastARIMA301$`Lo 95` | hasil_forecastARIMA301$`Point Forecast` > hasil_forecastARIMA301$`Hi 95`)
out## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSEBerdasarkan kedua plot yaitu plot interval prediksi 80% dan plot interval prediksi 95% menunjukkan hasil yang hampir sama hanya berbeda di rentang selangnya. Interval prediksi 95% memiliki rentang yang lebih luas dari pada interval prediksi 80%.