Stage 1: Model Identification Contoh berikut menunjukkan prosedur untuk mencocokkan seri data dengan komponen musiman akan menggunakan data HIV yang terjadi di Indonesia tahun 1990-2019.
library(tseries)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(forecast)
library(readxl)data_lopo <- read_excel("~/.A. KULIAH/Semester 7/Analis Runtun Waktu/Data lopo paper.xlsx")#check the structure of data set
str(data_lopo)## tibble [30 × 5] (S3: tbl_df/tbl/data.frame)
## $ No : num [1:30] 1 2 3 4 5 6 7 8 9 10 ...
## $ Negara : chr [1:30] "Indonesia" "Indonesia" "Indonesia" "Indonesia" ...
## $ Kode Negara : chr [1:30] "IDN" "IDN" "IDN" "IDN" ...
## $ Tahun : num [1:30] 1990 1991 1992 1993 1994 ...
## $ Jumlah Penderita: num [1:30] 2747 2192 1121 1582 1801 ...data_lopo.ts <- ts(data_lopo$`Jumlah Penderita`, start = c(1990, 1),
end = c(2019, 12), frequency = 12)options(scipen = 2) #Suppress the scientific value for ylab
plot(data_lopo.ts, type ="o", pch = 16,
ylab = "Data No. 1", panel.first = grid ())
Seri data diselidiki dengan memplot data infeksi HIV terhadap waktu. Grafik tersebut menunjukkan bahwa rangkaian tersebut memiliki puncak dan palung yang teratur tidak stasioner karena terdapat komponen musiman dalam rangkaian tersebut. Dari seri asli, puncak reguler dan melalui pola musiman yang ditunjukkan dalam series.
par(mfrow = c(2, 1)) # 2 graphs in a single device
acf(data_lopo.ts, 36, main = "ACF DATA no. 1 HIV",
panel.first = grid (), xaxt = "n")
axis(1, at = 0:36/12, labels=0:36) # xaxis re-label
pacf(data_lopo.ts, 36, main = "PACF DATA no. 1 HIV",
panel.first = grid (), xaxt = "n")
axis(1, at = 0:36/12, labels=0:36) # xaxis re-label
Berdasarkan plot ACF, terdapat pola seperti gelombang yang menegaskan
komponen musiman dalam rangkaian infeksi HIV dan sedikit perluasan pada
plot PACF. Pola ini mencirikan seri dengan efek musiman.
Karena ada komponen musiman pada seri, pembedaan musim urutan pertama dilakukan. Perhatikan bahwa kami menambahkan argumen tambahan dalam fungsi diff(), lag = 12 untuk melakukan pembedaan musiman.
data_lopo_sea.ts <- diff(data_lopo.ts, lag = 12)par(mfrow = c(1, 1)) # 1 graph in a single device
plot(data_lopo_sea.ts, type ="o", pch = 16,
ylab = "DATA No. 1 HIV",
main = "Seasonal Differencing", panel.first = grid ())
par (mfrow = c(2, 1)) # 2 graphs in a single device
acf(data_lopo_sea.ts, 36,
main = "ACF for DATA no.1 (Seasonal Diff)",
panel.first = grid (), xaxt = "n")
axis(1, at = 0:36/12, labels=0:36)
pacf(data_lopo_sea.ts, 36, main = "PACF for DATA no.1 (Seasonal Diff)", panel.first = grid (), xaxt = "n")
axis(1, at = 0:36/12, labels=0:36)
Penjelasan: Setelah melakukan pembedaan musiman pertama, plot ACF dan
PACF menunjukkan karakteristik peluruhan dan bergelombang, yang
menunjukkan bahwa deret tersebut masih belum stasioner. Plot ACF
menunjukkan pola yang menurun dan masih menunjukkan pola gelombang
(tetapi tidak begitu jelas dibandingkan dengan data asli) sementara PACF
menghasilkan beberapa lonjakan yang signifikan terutama pada lag
pertama.
Ini menunjukkan bahwa kita perlu melakukan pembedaan non-musiman pertama untuk seri tersebut.
data_lopo_nonsea.ts <- diff(data_lopo_sea.ts, differences = 1)par (mfrow = c(1, 1)) # 1 graph in a single device
plot (data_lopo_nonsea.ts, type ="o", pch = 16,
ylab = "DATa No.1",
main = "Non-Seasonal Differencing",
panel.first = grid ())
par(mfrow = c(2, 1)) # 2 graphs in a single device
acf(data_lopo_nonsea.ts, 36,
main = "ACF for DATA No.1 HIV (Non-Seasonal Diff)",
panel.first = grid (), xaxt = "n")
axis(1, at = 0:36/12, labels=0:36)
pacf(data_lopo_nonsea.ts, 36, main = "PACF for DATA No. HIV (Non-Seasonal Diff)", panel.first = grid (), xaxt = "n")
axis(1, at = 0:36/12, labels=0:36)
Penjelasan: Setelah melakukan non-seasonal differencing, time plot tidak
menunjukkan adanya trend, baik ACF maupun PACF menunjukkan beberapa
lonjakan yang signifikan tanpa adanya pola yang jelas.
Ini dapat menyimpulkan bahwa akhirnya, deret kami sekarang stasioner setelah pembedaan musiman dan nonmusiman pertama, sehingga D=1 dan d=1. Untuk menentukan urutan AR(p), MA(q), SAR(P), SMA(Q), akan diamati lonjakan pada periode lag tertentu.
untuk mengidentifikasi bagian musiman, salah satu kebutuhan untuk mengamati lonjakan pada lag 12 atau 24 atau 36 (untuk data bulanan) dan untuk mengidentifikasi bagian non-musiman adalah dengan mengamati lonjakan signifikan pada periode selain itu.
Berdasarkan ACF, dua lonjakan signifikan terdeteksi pada lag 1 dan lag 8. (Perhatikan bahwa, kami mengecualikan ACF pada lag = 0 sebagai bagian dari perhitungan lonjakan signifikan) Kedua lonjakan ini dapat digunakan untuk menentukan bagian MA non-musiman. Lonjakan signifikan lainnya diamati pada lag 12 untuk menyarankan bagian SMA musiman dari model. Dengan demikian, MA (2) dan SMA (1) disarankan berdasarkan plot ACF.
Selanjutnya, untuk mengidentifikasi urutan p dan P, kami mengamati PACF. Ada dua kemungkinan lonjakan signifikan yaitu pada lag 1 dan lag 8, menunjukkan AR(2) untuk AR non-musiman dan ada lonjakan signifikan pada lag 12, menunjukkan SAR(1) untuk AR musiman.
selanjutnya untuk memastikan bahwa model yang ditentukan dengan baik tidak terlewatkan, beberapa formulasi model akan diidentifikasi dan diestimasi.
Model berikut telah diidentifikasi: SARIMA(2, 1, 2)(1, 1, 1), SARIMA(1, 1, 2)(0, 1, 1), SARIMA(1, 1, 2)(1, 1, 1) dan SARIMA(2, 1, 3)(1, 1, 1).
# Model 1: ARIMA (2, 1, 2)(1, 1, 1)_30
model1_arima <- arima (data_lopo.ts, order = c(2, 1, 2),
seas = list(order = c(1, 1, 1), 12))
summary (model1_arima)##
## Call:
## arima(x = data_lopo.ts, order = c(2, 1, 2), seasonal = list(order = c(1, 1,
## 1), 12))
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1 sma1
## -0.2545 0.4932 0.5020 -0.4980 -0.2797 -1.0000
## s.e. 0.1825 0.1489 0.1871 0.1867 0.0531 0.0268
##
## sigma^2 estimated as 3239571: log likelihood = -3121.06, aic = 6256.12
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -26.22204 1767.084 876.5016 -5.934164 24.26688 0.969964
## ACF1
## Training set -0.04169811# Model 2: SARIMA (1, 1, 2)(0, 1, 1)_12
model2_sarima <- arima (data_lopo.ts, order = c(1, 1, 2),
seas = list(order = c(0, 1, 1), 12))
summary (model2_sarima)##
## Call:
## arima(x = data_lopo.ts, order = c(1, 1, 2), seasonal = list(order = c(0, 1,
## 1), 12))
##
## Coefficients:
## ar1 ma1 ma2 sma1
## 0.5515 -0.3917 0.0253 -0.9999
## s.e. 0.2476 0.2552 0.0654 0.0253
##
## sigma^2 estimated as 3812928: log likelihood = -3141.97, aic = 6293.94
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -19.69458 1917.091 821.9124 -4.405386 21.01595 0.9095539
## ACF1
## Training set -0.001980122# Model 3: ARIMA (1, 1, 1)
model3 <- arima(data_lopo.ts, order = c(1, 1, 1))
summary(model3)##
## Call:
## arima(x = data_lopo.ts, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## 0.6146 -0.5030
## s.e. 0.1815 0.1965
##
## sigma^2 estimated as 4318858: log likelihood = -3251.9, aic = 6509.81
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 24.47842 2075.298 839.7152 -5.023481 23.54406 0.929255 -0.01517494# Model 4: ARIMA (1, 1, 0)
model4 <- arima(data_lopo.ts, order = c(1, 1, 0))
summary(model4)##
## Call:
## arima(x = data_lopo.ts, order = c(1, 1, 0))
##
## Coefficients:
## ar1
## 0.1094
## s.e. 0.0524
##
## sigma^2 estimated as 4353543: log likelihood = -3253.33, aic = 6510.66
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 27.95523 2083.615 866.488 -7.516636 24.35049 0.9588826 -0.01171052Penjelasan: Berdasarkan nilai AIC terendah, kami memilih model 4, ARIMA (1, 1, 0) sebagai model terbaik. Selain itu, berdasarkan uji Ljung-Box, kami menerima hipotesis nol yang menunjukkan bahwa residunya adalah white noise.
Box.test(model4$residuals, type = "Ljung-Box")##
## Box-Ljung test
##
## data: model4$residuals
## X-squared = 0.049782, df = 1, p-value = 0.8234Selanjutnya dilanjutkan dengan penerapan ARIMA (1, 1, 0) untuk menghasilkan nilai forecast 3 tahun ke depan (2020-2023).
Stage 3: Model Application
# forecasting CPI data by using ARIMA (1, 1, 0) model
model4.for <- predict(model4, n.ahead = 3)
ts.plot(data_lopo.ts, model4.for$pred, ylim = c(1500, 20000),
ylab = "CPI", col = "blue")
grid ()
# calculating upper and lower limit
# additional plotting
U = model4.for$pred + model4.for$se
L = model4.for$pred - model4.for$se
xx = c(time (U), rev (time (U)))
yy = c(L, rev(U))
polygon(xx, yy, border = 8, col = gray (0.6, alpha=0.2))
lines(model4.for$pred, type = "p", col = "red")
model4.for## $pred
## Jan Feb Mar
## 2020 14038.70 14039.98 14040.12
##
## $se
## Jan Feb Mar
## 2020 2086.515 3116.349 3896.913# plotting forecast value for 2020-2023
plot(model4.for$pred, type = "o", ylim = c(1500, 20000), ylab ="CPI",
main = "Forecast Value for 1990-2024")
grid()
polygon(xx, yy, border = 8, col = gray (0.6, alpha=0.2))
text(model4.for$pred~time(model4.for$pred), labels = round (pred, 2),
data=model4.for, cex=0.7, font=2, pos=1)
Explanation: The forecast values for 2020-2023 are 112.39, 114.52, 116.34, 117.91 and 119.25.