Penerapan Model ARIMA
Rachmat Darmawan
3/21/2021
Penerapan pada Data Aktual
5.11 Tabel 8.2 berisi data penjualan produk farmasi. Diketahui datanya sebagai berikut :
#Membaca dan membuat plot data time series
library("readxl")## Warning: package 'readxl' was built under R version 4.0.3data1 <- read_xlsx("D:/Penjualan Produk Farmasi.xlsx")
data1 <- ts(data1)
plot.ts(data1, xlab="Periode Waktu",ylab="Penjualan", main="Time Series Plot Data Aktual n=120", cex.main=0.9)
points(data1)
Rangkaian proses pengolahan data time series ini menggunakan beberapa paket yaitu :
library("TTR")## Warning: package 'TTR' was built under R version 4.0.3library("TSA")## Warning: package 'TSA' was built under R version 4.0.3##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarlibrary("tseries")## Warning: package 'tseries' was built under R version 4.0.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary("graphics")
library("forecast")## Warning: package 'forecast' was built under R version 4.0.3## Registered S3 methods overwritten by 'forecast':
## method from
## fitted.Arima TSA
## plot.Arima TSA(a): 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.
Terdapat pada data aktual. Untuk memodelkan ARIMA, data yang digunakan adalah sebanyak 110 pengamatan pertama, selebihnya sebanyak 10 pengamatan terakhir digunakan untuk peramalan.
#Membuang 10 observasi terakhir
dataa <- head(data1,n=110)
dataa <- ts(dataa)
plot.ts(dataa,xlab="Periode Waktu",ylab="Penjualan", main="Time Series Plot Data Aktual n=110", cex.main=0.9)
points(dataa)
Konsep stasioneritas secara praktis dapat dideskripsikan jika suatu data deret waktu di plot dan kemudian tidak terdapat perubahan nilai tengah dari waktu ke waktu (stasioner terhadap rata-rata) dan apabila plot deret waktu tidak memperlihatkan perubahan ragam yang jelas dari waktu ke waktu. Berdasarkan time series plot tersebut terlihat datanya sudah stasioner, artinya waktu pengamatan tidak berpengaruh terhadap rata-rata dan ragam dari data. Ini berarti bahwa data deret waktu tersebut berfluktuasi disekitar rata-rata dan ragam yang konstan, dan dapat dikatakan bahwa data deret waktu tersebut stasioner dalam rata-rata dan ragam.
Karena datanya sudah stasioner maka akan dimodelkan dengan pemodelan ARIMA(p,d,q) dengan nilai d=0 sehingga nantinya akan membentuk model ARIMA(p,0,q) atau ARMA(p,q). Penentuan orde ke-p pada AR dan orde ke-q pada MA didasarkan cut off pada lag ke-p plot PACF untuk model AR(p) dan cut off pada lag ke-q plot ACF untuk model MA(q). selain plot PACF dan ACF, identifikasi model ARMA(p,q) dapat melalui EACF.
#identification model with acf for MA(q)
acf(dataa, lag.max=20, main="ACF Plot for MA(q)")
#identification model with pacf for AR(p)
pacf(dataa, lag.max=20, main="PACF Plot for AR(p)", ylab="PACF")
Pada plot ACF terlihat bahwa data turun secara eksponensial dan terjadi cut off setelah lag ke-3 sehingga dapat disimpulkan diperoleh model MA(3). Plot PACF juga terlihat bahwa data turun secara eksponensial dan cut off setelah lag ke-3 sehingga dapat disimpulkan diperoleh model AR(3). Kandidat model ARIMA yang akan diidentifikasi adalah model ARIMA(3,0,0) dan ARIMA(0,0,3)
#identification model with eacf for ARMA(p,q)
eacf(dataa)## 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 oTerlihat pada EACF terdapat simbol “o” dan “x”. cara membaca modelnya adalah dengan melihat simbol “o” paling ujung kiri atas kemudian sesuaikan dengan orde MA(q) secara horizontal dan orde AR(p) secara vertikal maka akan terbentuk model ARMA(p,q). Sehingga dapat diketahui kandidat model yang terbentuk adalah model ARIMA(1,0,1) dan ARIMA(3,0,1).
Kandidat-kandidat model tersebut akan diidentifikasi dan dipilih model yang terbaik berdasarkan nilai AIC minimum.
#candidate for MA(q) model
model003 <- arima(dataa, order = c(0,0,3), method = "ML")
model003##
## Call:
## arima(x = dataa, order = c(0, 0, 3), 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#candidate for AR(p) model
model300 <- arima(dataa, order = c(3,0,0), method = "ML")
model300##
## Call:
## arima(x = dataa, order = c(3, 0, 0), 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#candidate for ARMA(p,q) model
model101 <- arima(dataa, order = c(1,0,1), method = "ML")
model101##
## Call:
## arima(x = dataa, order = c(1, 0, 1), method = "ML")
##
## Coefficients:
## ar1 ma1 intercept
## 0.0672 0.0246 10375.1957
## s.e. 0.4398 0.4314 22.7734
##
## sigma^2 estimated as 47350: log likelihood = -748.18, aic = 1502.36model301 <- arima(dataa, order = c(3,0,1), method = "ML")
model301##
## Call:
## arima(x = dataa, order = c(3, 0, 1), 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.28Berdasarkan keempat model tersebut, yang memiliki nilai AIC minimum adalah model ARIMA(3,0,1) sehaingga model yang terbentuk adalah ini selanjutnya akan digunakan dalam pemodelan data deret waktu untuk 110 data pertama dan 10 data terakhir digunakan untuk peramalan untuk 10 periode berikutnya.
#selected model with AIC minimum
model <- fitted(model301)
cbind(data1,model)## Time Series:
## Start = 1
## End = 120
## Frequency = 1
## Penjualan Produk
## [1,] 10618.1
## [2,] 10537.9
## [3,] 10209.3
## [4,] 10553.0
## [5,] 9934.9
## [6,] 10534.5
## [7,] 10196.5
## [8,] 10511.8
## [9,] 10089.6
## [10,] 10371.2
## [11,] 10239.4
## [12,] 10472.4
## [13,] 10827.2
## [14,] 10640.8
## [15,] 10517.8
## [16,] 10154.2
## [17,] 9969.2
## [18,] 10260.4
## [19,] 10737.0
## [20,] 10430.0
## [21,] 10689.0
## [22,] 10430.4
## [23,] 10002.4
## [24,] 10135.7
## [25,] 10096.2
## [26,] 10288.7
## [27,] 10289.1
## [28,] 10589.9
## [29,] 10551.9
## [30,] 10208.3
## [31,] 10334.5
## [32,] 10480.1
## [33,] 10387.6
## [34,] 10202.6
## [35,] 10219.3
## [36,] 10382.7
## [37,] 10820.5
## [38,] 10358.7
## [39,] 10494.6
## [40,] 10497.6
## [41,] 10431.5
## [42,] 10447.8
## [43,] 10684.4
## [44,] 10176.5
## [45,] 10616.0
## [46,] 10627.7
## [47,] 10684.0
## [48,] 10246.7
## [49,] 10265.0
## [50,] 10090.4
## [51,] 9881.1
## [52,] 10449.7
## [53,] 10276.3
## [54,] 10175.2
## [55,] 10212.5
## [56,] 10395.5
## [57,] 10545.9
## [58,] 10635.7
## [59,] 10265.2
## [60,] 10551.6
## [61,] 10538.2
## [62,] 10286.2
## [63,] 10171.3
## [64,] 10393.1
## [65,] 10162.3
## [66,] 10164.5
## [67,] 10327.0
## [68,] 10365.1
## [69,] 10755.9
## [70,] 10463.6
## [71,] 10080.5
## [72,] 10479.6
## [73,] 9980.9
## [74,] 10039.2
## [75,] 10246.1
## [76,] 10368.0
## [77,] 10446.3
## [78,] 10535.3
## [79,] 10786.9
## [80,] 9975.8
## [81,] 10160.9
## [82,] 10422.1
## [83,] 10757.2
## [84,] 10463.8
## [85,] 10307.0
## [86,] 10134.7
## [87,] 10207.7
## [88,] 10488.0
## [89,] 10262.3
## [90,] 10785.9
## [91,] 10375.4
## [92,] 10123.4
## [93,] 10462.7
## [94,] 10205.5
## [95,] 10522.7
## [96,] 10253.2
## [97,] 10428.7
## [98,] 10615.8
## [99,] 10417.3
## [100,] 10445.4
## [101,] 10690.6
## [102,] 10271.8
## [103,] 10524.8
## [104,] 9815.0
## [105,] 10398.5
## [106,] 10553.1
## [107,] 10655.8
## [108,] 10199.1
## [109,] 10416.6
## [110,] 10391.3
## [111,] 10210.1
## [112,] 10352.5
## [113,] 10423.8
## [114,] 10519.3
## [115,] 10596.7
## [116,] 10650.0
## [117,] 10741.6
## [118,] 10246.0
## [119,] 10354.4
## [120,] 10155.4plot.ts(data1, xlab="Periode Waktu",ylab="Penjualan", main="Time Series Plot Data with ARIMA (3,0,1) Model", cex.main=0.9, ylim=c(9800,11100),col="black")
points(data1, pch=20,col="black")
#add line ARIMA model
lines(model, col="red")
points(model, pch=18, col="red")
legend("topright", c("Actual","ARIMA (3,0,1)"), lty=1, col=c("black", "red"), cex=0.56, bty="n", pch=c(16,16))
Berikut ini adalah nilai-nilai dugaan model ARIMA(3,0,1)
cbind(data1, model)## Time Series:
## Start = 1
## End = 120
## Frequency = 1
## Penjualan Produk
## [1,] 10618.1
## [2,] 10537.9
## [3,] 10209.3
## [4,] 10553.0
## [5,] 9934.9
## [6,] 10534.5
## [7,] 10196.5
## [8,] 10511.8
## [9,] 10089.6
## [10,] 10371.2
## [11,] 10239.4
## [12,] 10472.4
## [13,] 10827.2
## [14,] 10640.8
## [15,] 10517.8
## [16,] 10154.2
## [17,] 9969.2
## [18,] 10260.4
## [19,] 10737.0
## [20,] 10430.0
## [21,] 10689.0
## [22,] 10430.4
## [23,] 10002.4
## [24,] 10135.7
## [25,] 10096.2
## [26,] 10288.7
## [27,] 10289.1
## [28,] 10589.9
## [29,] 10551.9
## [30,] 10208.3
## [31,] 10334.5
## [32,] 10480.1
## [33,] 10387.6
## [34,] 10202.6
## [35,] 10219.3
## [36,] 10382.7
## [37,] 10820.5
## [38,] 10358.7
## [39,] 10494.6
## [40,] 10497.6
## [41,] 10431.5
## [42,] 10447.8
## [43,] 10684.4
## [44,] 10176.5
## [45,] 10616.0
## [46,] 10627.7
## [47,] 10684.0
## [48,] 10246.7
## [49,] 10265.0
## [50,] 10090.4
## [51,] 9881.1
## [52,] 10449.7
## [53,] 10276.3
## [54,] 10175.2
## [55,] 10212.5
## [56,] 10395.5
## [57,] 10545.9
## [58,] 10635.7
## [59,] 10265.2
## [60,] 10551.6
## [61,] 10538.2
## [62,] 10286.2
## [63,] 10171.3
## [64,] 10393.1
## [65,] 10162.3
## [66,] 10164.5
## [67,] 10327.0
## [68,] 10365.1
## [69,] 10755.9
## [70,] 10463.6
## [71,] 10080.5
## [72,] 10479.6
## [73,] 9980.9
## [74,] 10039.2
## [75,] 10246.1
## [76,] 10368.0
## [77,] 10446.3
## [78,] 10535.3
## [79,] 10786.9
## [80,] 9975.8
## [81,] 10160.9
## [82,] 10422.1
## [83,] 10757.2
## [84,] 10463.8
## [85,] 10307.0
## [86,] 10134.7
## [87,] 10207.7
## [88,] 10488.0
## [89,] 10262.3
## [90,] 10785.9
## [91,] 10375.4
## [92,] 10123.4
## [93,] 10462.7
## [94,] 10205.5
## [95,] 10522.7
## [96,] 10253.2
## [97,] 10428.7
## [98,] 10615.8
## [99,] 10417.3
## [100,] 10445.4
## [101,] 10690.6
## [102,] 10271.8
## [103,] 10524.8
## [104,] 9815.0
## [105,] 10398.5
## [106,] 10553.1
## [107,] 10655.8
## [108,] 10199.1
## [109,] 10416.6
## [110,] 10391.3
## [111,] 10210.1
## [112,] 10352.5
## [113,] 10423.8
## [114,] 10519.3
## [115,] 10596.7
## [116,] 10650.0
## [117,] 10741.6
## [118,] 10246.0
## [119,] 10354.4
## [120,] 10155.4(b): Forecast the last 10 observation
Setelah mendapatkan model terbaik, selanjutnya dilakukan peramalan untuk 10 waktu kedepan.
ARIMA_data <- Arima(dataa, order = c(3,0,1), method = "ML")peramalan <- forecast(ARIMA_data)
hasilperamalan <- as.data.frame(peramalan)
hasilperamalan## 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.29plot(peramalan)
Garis biru mewakili prediksi titik, area berbayang biru tua menunjukkan interval prediksi 80% dan area berbayang biru muda menunjukkan interval prediksi 95% untuk prediksi titik
Penilaian Akurasi dengan MAPE
Penilaian akurasi dilakukan dengan membandingkan data hasil ramalan dengan data aktual 10 amatan terakhir.
data_forecast1 <- hasilperamalan$`Point Forecast`
data_forecast1 <- as.data.frame(data_forecast1)
head(data_forecast1)## data_forecast1
## 1 10371.02
## 2 10361.97
## 3 10362.54
## 4 10367.76
## 5 10373.12
## 6 10376.24data3 <- data1[-111:-120,]
data2 <- data1[-1:-110,]
data_ <- as.data.frame(data2)
data_## data2
## 1 10210.1
## 2 10352.5
## 3 10423.8
## 4 10519.3
## 5 10596.7
## 6 10650.0
## 7 10741.6
## 8 10246.0
## 9 10354.4
## 10 10155.4akurasi.arima <- accuracy(data_forecast1$data_forecast1,data_$`data2`)
akurasi.arima## ME RMSE MAE MPE MAPE
## Test set 53.57804 193.8149 161.3917 0.4821229 1.538676Jika nilai MAPE kecil maka kesalahan hasil pendugaannya juga kecil. Peramalan penjualan menggunakan metode ARMA dinilai memiliki kemampuan ramalan yang sangat akurat karena hasil dari MAPE adalah sebesar 1.538676%
(c): In Exercise 4.12, you were asked to use simple exponential smoothing with λ=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?
data.ts<-ts(data3)
plot(data.ts, main="Plot Data Aktual Time Series")
points(data.ts)
# Pemulusan Eksponensial Sederhana
data1.ses<-HoltWinters(data.ts,alpha=0.1, beta=FALSE, gamma=FALSE)
data1.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 10393.24SSE <- data1.ses$SSE
SSE## [1] 5954375MSE <- SSE/110
MSE## [1] 54130.69hasil_forecast.sess <- forecast(data1.ses)
hasil_forecast.ses <- as.data.frame(hasil_forecast.sess)
hasil_forecast.ses## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 111 10393.24 10093.50 10692.98 9934.830 10851.65
## 112 10393.24 10092.01 10694.48 9932.543 10853.94
## 113 10393.24 10090.52 10695.96 9930.268 10856.22
## 114 10393.24 10089.04 10697.44 9928.004 10858.48
## 115 10393.24 10087.57 10698.92 9925.751 10860.73
## 116 10393.24 10086.10 10700.38 9923.509 10862.97
## 117 10393.24 10084.64 10701.84 9921.278 10865.21
## 118 10393.24 10083.19 10703.29 9919.057 10867.43
## 119 10393.24 10081.74 10704.74 9916.846 10869.64
## 120 10393.24 10080.30 10706.18 9914.645 10871.84plot(hasil_forecast.sess)
garis biru mewakili prediksi titik, area berbayang biru tua menunjukkan interval prediksi 80% dan area berbayang biru muda menunjukkan interval prediksi 95% untuk prediksi titik
data_forecast2 <- hasil_forecast.ses$`Point Forecast`
data_forecast2 <- as.data.frame(data_forecast2)
head(data_forecast2)## data_forecast2
## 1 10393.24
## 2 10393.24
## 3 10393.24
## 4 10393.24
## 5 10393.24
## 6 10393.24dataaktual <- as.data.frame(data2)
dataaktual## data2
## 1 10210.1
## 2 10352.5
## 3 10423.8
## 4 10519.3
## 5 10596.7
## 6 10650.0
## 7 10741.6
## 8 10246.0
## 9 10354.4
## 10 10155.4akurasi.ses <- accuracy(data_forecast2$data_forecast2,dataaktual$`data2`)
akurasi.ses## ME RMSE MAE MPE MAPE
## Test set 31.73833 189.8569 161.3 0.2723997 1.540696Jika nilai MAPE kecil maka kesalahan hasil pendugaannya juga kecil. Peramalan penjualan menggunakan metode ARMA dinilai memiliki kemampuan ramalan yang sangat akurat karena hasil dari MAPE adalah sebesar 1.540696%
Perbandingan Model ARIMA dan Pemulusan Eksponensial Sederhana.
Setelah dilakukan peramalan dengan menggunakan dua metode, akan diamati nilai MAPE dari masing-masing metode. Metode dengan nilai MAPE paling kecil akan terpilih menjadi model terbaik.
akurasi.arima## ME RMSE MAE MPE MAPE
## Test set 53.57804 193.8149 161.3917 0.4821229 1.538676akurasi.ses## ME RMSE MAE MPE MAPE
## Test set 31.73833 189.8569 161.3 0.2723997 1.540696Hasil dari kedua metode tersebut menghasilkan nilai MAPE 1.538676% untuk metode ARIMA dan 1.540696% untuk metode Simple Exponential Smoothing. Dengan demikian, metode yang mempunyai MAPE lebih kecil adalah metode ARIMA(3,0,1) dengan nilai \(p=3,d=0,q=1\)
(d): How would prediction intervals be obtained for the ARIMA forecasts?
fit <- Arima(data.ts, order=c(3,0,1), method="ML")
fit## Series: data.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.49hasilperamalan## 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.29# Interval Prediksi ARIMA
interval.prediction <- hasilperamalan[,-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.29plot(hasilperamalan$`Point Forecast`, type="n", ylim=range(hasilperamalan$`Lo 80`,hasilperamalan$`Hi 80`))
polygon(c(time(hasilperamalan$`Point Forecast`),rev(time(hasilperamalan$`Point Forecast`))), c(hasilperamalan$`Hi 80`,rev(hasilperamalan$`Lo 80`)),
col=rgb(0,0,0.6,0.2), border=FALSE)
lines(hasilperamalan$`Point Forecast`)
lines(fitted(fit),col='red')
out <- (hasilperamalan$`Point Forecast`< hasilperamalan$`Lo 80` | hasilperamalan$`Point Forecast` > hasilperamalan$`Hi 80`)plot(hasilperamalan$`Point Forecast`, type="n", ylim=range(hasilperamalan$`Lo 95`,hasilperamalan$`Hi 95`))
polygon(c(time(hasilperamalan$`Point Forecast`),rev(time(hasilperamalan$`Point Forecast`))), c(hasilperamalan$`Hi 95`,rev(hasilperamalan$`Lo 95`)),
col=rgb(0,0,0.6,0.2), border=FALSE)
lines(hasilperamalan$`Point Forecast`)
lines(fitted(fit),col='red')
out <- (hasilperamalan$`Point Forecast`< hasilperamalan$`Lo 95` | hasilperamalan$`Point Forecast` > hasilperamalan$`Hi 95`)Terlihat bahwa kedua plot interval prediksi 80% dan plot interval prediksi 95% menunjukkan hasil yang hampir sama dan hanya berbeda rentang selangnya. Interval prediksi 95% memiliki rentang yang lebih luas dari pada interval prediksi 80%.