Eksplorasi dan Pemulusan Data Time Series
Kelompok 5
P2-MPDW 2022
Anggota:
Hanung Safrizal (G1401201050)
Maulana Ahsan Fadillah (G1401201062)
Muhammad Irsyad Robbani (G1401201064)
Muhammad Raziv Zulfikar (G1401201083)
Dhea Puspita Adinda (G1401201090)Pendahuluan
Latar Belakang
Beras merupakan makanan pokok bagi sebagian besar masyarakat Indonesia. Konsumsi beras di Indonesia meningkat setiap tahunnya seiring dengan meningkatnya jumlah penduduk di Indonesia. Ketergantungan masyarakat Indonesia terhadap beras bisa menimbulkan masalah jika ketersediaan beras sudah tidak dapat tercukupi. Pemenuhan kebutuhan dan stabilitas harga beras merupakan isu yang tetap relevan dari waktu ke waktu. Jumlah penduduk yang besar menjadikan permasalahan beras bukan hanya pada jumlah ketersediaan dan harga saja, namun termasuk keberagaman kualitas dan jenis produk beras (BPS 2009). Pada awal tahun 2018 harga beras mengalami peningkatan. Kenaikan harga beras ini jika terus dibiarkan akan menyebabkan terjadinya inflasi yang berdampak pada melambatnya pertumbuhan ekonomi nasional serta dampak negatif lainnya. Dalam rangka perumusan kebijakan pengendalian inflasi maka data dan informasi terkait proyeksi keadaan pasar sangat dibutuhkan. Oleh karena itu, pemodelan harga beras di Indonesia sangat perlu dilakukan.
Tujuan
- Mengeksplorasi dan forecasting data time series harga beras premium di Indonesia per bulan pada tahun 2013-2022.
Tinjauan Pustaka
Metode Single Moving Average
Metode single moving average adalah metode peramalan yang menggunakan sejumlah data aktual permintaan yang baru untuk membangkitkan nilai ramalan untuk permintaan dimasa yang akan datang.(Naufal dan Andrean, 2017).
Metode Double Moving Average
Metode Double Moving Average adalah merupakan suatu metode peramalan yang dilakukan dengan menghitung nilai rata-rata bergerak sebanyak dua kali, lalu dilanjutkan dengan meramal menggunakan suatu persamaan. Metode DMA digunakan untuk mengatasi galat sistematis yang terjadi apabila rata-rata bergerak data dipakai pada data yang memiliki tren. Metode ini merupakan pengembangan dari metode SMA (Muzaki 2022).
Metode Single Exponential Smoothing
Metode Single Exponential Smoothing adalah metode yang menunjukan pembobotan menurun secara eksponensial terhadap nilai observasi yang lebih tua. Yaitu nilai yang lebih baru diberikan bobot yang relatif lebih besar dibanding nilai observasi yang lebih lama. Metode ini memberikan sebuah pembobotan eksponensial rata-rata bergerak dari semua nilai observasi sebelumnya (Hartono et.al, 2012).
Metode Double Exponential Smoothing
Metode Double Exponential Smoothing merupakan metode yang digunakan untuk meramalkan data yang mengalami trend kenaikan. Jika data yang digunakan semakin banyak dalam perhitungan peramalan, maka persentase error peramalannya akan semakin kecil dan sebaliknya (Hanief dan Purwanto 2017).
Winter Smoothing
Metode pemulusan eksponensial linear dari Winter’s digunakan untuk peramalan jika data memiliki komponen musiman. Metode winter didasarkan pada tiga persamaan pemulusan, yaitu stasioneritas, trend, dan musiman. Holt Winter sendiri memiliki dua metode yang berbeda, bergantung pada sifat musiman itu sendiri, yaitu aditif atau multiplikatif.
SSE
Sum Square Error (SSE) menjumlahkan kuadrat dari kesalahan.
MSE
Mean Squared Error (MSE) merata-ratakan kuadrat dari kesalahan menggunakan penyebut n tanpa memperhatikan derajat bebas model. MSE banyak dipakai dalam optimalisasi pembobotan karena memberikan ketelitian yang lebih baik.
RMSE
Root Mean Square Error (RMSE) diperoleh dengan mengakarkan hasil MSE.
MAD
Simpangan absolut rata-rata atau mean absolute deviation (MAD) merata-ratakan nilai absolut kesalahan peramalan dalam unit ukuran yang sama seperti data aktual.
MAPE
Mean absolute percentage error (MAPE) merata-ratakan persentase absolut kesalahan dan memberikan gambaran seberapa kesalahan peramalan dibanding dengan nilai sebenarnya. MAPE banyak dipakai dalam perbandingan data-data yang mempunyai skala interval waktu yang berbeda.
Package dan Data
1. Memanggil Library
Packages yang digunakan dalam analisis eksplorasi dan pemulusan data deret waktu, sebagai berikut:
library("tseries")## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary("forecast")
library("TTR")
library("TSA")## Registered S3 methods overwritten by 'TSA':
## method from
## fitted.Arima forecast
## plot.Arima forecast##
## Attaching package: 'TSA'## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarlibrary("graphics")
library("astsa")##
## Attaching package: 'astsa'## The following object is masked from 'package:forecast':
##
## gaslibrary("car")## Loading required package: carDatalibrary("readxl")
library("dplyr")##
## Attaching package: 'dplyr'## The following object is masked from 'package:car':
##
## recode## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union2. Memanggil Data
Data yang digunakan pada pengerjaan tugas analisis eksplorasi dan pemulusan deret waktua adalah data beras premium di Indonesia per bulan dari tahun 2013 hingga 2022. Data diperoleh dari web BPS Indonesia ( https://www.bps.go.id/indicator/36/500/1/harga-beras-di-penggilingan-menurut-kualitas.html). Pada data harga beras terdapat 3 peubah yang digunakan yaitu beras premium, medium, dan luar kualitas tetapi peubah yang digunakan dalam eksplorasi kali ini adalah data beras premium dengan jumlah amatan sebanyak 115 amatan.
data <- read_xlsx("C:/0 SEM5/MPDW/p1/Data Kelompok 15.xlsx")
data$Premium <- as.numeric(data$Premium)3. Split Data
Spliting data dilakukan untuk memisahkan data menjadi dua bagian yaitu data training dan data testing. Data training digunakan untuk menentukan model yang sesuai dengan keseluruhan data dan membuat model, sedangkan data testing digunakan untuk menguji performa dari model yang telah didapatkan. Jumlah data training sebanyak 88 amatan (Januari 2013 hingga April 2020) dan data testing sebanyak 27 amatan (September 2020 hingga Juli 2022).
train <- data[1:88, 3]
test <- data[89:115, 3]
data.premium <- ts(data[,3])
train.ts <- ts(train)
test.ts <- ts(test)Eksplorasi Data Analisis
Plot Data Time Series
Plot Time Series Harga Beras Premium
data.premium <- ts(data$Premium, start = 2013, frequency = 12)
ts.plot(data.premium, main = "Harga Beras Premium Tahun 2013 - 2022", ylab = "Harga", lwd = 1.5)
points(data.premium)
legend("topleft", c("Premium"), cex = 0.7,
col = c("black", "red", "blue"), lty = 1)
Single Moving Average (SMA)
Pemulusan SMA dengan n=4
data.sma<-SMA(train.ts, n=4)
data.sma## Time Series:
## Start = 1
## End = 88
## Frequency = 1
## [1] NA NA NA 7641.970 7578.912 7522.652 7584.505
## [8] 7669.647 7719.840 7794.297 7818.372 7872.230 7987.885 8102.180
## [15] 8170.402 8156.205 8106.365 8072.450 8081.307 8183.692 8258.962
## [22] 8316.285 8397.995 8570.225 8923.462 9270.422 9524.597 9551.942
## [29] 9298.530 9081.195 8924.512 8924.265 9107.827 9242.740 9397.527
## [36] 9531.725 9601.575 9683.997 9685.892 9551.942 9416.560 9308.970
## [43] 9259.602 9319.485 9301.647 9246.207 9216.955 9210.525 9290.722
## [50] 9359.675 9392.467 9388.222 9389.427 9398.400 9397.187 9425.222
## [57] 9433.822 9448.420 9487.187 9593.100 9812.930 10032.697 10121.150
## [64] 10037.292 9830.785 9604.862 9511.705 9494.982 9507.020 9548.832
## [71] 9611.592 9701.592 9836.355 9927.007 9937.897 9849.550 9687.292
## [78] 9564.192 9490.397 9506.727 9539.812 9575.727 9631.430 9708.495
## [85] 9818.062 9923.470 10008.537 10053.520data.ramal<-c(NA,data.sma)
data.ramal #forecast 1 periode ke depan## [1] NA NA NA NA 7641.970 7578.912 7522.652
## [8] 7584.505 7669.647 7719.840 7794.297 7818.372 7872.230 7987.885
## [15] 8102.180 8170.402 8156.205 8106.365 8072.450 8081.307 8183.692
## [22] 8258.962 8316.285 8397.995 8570.225 8923.462 9270.422 9524.597
## [29] 9551.942 9298.530 9081.195 8924.512 8924.265 9107.827 9242.740
## [36] 9397.527 9531.725 9601.575 9683.997 9685.892 9551.942 9416.560
## [43] 9308.970 9259.602 9319.485 9301.647 9246.207 9216.955 9210.525
## [50] 9290.722 9359.675 9392.467 9388.222 9389.427 9398.400 9397.187
## [57] 9425.222 9433.822 9448.420 9487.187 9593.100 9812.930 10032.697
## [64] 10121.150 10037.292 9830.785 9604.862 9511.705 9494.982 9507.020
## [71] 9548.832 9611.592 9701.592 9836.355 9927.007 9937.897 9849.550
## [78] 9687.292 9564.192 9490.397 9506.727 9539.812 9575.727 9631.430
## [85] 9708.495 9818.062 9923.470 10008.537 10053.520data.gab<-cbind(aktual=c(train.ts,rep(NA,5)),pemulusan=c(data.sma,rep(NA,5)),ramalan=c(data.ramal,rep(data.ramal[length(data.ramal)],4)))
data.gab #forecast 5 periode ke depan## aktual pemulusan ramalan
## [1,] 7797.63 NA NA
## [2,] 7773.26 NA NA
## [3,] 7576.27 NA NA
## [4,] 7420.72 7641.970 NA
## [5,] 7545.40 7578.912 7641.970
## [6,] 7548.22 7522.652 7578.912
## [7,] 7823.68 7584.505 7522.652
## [8,] 7761.29 7669.647 7584.505
## [9,] 7746.17 7719.840 7669.647
## [10,] 7846.05 7794.297 7719.840
## [11,] 7919.98 7818.372 7794.297
## [12,] 7976.72 7872.230 7818.372
## [13,] 8208.79 7987.885 7872.230
## [14,] 8303.23 8102.180 7987.885
## [15,] 8192.87 8170.402 8102.180
## [16,] 7919.93 8156.205 8170.402
## [17,] 8009.43 8106.365 8156.205
## [18,] 8167.57 8072.450 8106.365
## [19,] 8228.30 8081.307 8072.450
## [20,] 8329.47 8183.692 8081.307
## [21,] 8310.51 8258.962 8183.692
## [22,] 8396.86 8316.285 8258.962
## [23,] 8555.14 8397.995 8316.285
## [24,] 9018.39 8570.225 8397.995
## [25,] 9723.46 8923.462 8570.225
## [26,] 9784.70 9270.422 8923.462
## [27,] 9571.84 9524.597 9270.422
## [28,] 9127.77 9551.942 9524.597
## [29,] 8709.81 9298.530 9551.942
## [30,] 8915.36 9081.195 9298.530
## [31,] 8945.11 8924.512 9081.195
## [32,] 9126.78 8924.265 8924.512
## [33,] 9444.06 9107.827 8924.265
## [34,] 9455.01 9242.740 9107.827
## [35,] 9564.26 9397.527 9242.740
## [36,] 9663.57 9531.725 9397.527
## [37,] 9723.46 9601.575 9531.725
## [38,] 9784.70 9683.997 9601.575
## [39,] 9571.84 9685.892 9683.997
## [40,] 9127.77 9551.942 9685.892
## [41,] 9181.93 9416.560 9551.942
## [42,] 9354.34 9308.970 9416.560
## [43,] 9374.37 9259.602 9308.970
## [44,] 9367.30 9319.485 9259.602
## [45,] 9110.58 9301.647 9319.485
## [46,] 9132.58 9246.207 9301.647
## [47,] 9257.36 9216.955 9246.207
## [48,] 9341.58 9210.525 9216.955
## [49,] 9431.37 9290.722 9210.525
## [50,] 9408.39 9359.675 9290.722
## [51,] 9388.53 9392.467 9359.675
## [52,] 9324.60 9388.222 9392.467
## [53,] 9436.19 9389.427 9388.222
## [54,] 9444.28 9398.400 9389.427
## [55,] 9383.68 9397.187 9398.400
## [56,] 9436.74 9425.222 9397.187
## [57,] 9470.59 9433.822 9425.222
## [58,] 9502.67 9448.420 9433.822
## [59,] 9538.75 9487.187 9448.420
## [60,] 9860.39 9593.100 9487.187
## [61,] 10349.91 9812.930 9593.100
## [62,] 10381.74 10032.697 9812.930
## [63,] 9892.56 10121.150 10032.697
## [64,] 9524.96 10037.292 10121.150
## [65,] 9523.88 9830.785 10037.292
## [66,] 9478.05 9604.862 9830.785
## [67,] 9519.93 9511.705 9604.862
## [68,] 9458.07 9494.982 9511.705
## [69,] 9572.03 9507.020 9494.982
## [70,] 9645.30 9548.832 9507.020
## [71,] 9770.97 9611.592 9548.832
## [72,] 9818.07 9701.592 9611.592
## [73,] 10111.08 9836.355 9701.592
## [74,] 10007.91 9927.007 9836.355
## [75,] 9814.53 9937.897 9927.007
## [76,] 9464.68 9849.550 9937.897
## [77,] 9462.05 9687.292 9849.550
## [78,] 9515.51 9564.192 9687.292
## [79,] 9519.35 9490.397 9564.192
## [80,] 9530.00 9506.727 9490.397
## [81,] 9594.39 9539.812 9506.727
## [82,] 9659.17 9575.727 9539.812
## [83,] 9742.16 9631.430 9575.727
## [84,] 9838.26 9708.495 9631.430
## [85,] 10032.66 9818.062 9708.495
## [86,] 10080.80 9923.470 9818.062
## [87,] 10082.43 10008.537 9923.470
## [88,] 10018.19 10053.520 10008.537
## [89,] NA NA 10053.520
## [90,] NA NA 10053.520
## [91,] NA NA 10053.520
## [92,] NA NA 10053.520
## [93,] NA NA 10053.520Hasil peramalan dengan menggunakan N=5 pada metode SMA adalah sebesar 10053.520
Plot time series
ts.plot(data.gab[,1], xlab="Time Period ", ylab="Sales", main= "SMA N=4 Data Premium")
points(data.gab[,1])
lines(data.gab[,2],col="green",lwd=2)
lines(data.gab[,3],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"), lty=8, col=c("black","green","red"), cex=0.8)
Menghitung nilai keakuratan
data.premium=as.numeric(data.premium)error.sma = test.ts-data.ramal[1:length(test.ts)]
SSE.sma = sum(error.sma[5:length(test.ts)]^2)
MSE.sma = mean(error.sma[5:length(test.ts)]^2)
MAPE.sma = mean(abs((error.sma[5:length(test.ts)]/test.ts[5:length(test.ts)])*100))
akurasi.sma <- matrix(c(SSE.sma, MSE.sma, MAPE.sma))
row.names(akurasi.sma)<- c("SSE", "MSE", "MAPE")
colnames(akurasi.sma) <- c("Akurasi m = 4")
akurasi.sma## Akurasi m = 4
## SSE 6.172562e+07
## MSE 2.683723e+06
## MAPE 1.616197e+01Double Moving Average (DMA)
Pemulusan DMA dengan n=4
dma <- SMA(data.sma, n = 4)
At <- 2*data.sma - dma
Bt <- 2/(4-1)*(data.sma - dma)
data.dma<- At+Bt
data.ramal2<- c(NA, data.dma)
t = 1:5
f = c()
for (i in t) {
f[i] = At[length(At)] + Bt[length(Bt)]*(i)
}
data.gab2 <- cbind(aktual = c(train.ts,rep(NA,5)), pemulusan1 = c(data.sma,rep(NA,5)),pemulusan2 = c(data.dma, rep(NA,5)),At = c(At, rep(NA,5)), Bt = c(Bt,rep(NA,5)),ramalan = c(data.ramal2, f[-1]))
data.gab2## aktual pemulusan1 pemulusan2 At Bt ramalan
## [1,] 7797.63 NA NA NA NA NA
## [2,] 7773.26 NA NA NA NA NA
## [3,] 7576.27 NA NA NA NA NA
## [4,] 7420.72 7641.970 NA NA NA NA
## [5,] 7545.40 7578.912 NA NA NA NA
## [6,] 7548.22 7522.652 NA NA NA NA
## [7,] 7823.68 7584.505 7588.663 7587.000 1.663333 NA
## [8,] 7761.29 7669.647 7804.178 7750.366 53.812083 7588.663
## [9,] 7746.17 7719.840 7879.305 7815.519 63.785833 7804.178
## [10,] 7846.05 7794.297 7964.672 7896.522 68.150000 7879.305
## [11,] 7919.98 7818.372 7931.428 7886.206 45.222083 7964.672
## [12,] 7976.72 7872.230 7990.638 7943.275 47.363333 7931.428
## [13,] 8208.79 7987.885 8187.366 8107.574 79.792500 7990.638
## [14,] 8303.23 8102.180 8363.869 8259.193 104.675417 8187.366
## [15,] 8192.87 8170.402 8399.116 8307.631 91.485417 8363.869
## [16,] 7919.93 8156.205 8242.933 8208.242 34.691250 8399.116
## [17,] 8009.43 8106.365 8060.660 8078.942 -18.282083 8242.933
## [18,] 8167.57 8072.450 7982.607 8018.544 -35.937083 8060.660
## [19,] 8228.30 8081.307 8043.350 8058.533 -15.182917 7982.607
## [20,] 8329.47 8183.692 8304.924 8256.431 48.492500 8043.350
## [21,] 8310.51 8258.962 8442.061 8368.822 73.239583 8304.924
## [22,] 8396.86 8316.285 8493.324 8422.508 70.815417 8442.061
## [23,] 8555.14 8397.995 8579.264 8506.756 72.507500 8493.324
## [24,] 9018.39 8570.225 8877.489 8754.583 122.905417 8579.264
## [25,] 9723.46 8923.462 9542.580 9294.933 247.647083 8877.489
## [26,] 9784.70 9270.422 10070.250 9750.319 319.930833 9542.580
## [27,] 9571.84 9524.597 10278.632 9977.018 301.613750 10070.250
## [28,] 9127.77 9551.942 9942.503 9786.279 156.224167 10278.632
## [29,] 8709.81 9298.530 9110.458 9185.687 -75.228750 9942.503
## [30,] 8915.36 9081.195 8609.743 8798.324 -188.580833 9110.458
## [31,] 8945.11 8924.512 8441.958 8634.980 -193.021667 8609.743
## [32,] 9126.78 8924.265 8702.831 8791.404 -88.573750 8441.958
## [33,] 9444.06 9107.827 9271.790 9206.205 65.585000 8702.831
## [34,] 9455.01 9242.740 9564.246 9435.644 128.602500 9271.790
## [35,] 9564.26 9397.527 9779.923 9626.965 152.958333 9564.246
## [36,] 9663.57 9531.725 9884.675 9743.495 141.180000 9779.923
## [37,] 9723.46 9601.575 9865.214 9759.758 105.455417 9884.675
## [38,] 9784.70 9683.997 9901.150 9814.289 86.860833 9865.214
## [39,] 9571.84 9685.892 9786.051 9745.987 40.063333 9901.150
## [40,] 9127.77 9551.942 9420.427 9473.033 -52.606250 9786.051
## [41,] 9181.93 9416.560 9136.496 9248.522 -112.025417 9420.427
## [42,] 9354.34 9308.970 9005.851 9127.099 -121.247500 9136.496
## [43,] 9374.37 9259.602 9051.825 9134.936 -83.110833 9005.851
## [44,] 9367.30 9319.485 9308.369 9312.816 -4.446250 9051.825
## [45,] 9110.58 9301.647 9308.683 9305.869 2.814167 9308.369
## [46,] 9132.58 9246.207 9186.994 9210.679 -23.685417 9308.683
## [47,] 9257.36 9216.955 9126.757 9162.836 -36.079167 9186.994
## [48,] 9341.58 9210.525 9155.010 9177.216 -22.205833 9126.757
## [49,] 9431.37 9290.722 9373.422 9340.342 33.080000 9155.010
## [50,] 9408.39 9359.675 9510.018 9449.881 60.137083 9373.422
## [51,] 9388.53 9392.467 9524.334 9471.587 52.746667 9510.018
## [52,] 9324.60 9388.222 9438.974 9418.673 20.300417 9524.334
## [53,] 9436.19 9389.427 9401.060 9396.407 4.652917 9438.974
## [54,] 9444.28 9398.400 9408.851 9404.671 4.180417 9401.060
## [55,] 9383.68 9397.187 9403.651 9401.066 2.585417 9408.851
## [56,] 9436.74 9425.222 9462.994 9447.886 15.108750 9403.651
## [57,] 9470.59 9433.822 9467.430 9453.987 13.442917 9462.994
## [58,] 9502.67 9448.420 9485.515 9470.677 14.837917 9467.430
## [59,] 9538.75 9487.187 9551.395 9525.712 25.682917 9485.515
## [60,] 9860.39 9593.100 9763.879 9695.567 68.311667 9551.395
## [61,] 10349.91 9812.930 10192.131 10040.451 151.680417 9763.879
## [62,] 10381.74 10032.697 10534.729 10333.916 200.812500 10192.131
## [63,] 9892.56 10121.150 10506.451 10352.331 154.120417 10534.729
## [64,] 9524.96 10037.292 10097.751 10073.567 24.183333 10506.451
## [65,] 9523.88 9830.785 9539.625 9656.089 -116.464167 10097.751
## [66,] 9478.05 9604.862 9115.429 9311.203 -195.773333 9539.625
## [67,] 9519.93 9511.705 9120.945 9277.249 -156.304167 9115.429
## [68,] 9458.07 9494.982 9302.314 9379.381 -77.067500 9120.945
## [69,] 9572.03 9507.020 9469.316 9484.398 -15.081667 9302.314
## [70,] 9645.30 9548.832 9604.162 9582.030 22.131667 9469.316
## [71,] 9770.97 9611.592 9729.902 9682.578 47.323750 9604.162
## [72,] 9818.07 9701.592 9883.814 9810.926 72.888750 9729.902
## [73,] 10111.08 9836.355 10105.958 9998.117 107.841250 9883.814
## [74,] 10007.91 9927.007 10190.125 10084.878 105.247083 10105.958
## [75,] 9814.53 9937.897 10083.205 10025.082 58.122917 10190.125
## [76,] 9464.68 9849.550 9785.962 9811.397 -25.435000 10083.205
## [77,] 9462.05 9687.292 9415.385 9524.148 -108.762917 9785.962
## [78,] 9515.51 9564.192 9238.291 9368.652 -130.360417 9415.385
## [79,] 9519.35 9490.397 9227.963 9332.937 -104.973750 9238.291
## [80,] 9530.00 9506.727 9414.353 9451.302 -36.950000 9227.963
## [81,] 9594.39 9539.812 9564.029 9554.342 9.686667 9414.353
## [82,] 9659.17 9575.727 9654.996 9623.289 31.707500 9564.029
## [83,] 9742.16 9631.430 9744.773 9699.436 45.337083 9654.996
## [84,] 9838.26 9708.495 9866.210 9803.124 63.085833 9744.773
## [85,] 10032.66 9818.062 10042.452 9952.696 89.755833 9866.210
## [86,] 10080.80 9923.470 10178.646 10076.576 102.070417 10042.452
## [87,] 10082.43 10008.537 10248.365 10152.434 95.930833 10178.646
## [88,] 10018.19 10053.520 10224.558 10156.143 68.415000 10248.365
## [89,] NA NA NA NA NA 10224.558
## [90,] NA NA NA NA NA 10292.973
## [91,] NA NA NA NA NA 10361.388
## [92,] NA NA NA NA NA 10429.803
## [93,] NA NA NA NA NA 10498.218Hasil peramalan dengan menggunakan N=5 selama 5 periode kedepan pada metode DMA terus naik dari 10224.558 sampai 10498.218
Plot time series
data.gab2 <- ts(data.gab2, start = 2013, frequen = 12)
ts.plot(data.gab2[,1], xlab="Time Period ", ylab="Sales", main= "DMA N=4 Data Premium")
points(data.gab2[,1])
lines(data.gab2[,3],col="green",lwd=2)
lines(data.gab2[,6],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"), lty=8, col=c("black","green","red"), cex=0.8)
Menghitung nilai keakuratan
error.dma = test.ts-data.ramal2[1:length(test.ts)]
SSE.dma = sum(error.dma[8:length(test.ts)]^2)
MSE.dma = mean(error.dma[8:length(test.ts)]^2)
MAPE.dma = mean(abs((error.dma[8:length(test.ts)]/test.ts[8:length(test.ts)])*100))
akurasi.dma <- matrix(c(SSE.dma, MSE.dma, MAPE.dma))
row.names(akurasi.dma)<- c("SSE", "MSE", "MAPE")
colnames(akurasi.dma) <- c("Akurasi m = 4")
akurasi.dma## Akurasi m = 4
## SSE 4.039817e+07
## MSE 2.019909e+06
## MAPE 1.377331e+01akurasi.sma## Akurasi m = 4
## SSE 6.172562e+07
## MSE 2.683723e+06
## MAPE 1.616197e+01Dengan metode SMA, diketahui SSE, MSE, dan MAPE bernilai lebih kecil Metode DMA lebih baik digunakan
Looping for Best M value
SMA
m = 2:30
akurasi.full <- c()
data.premium <- as.numeric(train.ts)
for(i in m){
data.sma <- SMA(train.ts, n = i)
data.ramal <- c(NA, data.sma)
error.sma = train.ts - data.ramal[1:length(train.ts)]
SSE.sma = sum(error.sma[(i+1):length(train.ts)]^2)
MSE.sma = mean(error.sma[(i+1):length(train.ts)]^2)
MAPE.sma = mean(abs((error.sma[(i+1):length(train.ts)]/train.ts[(i+1):length(train.ts)])*100))
tabel <- matrix(c(SSE.sma, MSE.sma, MAPE.sma))
colnames(tabel) <- paste("M =", i)
rownames(tabel) <- c("SSE", "MSE", "MAPE")
akurasi.full <- cbind(akurasi.full, tabel)
}
akurasi.full## M = 2 M = 3 M = 4 M = 5 M = 6
## SSE 5.404590e+06 7.220289e+06 8.430645e+06 9.231573e+06 9.815792e+06
## MSE 6.284407e+04 8.494458e+04 1.003648e+05 1.112238e+05 1.197048e+05
## MAPE 1.958499e+00 2.288197e+00 2.534714e+00 2.699857e+00 2.805653e+00
## M = 7 M = 8 M = 9 M = 10 M = 11
## SSE 1.026683e+07 1.066054e+07 1.099992e+07 1.117805e+07 1.122272e+07
## MSE 1.267510e+05 1.332567e+05 1.392395e+05 1.433083e+05 1.457496e+05
## MAPE 2.898577e+00 2.940943e+00 2.962072e+00 2.939377e+00 2.917784e+00
## M = 12 M = 13 M = 14 M = 15 M = 16
## SSE 1.131593e+07 1.151234e+07 1.196400e+07 1.275476e+07 1.373600e+07
## MSE 1.488938e+05 1.534979e+05 1.616757e+05 1.747228e+05 1.907777e+05
## MAPE 2.917183e+00 2.912221e+00 2.964623e+00 3.068680e+00 3.212885e+00
## M = 17 M = 18 M = 19 M = 20 M = 21
## SSE 1.466445e+07 1.550042e+07 1.632269e+07 1.709597e+07 1.787718e+07
## MSE 2.065416e+05 2.214345e+05 2.365608e+05 2.514113e+05 2.668237e+05
## MAPE 3.331594e+00 3.416996e+00 3.513380e+00 3.596031e+00 3.668926e+00
## M = 22 M = 23 M = 24 M = 25 M = 26
## SSE 1.849690e+07 1.883418e+07 1.837787e+07 1.609895e+07 1.388503e+07
## MSE 2.802561e+05 2.897566e+05 2.871542e+05 2.555389e+05 2.239521e+05
## MAPE 3.731977e+00 3.757461e+00 3.714365e+00 3.577343e+00 3.454139e+00
## M = 27 M = 28 M = 29 M = 30
## SSE 1.258754e+07 1.248655e+07 1.304894e+07 1.341710e+07
## MSE 2.063530e+05 2.081092e+05 2.211685e+05 2.313292e+05
## MAPE 3.379226e+00 3.383222e+00 3.487242e+00 3.565869e+00Nilai error yang terkecil merupakan n optimal, terlihat error yang terkecil ada pada N=2
DMA
m = 2:30
akurasi.full2 <- c()
for(i in m){
data.sma <- SMA(train.ts, n = i)
dma <- SMA(data.sma, n = i)
At <- 2*data.sma - dma
Bt <- 2/(i - 1) * (data.sma - dma)
data.dma <- At + Bt
data.ramal2 <- c(NA, data.dma)
error.dma <- test.ts - data.ramal2[1:length(test.ts)]
SSE.dma <- sum(error.dma[(i*2):length(test.ts)]^2)
MSE.dma <- mean(error.dma[(i*2):length(test.ts)]^2)
MAPE.dma <- mean(abs(error.dma[(i*2):115]/test.ts[(i*2):115]) * 100)
tabel2 <- matrix(c(SSE.dma, MSE.dma, MAPE.dma))
colnames(tabel2) <- paste("M =", i)
rownames <- c("SSE", "MSE", "MAPE")
akurasi.full2 <- cbind(akurasi.full2, tabel2)
}
akurasi.full2 <- as.data.frame(akurasi.full2)
akurasi.full2## M = 2 M = 3 M = 4 M = 5 M = 6 M = 7 M = 8 M = 9
## 1 62652793 50900812 40398173 32035511 25409070 20262264 16816485 13675465
## 2 2610533 2313673 2019909 1779751 1588067 1447305 1401374 1367546
## 3 NA NA NA NA NA NA NA NA
## M = 10 M = 11 M = 12 M = 13 M = 14 M = 15 M = 16 M = 17 M = 18 M = 19
## 1 10495589 6475419 2618966.1 719520.3 NA NA NA NA NA NA
## 2 1311949 1079236 654741.5 359760.2 NA NA NA NA NA NA
## 3 NA NA NA NA NA NA NA NA NA NA
## M = 20 M = 21 M = 22 M = 23 M = 24 M = 25 M = 26 M = 27 M = 28 M = 29 M = 30
## 1 NA NA NA NA NA NA NA NA NA NA NA
## 2 NA NA NA NA NA NA NA NA NA NA NA
## 3 NA NA NA NA NA NA NA NA NA NA NANilai error yang terkecil merupakan n optimal, terlihat error yang terkecil ada pada N=11
Perbandingan SMA dan DMA
comp <- cbind(akurasi.full[,1], akurasi.full2[,10])
colnames(comp) <- c("SMA","DMA")
rownames(comp) <- c("SSE","MSE","MAPE")
comp <- as.data.frame(comp)
comp## SMA DMA
## SSE 5.404590e+06 6475419
## MSE 6.284407e+04 1079236
## MAPE 1.958499e+00 NAJika dibandingkan antara SMA dan DMA pada nilai n optimum, metode DMA menghasilkan error yang lebih kecil dibanding SMA pada ketiga ukuran kebaikan diatas.
Pemulusan Eksponensial
Import data
dataekspo<- read_xlsx("C:/0 SEM5/MPDW/p1/Data Kelompok 15.xlsx")
dataekspo$Premium <- as.numeric(dataekspo$Premium)
str(dataekspo)## tibble [115 × 7] (S3: tbl_df/tbl/data.frame)
## $ Tahun : num [1:115] 2013 2013 2013 2013 2013 ...
## $ Bulan : chr [1:115] "Januari" "Februari" "Maret" "April" ...
## $ Premium : num [1:115] 7798 7773 7576 7421 7545 ...
## $ Medium : chr [1:115] "7697.37" "7645.05" "7503.27" "7290.96" ...
## $ Luar Kualitas: chr [1:115] "7545.32" "7328.44" "7033.14" "6870.91" ...
## $ IHK : num [1:115] 1.03 1.79 2.43 2.32 2.3 3.35 6.75 7.94 7.57 7.66 ...
## $ Harga Gabah : chr [1:115] "4333.19" "4265.58" "3783.15" "3669.04" ...Membentuk objek time series
training<-dataekspo[1:88,3]
testing<-dataekspo[89:115,3]
dataekspo.ts<-ts(dataekspo$Premium, start = 2013, frequency = 12)
training.ts<-ts(training, start = 2013, frequency = 12)
testing.ts<-ts(testing, start = 2020, frequency = 12)Eksplorasi
dataekspo.ts <- as.numeric(dataekspo.ts)
plot(dataekspo.ts, col="red",main="Plot seluruh data")
points(dataekspo.ts)
plot(training.ts, col="blue",main="Plot data training")
points(training.ts)
plot(testing.ts, col="blue",main="Plot data testing")
points(testing.ts)
Single Exponential
Fungsi Holtwinter
ses1<- HoltWinters(training.ts, gamma = FALSE, beta = FALSE, alpha = 0.2)
plot(ses1)
Fungsi Holtwinter optimal
sesopt<- HoltWinters(training.ts, gamma = FALSE, beta = FALSE)
sesopt## Holt-Winters exponential smoothing without trend and without seasonal component.
##
## Call:
## HoltWinters(x = training.ts, beta = FALSE, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 0.9999339
## beta : FALSE
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 10018.19plot(sesopt)
Akurasi Data Training
SSE
ses1<- HoltWinters(training.ts, gamma = FALSE, beta = FALSE, alpha = 0.2)
sse1.train <- ses1$SSE
sse1.train## [1] 9033633sseopt.train <- sesopt$SSE
sseopt.train## [1] 3236764MSE
mse1.train <- sse1.train/length(training.ts)
mse1.train## [1] 102654.9mseopt.train <- sseopt.train/length(training.ts)
mseopt.train## [1] 36781.41RMSE
rmse1.train <- sqrt(mse1.train)
rmse1.train## [1] 320.3981rmseopt.train <- sqrt(mseopt.train)
rmseopt.train## [1] 191.7848akurasi.ses1 <- matrix(c(sse1.train, mse1.train, rmse1.train, sseopt.train, mseopt.train, rmseopt.train), nrow=3, ncol=2)
row.names(akurasi.ses1)<- c("SSE", "MSE", "RMSE")
colnames(akurasi.ses1) <- c("alpha 0.2", "alpha 0.99")
akurasi.ses1## alpha 0.2 alpha 0.99
## SSE 9033633.2419 3236764.3104
## MSE 102654.9232 36781.4126
## RMSE 320.3981 191.7848Metode SES memiliki parameter pemulusan yaitu alpha yang bernilai antara 0 dan 1. Pada awal pemulusan, digunakan alpha 0.2. Lalu dilakukan pemulusan dengan alpha yang optimal secara otomatis dengan program, didapatkan alpha sebesar 0.99 . Parameter pemulusan yang dipilih adalah alpha 0.99 karena memiliki error yang lebih kecil.
Double Exponential
des1<- HoltWinters(training.ts, gamma = FALSE, beta = 0.2, alpha = 0.2)
plot(des1)
desopt<- HoltWinters(training.ts, gamma = FALSE)
desopt## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = training.ts, gamma = FALSE)
##
## Smoothing parameters:
## alpha: 1
## beta : 0.01901135
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 10018.19000
## b 12.97931plot(desopt)
Forecasting
ramalandes1<- forecast(des1, h=10)
ramalandes1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2020 9975.226 9555.017 10395.44 9332.571 10617.88
## Jun 2020 10018.687 9586.545 10450.83 9357.782 10679.59
## Jul 2020 10062.148 9614.275 10510.02 9377.185 10747.11
## Aug 2020 10105.609 9637.985 10573.23 9390.440 10820.78
## Sep 2020 10149.070 9657.586 10640.55 9397.410 10900.73
## Oct 2020 10192.531 9673.100 10711.96 9398.130 10986.93
## Nov 2020 10235.992 9684.635 10787.35 9392.765 11079.22
## Dec 2020 10279.453 9692.361 10866.54 9381.573 11177.33
## Jan 2021 10322.914 9696.477 10949.35 9364.861 11280.97
## Feb 2021 10366.375 9697.198 11035.55 9342.957 11389.79ramalandesopt<- forecast(desopt, h=10)
ramalandesopt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2020 10031.17 9780.149 10282.19 9647.267 10415.07
## Jun 2020 10044.15 9685.762 10402.54 9496.043 10592.25
## Jul 2020 10057.13 9614.031 10500.23 9379.469 10734.79
## Aug 2020 10070.11 9553.640 10586.57 9280.238 10859.98
## Sep 2020 10083.09 9500.251 10665.92 9191.716 10974.46
## Oct 2020 10096.07 9451.661 10740.47 9110.534 11081.60
## Nov 2020 10109.05 9406.575 10811.51 9034.710 11183.38
## Dec 2020 10122.02 9364.158 10879.89 8962.967 11281.08
## Jan 2021 10135.00 9323.833 10946.17 8894.426 11375.58
## Feb 2021 10147.98 9285.188 11010.78 8828.451 11467.52Akurasi data training
SSE
sse.des.train <- des1$SSE
sse.des.train## [1] 9172032sseopt.des.train <- desopt$SSE
sseopt.des.train## [1] 3305982MSE
mse.des.train <- sse1.train/length(training.ts)
mse.des.train## [1] 102654.9mseopt.des.train <- sseopt.train/length(training.ts)
mseopt.des.train## [1] 36781.41RMSE
rmse.des.train <- sqrt(mse.des.train)
rmse.des.train## [1] 320.3981rmseopt.des.train <- sqrt(mseopt.des.train)
rmseopt.des.train## [1] 191.7848akurasi.des1 <- matrix(c(sse.des.train,mse.des.train,rmse.des.train,sseopt.des.train,mseopt.des.train,rmseopt.des.train), nrow=3, ncol=2)
row.names(akurasi.des1)<- c("SSE", "MSE", "RMSE")
colnames(akurasi.des1) <- c("alpha=0.2, beta=0.2", "alpha=1, beta=0.019")
akurasi.des1## alpha=0.2, beta=0.2 alpha=1, beta=0.019
## SSE 9172032.2713 3305981.8294
## MSE 102654.9232 36781.4126
## RMSE 320.3981 191.7848Metode DES memiliki dua parameter pemulusan yaitu alpha dan beta yang bernilai antara 0 dan 1. Pada awal pemulusan, digunakan alpha dan beta sebesar 0.2. Lalu dilakukan pemulusan dengan alpha dan beta yang optimal secara otomatis dengan program, didapatkan alpha sebesar 1 dan beta sebesar 0.019. Parameter pemulusan yang dipilih adalah alpha = 1 dan beta = 0.0181 karena memiliki error yang lebih kecil.
Perbandingan Akurasi Test Single Exponential dan Double Exponential
ramalanopt<- forecast(sesopt, h=10)
selisihses<-ramalanopt$mean-testing.ts
selisihses## Jan Feb Mar Apr May Jun Jul Aug Sep
## 2020 147.0742 205.2942 303.6942 230.3342 238.0042
## 2021 391.1142 481.0342
## Oct Nov Dec
## 2020 246.0942 411.2242 468.4342
## 2021SSEtestingses<-sum(selisihses^2)
str(testing.ts)## Time-Series [1:27, 1] from 2020 to 2022: 9827 9919 9932 9963 9871 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr "Premium"selisihdes<-ramalandesopt$mean-testing.ts
selisihdes## Jan Feb Mar Apr May Jun Jul Aug Sep
## 2020 160.0493 231.2486 342.6279 282.2472 302.8965
## 2021 507.9238 610.8231
## Oct Nov Dec
## 2020 323.9659 502.0752 572.2645
## 2021SSEtestingdes<-sum(selisihdes^2)
akurasi <- matrix(c(SSEtestingses, SSEtestingdes), nrow=1, ncol=2)
row.names(akurasi)<- "SSE"
colnames(akurasi) <- c("SES", "DES")
akurasi## SES DES
## SSE 1099169 1683507Dibandingkan degan metode moving average, metode exponential memiliki error yang jauh lebih kecil, hal ini menandakan bahwa adanya perbedaan pengaruh harga beras antar periode sehingga metode exponential lebih cocok digunakan.
Pemulusan Winter
Import data
winter <- data
str(winter)## tibble [115 × 7] (S3: tbl_df/tbl/data.frame)
## $ Tahun : num [1:115] 2013 2013 2013 2013 2013 ...
## $ Bulan : chr [1:115] "Januari" "Februari" "Maret" "April" ...
## $ Premium : num [1:115] 7798 7773 7576 7421 7545 ...
## $ Medium : chr [1:115] "7697.37" "7645.05" "7503.27" "7290.96" ...
## $ Luar Kualitas: chr [1:115] "7545.32" "7328.44" "7033.14" "6870.91" ...
## $ IHK : num [1:115] 1.03 1.79 2.43 2.32 2.3 3.35 6.75 7.94 7.57 7.66 ...
## $ Harga Gabah : chr [1:115] "4333.19" "4265.58" "3783.15" "3669.04" ...Membagi data menjadi training dan testing
training<-winter[1:88,3]
testing<-winter[89:115,3]Membentuk objek time series
winter.ts<-ts(winter$Premium,start = 2013, frequency = 12,)
training.ts<-ts(training,start = 2013, frequency = 12)
testing.ts<-ts(testing, start = 2020, frequency = 12)Membuat plot time series
plot(winter.ts, col="red")
points(winter.ts)
plot(training.ts, col="blue")
points(training.ts)
Winter Aditif
Pemulusan
aditif <- HoltWinters(training.ts, seasonal = "additive")
aditif ## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = training.ts, seasonal = "additive")
##
## Smoothing parameters:
## alpha: 0.9514662
## beta : 0
## gamma: 1
##
## Coefficients:
## [,1]
## a 10245.978564
## b 45.393039
## s1 -212.784340
## s2 -111.756023
## s3 24.000143
## s4 -7.181151
## s5 -67.295108
## s6 -43.632850
## s7 4.759410
## s8 44.169309
## s9 234.497408
## s10 240.376417
## s11 76.866665
## s12 -227.788564Forecasting
ramalan1 <- forecast(aditif, h=27)
ramalan1## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2020 10078.59 9872.835 10284.34 9763.916 10393.26
## Jun 2020 10225.01 9941.004 10509.01 9790.661 10659.36
## Jul 2020 10406.16 10061.218 10751.10 9878.618 10933.70
## Aug 2020 10420.37 10023.749 10816.99 9813.790 11026.95
## Sep 2020 10405.65 9963.345 10847.95 9729.203 11082.09
## Oct 2020 10474.70 9991.013 10958.40 9734.962 11214.45
## Nov 2020 10568.49 10046.683 11090.30 9770.455 11366.52
## Dec 2020 10653.29 10095.971 11210.61 9800.943 11505.64
## Jan 2021 10889.01 10298.309 11479.72 9985.610 11792.42
## Feb 2021 10940.29 10317.987 11562.58 9988.562 11892.01
## Mar 2021 10822.17 10169.803 11474.53 9824.462 11819.88
## Apr 2021 10562.91 9881.801 11244.01 9521.245 11604.57
## May 2021 10623.30 9911.799 11334.81 9535.151 11711.46
## Jun 2021 10769.73 10031.779 11507.67 9641.135 11898.32
## Jul 2021 10950.87 10187.403 11714.35 9783.246 12118.50
## Aug 2021 10965.09 10176.915 11753.26 9759.683 12170.49
## Sep 2021 10950.37 10138.246 11762.48 9708.336 12192.39
## Oct 2021 11019.42 10184.039 11854.80 9741.815 12297.03
## Nov 2021 11113.21 10255.192 11971.22 9800.988 12425.42
## Dec 2021 11198.01 10317.945 12078.07 9852.068 12543.95
## Jan 2022 11433.73 10532.156 12335.30 10054.891 12812.57
## Feb 2022 11485.00 10562.418 12407.59 10074.032 12895.97
## Mar 2022 11366.89 10423.760 12310.01 9924.499 12809.27
## Apr 2022 11107.62 10144.394 12070.85 9634.491 12580.75
## May 2022 11168.02 10183.061 12152.98 9661.655 12674.39
## Jun 2022 11314.44 10310.216 12318.67 9778.611 12850.27
## Jul 2022 11495.59 10472.462 12518.72 9930.850 13060.33Akurasi data training
SSE
sse1.train <- aditif$SSE
sse1.train## [1] 1944917Winter Multiplikatif
Pemulusan
multi <- HoltWinters(training.ts,seasonal = "multiplicative")
multi## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = training.ts, seasonal = "multiplicative")
##
## Smoothing parameters:
## alpha: 0.9318329
## beta : 0
## gamma: 1
##
## Coefficients:
## [,1]
## a 1.028539e+04
## b 4.539304e+01
## s1 9.732413e-01
## s2 9.851614e-01
## s3 1.000882e+00
## s4 9.992562e-01
## s5 9.922297e-01
## s6 9.941485e-01
## s7 1.000243e+00
## s8 1.005849e+00
## s9 1.028732e+00
## s10 1.029457e+00
## s11 1.010173e+00
## s12 9.740216e-01Forecasting
ramalan2 <- forecast(multi, h=27)
ramalan2## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## May 2020 10054.34 9845.032 10263.65 9734.230 10374.45
## Jun 2020 10222.21 9934.474 10509.94 9782.157 10662.25
## Jul 2020 10430.76 10079.609 10781.92 9893.719 10967.81
## Aug 2020 10459.17 10057.965 10860.38 9845.578 11072.77
## Sep 2020 10430.67 9987.007 10874.33 9752.147 11109.19
## Oct 2020 10495.97 10010.565 10981.37 9753.608 11238.33
## Nov 2020 10605.72 10079.897 11131.53 9801.546 11409.88
## Dec 2020 10710.82 10147.285 11274.35 9848.969 11572.67
## Jan 2021 11001.19 10392.946 11609.43 10070.961 11931.41
## Feb 2021 11055.66 10416.513 11694.81 10078.167 12033.16
## Mar 2021 10894.43 10237.457 11551.39 9889.679 11899.17
## Apr 2021 10548.76 9919.530 11177.98 9586.438 11511.07
## May 2021 10584.48 9917.499 11251.47 9564.419 11604.55
## Jun 2021 10758.84 10056.280 11461.40 9684.368 11833.31
## Jul 2021 10975.96 10236.272 11715.65 9844.705 12107.21
## Aug 2021 11003.49 10239.653 11767.32 9835.304 12171.67
## Sep 2021 10971.15 10187.911 11754.39 9773.287 12169.02
## Oct 2021 11037.50 10228.894 11846.10 9800.845 12274.15
## Nov 2021 11150.56 10314.036 11987.09 9871.204 12429.92
## Dec 2021 11258.72 10395.263 12122.18 9938.177 12579.26
## Jan 2022 11561.56 10657.466 12465.64 10178.870 12944.24
## Feb 2022 11616.43 10690.923 12541.93 10200.992 13031.86
## Mar 2022 11444.68 10515.580 12373.79 10023.741 12865.63
## Apr 2022 11079.32 10189.462 11969.18 9718.399 12440.24
## May 2022 11114.62 10194.882 12034.37 9708.000 12521.25
## Jun 2022 11295.47 10344.479 12246.47 9841.054 12749.89
## Jul 2022 11521.16 10535.780 12506.53 10014.153 13028.16Akurasi data training
SSE
sse2.train <- multi$SSE
sse2.train## [1] 2011779Akurasi data testing
selisih1<-as.numeric(ramalan1$mean)-as.numeric(testing.ts)
selisih1## [1] 251.9573 306.0586 474.4478 457.7196 534.5287 661.8039 853.9892
## [8] 865.4322 1108.8233 1168.1854 1215.1987 1013.1465 996.2237 1232.5651
## [15] 1549.2643 1465.7260 1494.8051 1569.9504 1573.9157 1525.4687 1609.4998
## [22] 1658.1218 1580.2551 1530.8729 1655.3902 1817.0416 1867.0208SSEtesting1<-sum(selisih1^2)
selisih2<-as.numeric(ramalan2$mean)-as.numeric(testing.ts)
selisih2## [1] 227.7125 303.2557 499.0526 496.5242 559.5484 683.0679 891.2152
## [8] 922.9575 1220.9981 1283.5635 1287.4558 998.9956 957.4031 1221.6794
## [15] 1574.3497 1504.1255 1515.5922 1588.0269 1611.2741 1586.1802 1737.3256
## [22] 1789.5455 1658.0538 1502.5712 1601.9936 1798.0730 1892.5868SSEtesting2<-sum(selisih2^2)
akurasi <- matrix(c(SSEtesting1, SSEtesting2), nrow=1, ncol=2)
row.names(akurasi)<- "SSE"
colnames(akurasi) <- c("Aditif", "Multiplikatif")
akurasi## Aditif Multiplikatif
## SSE 44163120 46549188Kesimpulan
- Metode Exponential lebih cocok digunakan dibandingkan dengan metode moving average. Adanya perbedaan pengaruh harga beras antar periode membuat metode exponential lebih cocok digunakan.
- Dalam Metode Eksponensial, metode Double Exponential Smoothing lebih cocok digunakan dibandingkan Single Exponential Smoothing.
- Pada metode Winter, multiplikative lebih cocok digunakan dibandingkan additive.
Kestasioneran Data
Library
library("tseries")
library("forecast")
library("TTR")
library("TSA")
library("graphics")
library("astsa")
library("car")Import Data
data.premium <- ts(data$Premium, start = 2013, frequency = 12)
ts.plot(data.premium, main = "Harga Beras Premium Tahun 2013 - 2022", ylab = "Harga", lwd = 1.5)
points(data.premium)
legend("topleft", c("Premium"), cex = 0.7,
col = c("black", "red", "blue"), lty = 1)
Evaluasi periode musiman
#Plot tahunan
library(forecast)
seasonplot(data.premium,12,main="Xt", ylab="Xt",year.labels = TRUE, col=rainbow(18))
monthplot(data.premium,ylab="Xt")
boxplot(data$Premium~data$Bulan,ylab="Xt",xlab="Bulan")
Berdasarkan Boxplot terlihat bahwa ragam data dari tahun ke tahun cenderung tidak homogen, sehingga untuk memastikan, selanjutnya dilakukan pengujian dengan fligner test.
Plot Transformasi Log
logprem <- log10(data[,"Premium"])
ts.plot(logprem, type="l", ylab="Log Premium")
Transformasi logaritma belum berhasil mengatasi ketidakstasioneran dalam ragam dan rataan.
Uji Homogenitas
library(car)
fligner.test(`Premium` ~ Tahun, data=data)##
## Fligner-Killeen test of homogeneity of variances
##
## data: Premium by Tahun
## Fligner-Killeen:med chi-squared = 17.725, df = 9, p-value = 0.0385fligner.test(log10(Premium) ~ Bulan, data=data)##
## Fligner-Killeen test of homogeneity of variances
##
## data: log10(Premium) by Bulan
## Fligner-Killeen:med chi-squared = 1.8348, df = 11, p-value = 0.999Berdasarkan hasil uji homogenitas pada data yang belum dilakukan transformasi logaritma di atas, diperoleh nilai p-value=0.0385 < 0.05 sehingga dapat disimpulkan bahwa data harga beras premium tidak homogen.
Oleh karena itu, berdasarkan hasil plot data dan pengujian homogenitas dapat disimpulkan bahwa data harga beras premium tidak stasioner dalam ragam dan rataan.
Sedangkan berdasarkan hasil uji homogenitas pada data yang telah dilakukan transformasi logaritma di atas, diperoleh nilai p-value=0.999 > 0.05 sehingga dapat disimpulkan bahwa data harga beras premium sudah homogen.
Pembagian Data Menjadi Training-Testing
log.train <- subset(log10(data[,"Premium"]),start=1,end=88)
log.test <- subset(log10(data[,"Premium"]),start=89,end=115)Untuk analisis berikutnya, data yang digunakan adalah yang telah ditransformasi logaritma (sebut saja Yt)
ACF Yt
acf(log.train,main="ACF",lag.max=50,xaxt="n")
axis(1, at=0:50, labels=0:50)
Plot ACF data setelah transformasi logaritma menunjukkan pola nonstasioner.
ACF Musiman Yt
acf0<-acf(log.train,main="ACF",lag.max=50,xaxt="n")
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 data differencing non musiman
log.train <- ts(log.train)
d1 <- diff(log.train)
ts.plot(d1, type="l", ylab="d1 Beras Premium")
differencing non musiman d = 1 berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen non musimannya.
Plot ACF Data d1
acf1 <- acf(d1,lag.max=50,xaxt="n", main="ACF d1")
axis(1, at=0:50, labels=0:50)
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 differencing non musiman d = 1 mengkonfirmasi kestasioneran komponen non musiman
Plot Data Differencing Musiman d1
D1 <- diff(log.train,12)
ts.plot(D1, type="l", ylab="D1 Beras Premium")
Plot ACF Data D1
D1 <- diff(log.train,12)
acf2 <- acf(D1,lag.max=50,xaxt="n", main="ACF d1")
ts.plot(D1, type="l", ylab="D1 Harga Beras Premium")
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")
differencing musiman D = 1 berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen musimannya
Plot Data d1D1
d1D1 <- diff(D1)
ts.plot(d1D1, type="l", ylab="d1 D1 Harga Beras Premium")
# Correlogram Data d1D1
acf3 <- acf(d1D1,lag.max=50,xaxt="n", main="ACF d1D1")
axis(1, at=0:50, labels=0:50)
pacf3 <- pacf(d1D1,lag.max=50,xaxt="n", main="PACF d1D1")
axis(1, at=0:50, labels=0:50)
Kedua komponen telah stasioner. ARIMA(0,1,4)
Correlogram Musiman Data d1D1
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 = "ACF", xlab = "Lag")
Identifikasi komponen musiman adalah ARIMA()12, sehingga model tentatif adalah ARIMA(0,1,4)?(0,1,1)12.
Pengepasan ARIMA (0,1,4)(0,1,1)12
tmodel <- Arima(log.train,order=c(0,1,4),seasonal=c(0,1,1))
summary(tmodel)## Series: log.train
## ARIMA(0,1,4)
##
## Coefficients:
## ma1 ma2 ma3 ma4
## 0.4381 0.0439 -0.2887 -0.3129
## s.e. 0.0881 0.0920 0.0887 0.0823
##
## sigma^2 = 5.066e-05: log likelihood = 403.85
## AIC=-797.7 AICc=-797.15 BIC=-784.02
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0009743268 0.006960857 0.005191812 0.0245862 0.1310721 0.90942
## ACF1
## Training set -0.01831919Diagnosa Sisaan
checkresiduals(tmodel)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,4)
## Q* = 6.0052, df = 6, p-value = 0.4226
##
## Model df: 4. Total lags used: 10Diagnosa Sisaan
tsdisplay(residuals(tmodel), lag.max=45, main='ARIMA(0,1,2)(0,1,1)12 Model Residuals')
Diagnosa Sisaan
library(tseries)
jarque.bera.test(residuals(tmodel))##
## Jarque Bera Test
##
## data: residuals(tmodel)
## X-squared = 9.5284, df = 2, p-value = 0.00853Validasi Model
ramalan_arima = forecast(tmodel, 12)plot(ramalan_arima)
Validasi Model
Plot ACF dan Proses Differencing
ACF
acf(train.ts[,1], lag.max = 50)
Berdasarkan plot ACF di atas terlihat bahwa nilai autokorelasi beras premium turun secara perlahan yang mengindikasikan bahwa pola data tidak stasioner.
ACF Musiman
acf0.Zt<-acf(train.ts, main="ACF", lag.max=50, xaxt="n")
acf0.Zt$lag <- acf0.Zt$lag * 12
acf0.Zt.1 <- as.data.frame(cbind(acf0.Zt$acf,acf0.Zt$lag))
acf0.Zt.2 <- acf0.Zt.1[which(acf0.Zt.1$V2%%12==0),]
barplot(height = acf0.Zt.2$V1, names.arg = acf0.Zt.2$V2, ylab = "ACF", xlab = "Lag")
Berdasarkan plot ACF di atas terlihat bahwa…
Differencing
# Plot data differencing non musiman
d1.Zt <- diff(train.ts)
ts.plot(d1.Zt, type="l", ylab="d1 Zt")
Plot data hasil differencing untuk komponen non musiman di atas menunjukkan bahwa, dengan differencing ordo 1 berhasil mengatasi ketidakstasioneran dalam rataan untuk komponen non musiman.
Plot ACF Data Differencing
acf(d1.Zt[,1], lag.max = 50)
acf1.Zt <- acf(d1.Zt,lag.max=50,xaxt="n", main="ACF d1.Zt")
acf1.Zt$lag <- acf1.Zt$lag * 12
acf1.Zt.1 <- as.data.frame(cbind(acf1.Zt$acf,acf1.Zt$lag))
acf1.Zt.2 <- acf1.Zt.1[which(acf1.Zt.1$V2%%12==0),]
barplot(height = acf1.Zt.2$V1, names.arg=acf1.Zt.2$V2, ylab="ACF", xlab="Lag")
ACF untuk data differencing beras premium untuk komponen non musiman menunjukkan bahwa komponen non musiman telah stasioner. Sehingga, selanjutnya dapat dilakukan identifikasi model ARIMA Musiman.
Identifikasi Model (Komponen Non Musiman)
# Correlogram Data d1.Zt
acf3.Zt<- acf(d1.Zt,lag.max=50,xaxt="n", main="ACF d1.Zt")
axis(1, at=0:50, labels=0:50)
#PACF
pacf3.Zt <- pacf(d1.Zt,lag.max=50,xaxt="n", main="PACF d1.Zt")
axis(1, at=0:50, labels=0:50)
#EACF
eacf(d1.Zt)## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o x x o o o o o o x x x o
## 1 x o x x o o o o o o o x x o
## 2 x o o x o o o o o o o o o o
## 3 x o o x o o o o o o o o o o
## 4 o o o x o o o o o o o o o o
## 5 o x x o o o o o o o o o o o
## 6 x o x o o o o o o o o o o o
## 7 x o x o o o o x o o o o o oBerdasarkan plot ACF (cut off setelah lag ke-1), PACF (cut off setelah lag ke-2), dan EACF, diperoleh kandidat model untuk komponen non musiman:
ARIMA (0,1,4)
ARIMA (5,1,3)
ARIMA (7,1,8)
ARIMA (0,1,1)
ARIMA (2,1,0)
Identifikasi Model (Komponen Musiman)
# Correlogram Data d1.Zt
acf3.Zt<- acf(d1.Zt,lag.max=120,xaxt="n", main="ACF d1.Zt")
axis(1, at=0:120, labels=0:120)
#PACF
pacf3.Zt <- pacf(d1.Zt,lag.max=120,xaxt="n", main="PACF d1.Zt")
axis(1, at=0:120, labels=0:120)
# Correlogram Musiman Data d1.Zt
acf3.Zt$lag <- acf3.Zt$lag * 12
acf3.Zt.1 <- as.data.frame(cbind(acf3.Zt$acf,acf3.Zt$lag))
acf3.Zt.2 <- acf3.Zt.1[which(acf3.Zt.1$V2%%12==0),]
barplot(height = acf3.Zt.2$V1, names.arg=acf3.Zt.2$V2, ylab="ACF", xlab="Lag")
# PACF
pacf3.Zt$lag <- pacf3.Zt$lag * 12
pacf3.Zt.1 <- as.data.frame(cbind(pacf3.Zt$acf,pacf3.Zt$lag))
pacf3.Zt.2 <- pacf3.Zt.1[which(pacf3.Zt.1$V2%%12==0),]
barplot(height = pacf3.Zt.2$V1, names.arg=pacf3.Zt.2$V2, ylab="PACF", xlab="Lag")
ARIMA (0,0,2)12
ARIMA (1,0,1)12
ARIMA (1,0,2)12
Pendugaan Parameter Model ARIMA Musiman
Model tentatif yang diperoleh berdasarkan hasil identifikasi menggunakan plot ACF, PACF, dan EACF adalah: ARIMA (0,1,4)x(0,0,2)12
ARIMA (5,1,3)x(0,0,2)12
ARIMA (7,1,8)x(0,0,2)12
ARIMA (0,1,1)x(0,0,2)12
ARIMA (2,1,0)x(0,0,2)12
ARIMA (0,1,4)x(1,0,1)12
ARIMA (5,1,3)x(1,0,1)12
ARIMA (7,1,8)x(1,0,1)12
ARIMA (0,1,1)x(1,0,1)12
ARIMA (2,1,0)x(1,0,1)12
ARIMA (0,1,4)x(1,0,2)12
ARIMA (5,1,3)x(1,0,2)12
ARIMA (7,1,8)x(1,0,2)12
ARIMA (0,1,1)x(1,0,2)12
ARIMA (2,1,0)x(1,0,2)12
Model untuk komponen musiman ARIMA (0,0,2)12
#ARIMA(0,1,4)x(0,0,2)12
modela <- Arima(train.ts,order=c(0,1,4),seasonal=list(order=c(0,0,2),period=12))
summary(modela)## Series: train.ts
## ARIMA(0,1,4)(0,0,2)[12]
##
## Coefficients:
## ma1 ma2 ma3 ma4 sma1 sma2
## 0.4232 0.0275 -0.2798 -0.3499 0.2778 0.0732
## s.e. 0.1038 0.1076 0.1057 0.0951 0.1165 0.1214
##
## sigma^2 = 24506: log likelihood = -560.84
## AIC=1135.68 AICc=1137.09 BIC=1152.94
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 23.38814 150.1901 111.1859 0.2528414 1.224487 0.8326371
## ACF1
## Training set -0.03094113#ARIMA(5,1,3)x(0,0,2)12
modelb <- Arima(train.ts,order=c(5,1,3),seasonal=list(order=c(0,0,2),period=12))
summary(modelb)## Series: train.ts
## ARIMA(5,1,3)(0,0,2)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3 sma1
## 0.1842 -0.2464 0.4393 -0.4857 0.1750 0.2591 0.2453 -0.6750 0.2934
## s.e. 0.2899 0.1567 0.1820 0.1037 0.1693 0.2648 0.2472 0.2324 0.1215
## sma2
## 0.1234
## s.e. 0.1308
##
## sigma^2 = 23283: log likelihood = -556.96
## AIC=1135.93 AICc=1139.45 BIC=1163.05
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 20.66688 142.7323 111.3862 0.2239564 1.22949 0.8341367
## ACF1
## Training set 0.0002868045#ARIMA(7,1,8)x(0,0,2)12
modelc <- Arima(train.ts,order=c(7,1,8),seasonal=list(order=c(0,0,2),period=12))
summary(modelc)## Series: train.ts
## ARIMA(7,1,8)(0,0,2)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ma1
## 0.7363 -0.0724 -0.3276 -0.1691 0.2922 0.2848 -0.5584 -0.3756
## s.e. 0.2579 0.3184 0.3217 0.2773 0.3295 0.3135 0.1914 0.2787
## ma2 ma3 ma4 ma5 ma6 ma7 ma8 sma1 sma2
## -0.2374 0.0397 0.1948 0.0060 -0.5676 0.3311 0.6113 0.1667 0.0071
## s.e. 0.2656 0.2632 0.3166 0.2448 0.2440 0.2031 0.1444 0.1444 0.1569
##
## sigma^2 = 20391: log likelihood = -551.09
## AIC=1138.17 AICc=1148.23 BIC=1182.56
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 18.65106 127.3572 101.1576 0.2021296 1.116991 0.7575376
## ACF1
## Training set -0.02963039#ARIMA(0,1,1)x(0,0,2)12
modeld <- Arima(train.ts,order=c(0,1,1),seasonal=list(order=c(0,0,2),period=12))
summary(modeld)## Series: train.ts
## ARIMA(0,1,1)(0,0,2)[12]
##
## Coefficients:
## ma1 sma1 sma2
## 0.4459 0.2772 0.0889
## s.e. 0.0970 0.1103 0.1147
##
## sigma^2 = 26905: log likelihood = -566.23
## AIC=1140.46 AICc=1140.95 BIC=1150.33
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 13.99828 160.2572 120.8831 0.151143 1.329644 0.905256 0.003017726#ARIMA(2,1,0)x(0,0,2)12
modele <- Arima(train.ts,order=c(2,1,0),seasonal=list(order=c(0,0,2),period=12))
summary(modele)## Series: train.ts
## ARIMA(2,1,0)(0,0,2)[12]
##
## Coefficients:
## ar1 ar2 sma1 sma2
## 0.4694 -0.2523 0.2467 0.0855
## s.e. 0.1064 0.1048 0.1121 0.1149
##
## sigma^2 = 26550: log likelihood = -565.07
## AIC=1140.14 AICc=1140.89 BIC=1152.47
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 15.8629 158.2453 118.8865 0.1699197 1.309519 0.8903041 -0.05123798printstatarima.Zt <- 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.Zt(modela)##
## Coefficients:
## s.e. t sign.
## ma1 0.4232 0.1038 4.0771 0.0001
## ma2 0.0275 0.1076 0.2556 0.7989
## ma3 -0.2798 0.1057 -2.6471 0.0096
## ma4 -0.3499 0.0951 -3.6793 0.0004
## sma1 0.2778 0.1165 2.3845 0.0193
## sma2 0.0732 0.1214 0.6030 0.5481printstatarima.Zt(modelb)##
## Coefficients:
## s.e. t sign.
## ar1 0.1842 0.2899 0.6354 0.5268
## ar2 -0.2464 0.1567 -1.5724 0.1195
## ar3 0.4393 0.1820 2.4137 0.0179
## ar4 -0.4857 0.1037 -4.6837 0.0000
## ar5 0.1750 0.1693 1.0337 0.3042
## ma1 0.2591 0.2648 0.9785 0.3306
## ma2 0.2453 0.2472 0.9923 0.3238
## ma3 -0.6750 0.2324 -2.9045 0.0047
## sma1 0.2934 0.1215 2.4148 0.0178
## sma2 0.1234 0.1308 0.9434 0.3481printstatarima.Zt(modelc)##
## Coefficients:
## s.e. t sign.
## ar1 0.7363 0.2579 2.8550 0.0054
## ar2 -0.0724 0.3184 -0.2274 0.8207
## ar3 -0.3276 0.3217 -1.0183 0.3113
## ar4 -0.1691 0.2773 -0.6098 0.5436
## ar5 0.2922 0.3295 0.8868 0.3776
## ar6 0.2848 0.3135 0.9085 0.3661
## ar7 -0.5584 0.1914 -2.9175 0.0045
## ma1 -0.3756 0.2787 -1.3477 0.1813
## ma2 -0.2374 0.2656 -0.8938 0.3739
## ma3 0.0397 0.2632 0.1508 0.8805
## ma4 0.1948 0.3166 0.6153 0.5400
## ma5 0.0060 0.2448 0.0245 0.9805
## ma6 -0.5676 0.2440 -2.3262 0.0223
## ma7 0.3311 0.2031 1.6302 0.1067
## ma8 0.6113 0.1444 4.2334 0.0001
## sma1 0.1667 0.1444 1.1544 0.2515
## sma2 0.0071 0.1569 0.0453 0.9640printstatarima.Zt(modeld)##
## Coefficients:
## s.e. t sign.
## ma1 0.4459 0.0970 4.5969 0.0000
## sma1 0.2772 0.1103 2.5131 0.0138
## sma2 0.0889 0.1147 0.7751 0.4404printstatarima.Zt(modele)##
## Coefficients:
## s.e. t sign.
## ar1 0.4694 0.1064 4.4117 0.0000
## ar2 -0.2523 0.1048 -2.4074 0.0182
## sma1 0.2467 0.1121 2.2007 0.0304
## sma2 0.0855 0.1149 0.7441 0.4588Berdasarkan hasil di atas, dengan mempertimbangkan uji signifikansi penduga parameter dan ukuran keakuratan ramalan dapat disimpulkan bahwa dengan menggunakan model musiman (0,0,2)12 model terbaik adalah ARIMA model C
Akurasi Model
akurasiA<-data.frame("Model ARIMA"=c("ARIMA (0,1,1)x(0,0,2)12", "ARIMA(3,1,0)x(0,0,2)12", "ARIMA(3,1,1)x(0,0,2)12", "ARIMA (1,1,2)x(0,0,2)12"),
"MAD"=c(27.95632,27.93469, 27.76326,27.84618 ),
"MSE"=c(36.37349^2,36.34116^2, 36.09141^2, 36.28828^2),
"MAPE"=c(6.441591, 6.450431, 6.358971, 6.38023),
"AIC"=c(2405.97,2409.52,2408.34,2408.92))
akurasiA## Model.ARIMA MAD MSE MAPE AIC
## 1 ARIMA (0,1,1)x(0,0,2)12 27.95632 1323.031 6.441591 2405.97
## 2 ARIMA(3,1,0)x(0,0,2)12 27.93469 1320.680 6.450431 2409.52
## 3 ARIMA(3,1,1)x(0,0,2)12 27.76326 1302.590 6.358971 2408.34
## 4 ARIMA (1,1,2)x(0,0,2)12 27.84618 1316.839 6.380230 2408.92Model untuk komponen musiman ARIMA (1,0,1)12
modela2 <- Arima(train.ts,order=c(0,1,4),seasonal=list(order=c(1,0,1),period=12))
summary(modela2)## Series: train.ts
## ARIMA(0,1,4)(1,0,1)[12]
##
## Coefficients:
## ma1 ma2 ma3 ma4 sar1 sma1
## 0.3433 0.0174 -0.1953 -0.3043 0.9996 -0.9779
## s.e. 0.1035 0.1071 0.1028 0.0973 0.0023 0.0660
##
## sigma^2 = 17338: log likelihood = -554.73
## AIC=1123.46 AICc=1124.88 BIC=1140.72
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3.661292 126.3282 96.74171 0.04381646 1.065633 0.7244686
## ACF1
## Training set 0.005728288#ARIMA(5,1,3)x(1,0,1)12
modelb2 <- Arima(train.ts,order=c(5,1,3),seasonal=list(order=c(1,0,1),period=12))
summary(modelb2)## Series: train.ts
## ARIMA(5,1,3)(1,0,1)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## -0.0590 -0.3089 0.4069 -0.4669 0.0581 0.4473 0.3928 -0.5644
## s.e. 0.5162 0.1932 0.2426 0.0987 0.2726 0.6217 0.4465 0.3030
## sar1 sma1
## 0.9987 -0.9536
## s.e. 0.0085 0.1399
##
## sigma^2 = 15847: log likelihood = -549.66
## AIC=1121.32 AICc=1124.84 BIC=1148.44
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.155707 117.7543 94.3834 0.01732897 1.042683 0.706808 0.01425039#ARIMA(7,1,8)x(1,0,1)12
modelc2 <- Arima(train.ts,order=c(7,1,8),seasonal=list(order=c(1,0,1),period=12))
summary(modelc2)## Series: train.ts
## ARIMA(7,1,8)(1,0,1)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ma1
## 0.6304 -0.1539 -0.0171 -0.2284 0.3283 0.2347 -0.4719 -0.3145
## s.e. 0.3098 0.4593 0.4599 0.3770 0.3928 0.4077 0.2431 0.3246
## ma2 ma3 ma4 ma5 ma6 ma7 ma8 sar1
## -0.0820 -0.1953 0.1221 -0.1179 -0.4536 0.4224 0.5803 0.9996
## s.e. 0.3789 0.3742 0.2980 0.3627 0.3852 0.1494 0.2243 0.0022
## sma1
## -0.9801
## s.e. 0.0593
##
## sigma^2 = 15314: log likelihood = -546.45
## AIC=1128.89 AICc=1138.95 BIC=1173.28
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3.771469 110.3702 87.65015 0.04371762 0.9696608 0.6563848
## ACF1
## Training set 0.003824925#ARIMA(0,1,1)x(1,0,1)12
modeld2 <- Arima(train.ts,order=c(0,1,1),seasonal=list(order=c(1,0,1),period=12))
summary(modeld2)## Series: train.ts
## ARIMA(0,1,1)(1,0,1)[12]
##
## Coefficients:
## ma1 sar1 sma1
## 0.3757 0.9994 -0.9720
## s.e. 0.1048 0.0035 0.0817
##
## sigma^2 = 18440: log likelihood = -558.85
## AIC=1125.7 AICc=1126.19 BIC=1135.57
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3.105948 132.6735 102.3998 0.03541785 1.127001 0.7668407
## ACF1
## Training set -0.0007157275#ARIMA(2,1,0)x(1,0,1)12
modele2 <- Arima(train.ts,order=c(2,1,0),seasonal=list(order=c(1,0,1),period=12))
summary(modele2)## Series: train.ts
## ARIMA(2,1,0)(1,0,1)[12]
##
## Coefficients:
## ar1 ar2 sar1 sma1
## 0.3750 -0.1831 0.9996 -0.9775
## s.e. 0.1088 0.1072 0.0024 0.0653
##
## sigma^2 = 18502: log likelihood = -558.34
## AIC=1126.68 AICc=1127.42 BIC=1139.01
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3.210448 132.1029 100.8782 0.03617671 1.111796 0.7554458
## ACF1
## Training set -0.0209211printstatarima.Zt(modela2)##
## Coefficients:
## s.e. t sign.
## ma1 0.3433 0.1035 3.3169 0.0013
## ma2 0.0174 0.1071 0.1625 0.8713
## ma3 -0.1953 0.1028 -1.8998 0.0608
## ma4 -0.3043 0.0973 -3.1274 0.0024
## sar1 0.9996 0.0023 434.6087 0.0000
## sma1 -0.9779 0.0660 -14.8167 0.0000printstatarima.Zt(modelb2)##
## Coefficients:
## s.e. t sign.
## ar1 -0.0590 0.5162 -0.1143 0.9093
## ar2 -0.3089 0.1932 -1.5989 0.1135
## ar3 0.4069 0.2426 1.6772 0.0971
## ar4 -0.4669 0.0987 -4.7305 0.0000
## ar5 0.0581 0.2726 0.2131 0.8317
## ma1 0.4473 0.6217 0.7195 0.4738
## ma2 0.3928 0.4465 0.8797 0.3814
## ma3 -0.5644 0.3030 -1.8627 0.0659
## sar1 0.9987 0.0085 117.4941 0.0000
## sma1 -0.9536 0.1399 -6.8163 0.0000printstatarima.Zt(modelc2)##
## Coefficients:
## s.e. t sign.
## ar1 0.6304 0.3098 2.0349 0.0449
## ar2 -0.1539 0.4593 -0.3351 0.7384
## ar3 -0.0171 0.4599 -0.0372 0.9704
## ar4 -0.2284 0.3770 -0.6058 0.5462
## ar5 0.3283 0.3928 0.8358 0.4056
## ar6 0.2347 0.4077 0.5757 0.5663
## ar7 -0.4719 0.2431 -1.9412 0.0555
## ma1 -0.3145 0.3246 -0.9689 0.3353
## ma2 -0.0820 0.3789 -0.2164 0.8292
## ma3 -0.1953 0.3742 -0.5219 0.6031
## ma4 0.1221 0.2980 0.4097 0.6830
## ma5 -0.1179 0.3627 -0.3251 0.7459
## ma6 -0.4536 0.3852 -1.1776 0.2422
## ma7 0.4224 0.1494 2.8273 0.0058
## ma8 0.5803 0.2243 2.5872 0.0113
## sar1 0.9996 0.0022 454.3636 0.0000
## sma1 -0.9801 0.0593 -16.5278 0.0000printstatarima.Zt(modeld2)##
## Coefficients:
## s.e. t sign.
## ma1 0.3757 0.1048 3.5849 6e-04
## sar1 0.9994 0.0035 285.5429 0e+00
## sma1 -0.9720 0.0817 -11.8972 0e+00printstatarima.Zt(modele2)##
## Coefficients:
## s.e. t sign.
## ar1 0.3750 0.1088 3.4467 0.0009
## ar2 -0.1831 0.1072 -1.7080 0.0912
## sar1 0.9996 0.0024 416.5000 0.0000
## sma1 -0.9775 0.0653 -14.9694 0.0000Berdasarkan hasil di atas, dengan mempertimbangkan uji signifikansi penduga parameter dan ukuran keakuratan ramalan dapat disimpulkan bahwa dengan menggunakan model musiman (1,0,1)12 model terbaik adalah ARIMA model C
Model untuk komponen musiman ARIMA (1,0,2)12
modela3 <- Arima(train.ts,order=c(0,1,4),seasonal=list(order=c(1,0,2),period=12))
summary(modela3)## Series: train.ts
## ARIMA(0,1,4)(1,0,2)[12]
##
## Coefficients:
## ma1 ma2 ma3 ma4 sar1 sma1 sma2
## 0.3398 0.0347 -0.1807 -0.2830 0.9994 -1.1543 0.1840
## s.e. 0.1040 0.1106 0.1042 0.1029 0.0035 0.2284 0.2063
##
## sigma^2 = 16456: log likelihood = -554.31
## AIC=1124.62 AICc=1126.46 BIC=1144.34
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.7546349 122.3097 95.18991 -0.00285626 1.047509 0.7128477
## ACF1
## Training set 0.005553062#ARIMA(5,1,3)x(1,0,1)12
modelb3 <- Arima(train.ts,order=c(5,1,3),seasonal=list(order=c(1,0,1),period=12))
summary(modelb3)## Series: train.ts
## ARIMA(5,1,3)(1,0,1)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 ma2 ma3
## -0.0590 -0.3089 0.4069 -0.4669 0.0581 0.4473 0.3928 -0.5644
## s.e. 0.5162 0.1932 0.2426 0.0987 0.2726 0.6217 0.4465 0.3030
## sar1 sma1
## 0.9987 -0.9536
## s.e. 0.0085 0.1399
##
## sigma^2 = 15847: log likelihood = -549.66
## AIC=1121.32 AICc=1124.84 BIC=1148.44
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 1.155707 117.7543 94.3834 0.01732897 1.042683 0.706808 0.01425039tmodel <- Arima(train.ts, order=c(7,1,8),seasonal=c(1,0,2))
summary(tmodel)## Series: train.ts
## ARIMA(7,1,8)
##
## Coefficients:## Warning in sqrt(diag(x$var.coef)): NaNs produced## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ma1
## 0.6893 0.0214 -0.4828 -0.0071 0.0587 0.4618 -0.6756 -0.336
## s.e. NaN 0.1405 NaN 0.1265 0.1068 0.1233 0.0705 NaN
## ma2 ma3 ma4 ma5 ma6 ma7 ma8
## -0.3149 0.1561 0.0727 0.1744 -0.7132 0.4139 0.5982
## s.e. NaN NaN NaN NaN NaN NaN 0.1218
##
## sigma^2 = 19971: log likelihood = -552.07
## AIC=1136.14 AICc=1143.91 BIC=1175.6
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 23.11811 127.829 99.53925 0.2505431 1.099021 0.7454186 -0.02677237#ARIMA(7,1,8)x(1,0,1)12
modelc3 <- Arima(train.ts,order=c(7,1,8),seasonal=list(order=c(1,0,1),period=12))
summary(modelc3)## Series: train.ts
## ARIMA(7,1,8)(1,0,1)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7 ma1
## 0.6304 -0.1539 -0.0171 -0.2284 0.3283 0.2347 -0.4719 -0.3145
## s.e. 0.3098 0.4593 0.4599 0.3770 0.3928 0.4077 0.2431 0.3246
## ma2 ma3 ma4 ma5 ma6 ma7 ma8 sar1
## -0.0820 -0.1953 0.1221 -0.1179 -0.4536 0.4224 0.5803 0.9996
## s.e. 0.3789 0.3742 0.2980 0.3627 0.3852 0.1494 0.2243 0.0022
## sma1
## -0.9801
## s.e. 0.0593
##
## sigma^2 = 15314: log likelihood = -546.45
## AIC=1128.89 AICc=1138.95 BIC=1173.28
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3.771469 110.3702 87.65015 0.04371762 0.9696608 0.6563848
## ACF1
## Training set 0.003824925#ARIMA(0,1,1)x(1,0,1)12
modeld3 <-Arima(train.ts,order=c(0,1,1),seasonal=list(order=c(1,0,1),period=12))
summary(modeld3)## Series: train.ts
## ARIMA(0,1,1)(1,0,1)[12]
##
## Coefficients:
## ma1 sar1 sma1
## 0.3757 0.9994 -0.9720
## s.e. 0.1048 0.0035 0.0817
##
## sigma^2 = 18440: log likelihood = -558.85
## AIC=1125.7 AICc=1126.19 BIC=1135.57
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3.105948 132.6735 102.3998 0.03541785 1.127001 0.7668407
## ACF1
## Training set -0.0007157275#ARIMA(2,1,0)x(1,0,1)12
modele3 <- Arima(train.ts,order=c(2,1,0),seasonal=list(order=c(1,0,1),period=12))
summary(modele3)## Series: train.ts
## ARIMA(2,1,0)(1,0,1)[12]
##
## Coefficients:
## ar1 ar2 sar1 sma1
## 0.3750 -0.1831 0.9996 -0.9775
## s.e. 0.1088 0.1072 0.0024 0.0653
##
## sigma^2 = 18502: log likelihood = -558.34
## AIC=1126.68 AICc=1127.42 BIC=1139.01
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 3.210448 132.1029 100.8782 0.03617671 1.111796 0.7554458
## ACF1
## Training set -0.0209211printstatarima.Zt(modela3)##
## Coefficients:
## s.e. t sign.
## ma1 0.3398 0.1040 3.2673 0.0016
## ma2 0.0347 0.1106 0.3137 0.7545
## ma3 -0.1807 0.1042 -1.7342 0.0864
## ma4 -0.2830 0.1029 -2.7502 0.0072
## sar1 0.9994 0.0035 285.5429 0.0000
## sma1 -1.1543 0.2284 -5.0539 0.0000
## sma2 0.1840 0.2063 0.8919 0.3749printstatarima.Zt(modelb3)##
## Coefficients:
## s.e. t sign.
## ar1 -0.0590 0.5162 -0.1143 0.9093
## ar2 -0.3089 0.1932 -1.5989 0.1135
## ar3 0.4069 0.2426 1.6772 0.0971
## ar4 -0.4669 0.0987 -4.7305 0.0000
## ar5 0.0581 0.2726 0.2131 0.8317
## ma1 0.4473 0.6217 0.7195 0.4738
## ma2 0.3928 0.4465 0.8797 0.3814
## ma3 -0.5644 0.3030 -1.8627 0.0659
## sar1 0.9987 0.0085 117.4941 0.0000
## sma1 -0.9536 0.1399 -6.8163 0.0000printstatarima.Zt(modelc3)##
## Coefficients:
## s.e. t sign.
## ar1 0.6304 0.3098 2.0349 0.0449
## ar2 -0.1539 0.4593 -0.3351 0.7384
## ar3 -0.0171 0.4599 -0.0372 0.9704
## ar4 -0.2284 0.3770 -0.6058 0.5462
## ar5 0.3283 0.3928 0.8358 0.4056
## ar6 0.2347 0.4077 0.5757 0.5663
## ar7 -0.4719 0.2431 -1.9412 0.0555
## ma1 -0.3145 0.3246 -0.9689 0.3353
## ma2 -0.0820 0.3789 -0.2164 0.8292
## ma3 -0.1953 0.3742 -0.5219 0.6031
## ma4 0.1221 0.2980 0.4097 0.6830
## ma5 -0.1179 0.3627 -0.3251 0.7459
## ma6 -0.4536 0.3852 -1.1776 0.2422
## ma7 0.4224 0.1494 2.8273 0.0058
## ma8 0.5803 0.2243 2.5872 0.0113
## sar1 0.9996 0.0022 454.3636 0.0000
## sma1 -0.9801 0.0593 -16.5278 0.0000printstatarima.Zt(modeld3)##
## Coefficients:
## s.e. t sign.
## ma1 0.3757 0.1048 3.5849 6e-04
## sar1 0.9994 0.0035 285.5429 0e+00
## sma1 -0.9720 0.0817 -11.8972 0e+00printstatarima.Zt(modele3)##
## Coefficients:
## s.e. t sign.
## ar1 0.3750 0.1088 3.4467 0.0009
## ar2 -0.1831 0.1072 -1.7080 0.0912
## sar1 0.9996 0.0024 416.5000 0.0000
## sma1 -0.9775 0.0653 -14.9694 0.0000Berdasarkan hasil di atas, dengan mempertimbangkan uji signifikansi penduga parameter dan ukuran keakuratan ramalan dapat disimpulkan bahwa dengan menggunakan model musiman (1,0,2)12 model terbaik adalah ARIMA model C
Uji Signifikasi Model
lmtest::coeftest(modela)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.423213 0.103777 4.0781 4.54e-05 ***
## ma2 0.027496 0.107578 0.2556 0.7982680
## ma3 -0.279793 0.105699 -2.6471 0.0081190 **
## ma4 -0.349883 0.095096 -3.6793 0.0002339 ***
## sma1 0.277800 0.116517 2.3842 0.0171163 *
## sma2 0.073161 0.121431 0.6025 0.5468450
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modelb)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.18422 0.28988 0.6355 0.525098
## ar2 -0.24641 0.15666 -1.5729 0.115740
## ar3 0.43928 0.18197 2.4141 0.015776 *
## ar4 -0.48574 0.10366 -4.6860 2.786e-06 ***
## ar5 0.17504 0.16931 1.0338 0.301213
## ma1 0.25911 0.26484 0.9783 0.327904
## ma2 0.24526 0.24715 0.9924 0.321026
## ma3 -0.67498 0.23242 -2.9042 0.003682 **
## sma1 0.29345 0.12150 2.4152 0.015727 *
## sma2 0.12340 0.13082 0.9432 0.345559
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modelc)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.7362549 0.2579436 2.8543 0.004313 **
## ar2 -0.0723992 0.3184168 -0.2274 0.820134
## ar3 -0.3275657 0.3217352 -1.0181 0.308620
## ar4 -0.1691050 0.2773238 -0.6098 0.542011
## ar5 0.2922269 0.3294704 0.8870 0.375101
## ar6 0.2847810 0.3135484 0.9083 0.363745
## ar7 -0.5584299 0.1914346 -2.9171 0.003533 **
## ma1 -0.3755821 0.2787436 -1.3474 0.177848
## ma2 -0.2374302 0.2656439 -0.8938 0.371434
## ma3 0.0397402 0.2631888 0.1510 0.879980
## ma4 0.1947963 0.3166051 0.6153 0.538379
## ma5 0.0060212 0.2448104 0.0246 0.980378
## ma6 -0.5675911 0.2439870 -2.3263 0.020002 *
## ma7 0.3310655 0.2030680 1.6303 0.103034
## ma8 0.6113463 0.1443632 4.2348 2.288e-05 ***
## sma1 0.1667490 0.1444415 1.1544 0.248320
## sma2 0.0071020 0.1568562 0.0453 0.963886
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modeld)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.445869 0.096997 4.5967 4.292e-06 ***
## sma1 0.277215 0.110313 2.5130 0.01197 *
## sma2 0.088939 0.114741 0.7751 0.43827
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modele)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.469404 0.106394 4.4120 1.024e-05 ***
## ar2 -0.252273 0.104780 -2.4077 0.01606 *
## sma1 0.246685 0.112125 2.2001 0.02780 *
## sma2 0.085533 0.114916 0.7443 0.45669
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modela2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.3432814 0.1035476 3.3152 0.0009158 ***
## ma2 0.0173609 0.1071284 0.1621 0.8712610
## ma3 -0.1952835 0.1028278 -1.8991 0.0575472 .
## ma4 -0.3042935 0.0972600 -3.1287 0.0017561 **
## sar1 0.9996191 0.0023078 433.1421 < 2.2e-16 ***
## sma1 -0.9779300 0.0660428 -14.8075 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modelb2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.0589700 0.5162252 -0.1142 0.90905
## ar2 -0.3088773 0.1931600 -1.5991 0.10980
## ar3 0.4068936 0.2425937 1.6773 0.09349 .
## ar4 -0.4668936 0.0987004 -4.7304 2.241e-06 ***
## ar5 0.0580896 0.2725728 0.2131 0.83124
## ma1 0.4472785 0.6217367 0.7194 0.47189
## ma2 0.3927864 0.4465174 0.8797 0.37904
## ma3 -0.5643520 0.3030139 -1.8625 0.06254 .
## sar1 0.9986814 0.0085101 117.3524 < 2.2e-16 ***
## sma1 -0.9535675 0.1398792 -6.8171 9.291e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modelc2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.6304384 0.3097916 2.0350 0.041847 *
## ar2 -0.1539489 0.4593243 -0.3352 0.737501
## ar3 -0.0170652 0.4599415 -0.0371 0.970403
## ar4 -0.2284464 0.3770234 -0.6059 0.544567
## ar5 0.3282583 0.3928388 0.8356 0.403377
## ar6 0.2347275 0.4076786 0.5758 0.564773
## ar7 -0.4718809 0.2430594 -1.9414 0.052207 .
## ma1 -0.3144992 0.3245667 -0.9690 0.332554
## ma2 -0.0819744 0.3788595 -0.2164 0.828698
## ma3 -0.1952789 0.3741899 -0.5219 0.601760
## ma4 0.1220817 0.2980010 0.4097 0.682049
## ma5 -0.1178537 0.3626570 -0.3250 0.745202
## ma6 -0.4535501 0.3851725 -1.1775 0.238986
## ma7 0.4224193 0.1494440 2.8266 0.004704 **
## ma8 0.5802632 0.2242771 2.5873 0.009674 **
## sar1 0.9996339 0.0022074 452.8457 < 2.2e-16 ***
## sma1 -0.9800735 0.0593092 -16.5248 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modeld2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.3757266 0.1048430 3.5837 0.0003388 ***
## sar1 0.9994164 0.0034516 289.5520 < 2.2e-16 ***
## sma1 -0.9719659 0.0816741 -11.9005 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modele2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.3749838 0.1088248 3.4458 0.0005695 ***
## ar2 -0.1831115 0.1072048 -1.7081 0.0876266 .
## sar1 0.9995956 0.0023744 420.9864 < 2.2e-16 ***
## sma1 -0.9774620 0.0653309 -14.9617 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modela3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.3397531 0.1039854 3.2673 0.001086 **
## ma2 0.0347080 0.1106084 0.3138 0.753680
## ma3 -0.1807421 0.1041754 -1.7350 0.082745 .
## ma4 -0.2830164 0.1028793 -2.7510 0.005942 **
## sar1 0.9994195 0.0035369 282.5668 < 2.2e-16 ***
## sma1 -1.1543064 0.2284316 -5.0532 4.345e-07 ***
## sma2 0.1840136 0.2063055 0.8919 0.372421
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modelb3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.0589700 0.5162252 -0.1142 0.90905
## ar2 -0.3088773 0.1931600 -1.5991 0.10980
## ar3 0.4068936 0.2425937 1.6773 0.09349 .
## ar4 -0.4668936 0.0987004 -4.7304 2.241e-06 ***
## ar5 0.0580896 0.2725728 0.2131 0.83124
## ma1 0.4472785 0.6217367 0.7194 0.47189
## ma2 0.3927864 0.4465174 0.8797 0.37904
## ma3 -0.5643520 0.3030139 -1.8625 0.06254 .
## sar1 0.9986814 0.0085101 117.3524 < 2.2e-16 ***
## sma1 -0.9535675 0.1398792 -6.8171 9.291e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modelc3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.6304384 0.3097916 2.0350 0.041847 *
## ar2 -0.1539489 0.4593243 -0.3352 0.737501
## ar3 -0.0170652 0.4599415 -0.0371 0.970403
## ar4 -0.2284464 0.3770234 -0.6059 0.544567
## ar5 0.3282583 0.3928388 0.8356 0.403377
## ar6 0.2347275 0.4076786 0.5758 0.564773
## ar7 -0.4718809 0.2430594 -1.9414 0.052207 .
## ma1 -0.3144992 0.3245667 -0.9690 0.332554
## ma2 -0.0819744 0.3788595 -0.2164 0.828698
## ma3 -0.1952789 0.3741899 -0.5219 0.601760
## ma4 0.1220817 0.2980010 0.4097 0.682049
## ma5 -0.1178537 0.3626570 -0.3250 0.745202
## ma6 -0.4535501 0.3851725 -1.1775 0.238986
## ma7 0.4224193 0.1494440 2.8266 0.004704 **
## ma8 0.5802632 0.2242771 2.5873 0.009674 **
## sar1 0.9996339 0.0022074 452.8457 < 2.2e-16 ***
## sma1 -0.9800735 0.0593092 -16.5248 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modeld3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.3757266 0.1048430 3.5837 0.0003388 ***
## sar1 0.9994164 0.0034516 289.5520 < 2.2e-16 ***
## sma1 -0.9719659 0.0816741 -11.9005 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1lmtest::coeftest(modele3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.3749838 0.1088248 3.4458 0.0005695 ***
## ar2 -0.1831115 0.1072048 -1.7081 0.0876266 .
## sar1 0.9995956 0.0023744 420.9864 < 2.2e-16 ***
## sma1 -0.9774620 0.0653309 -14.9617 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1aic_model <- data.frame(
"Nama Model" = c("Model a", "Model b","Model c", "Model d","Model e", "Model a2","Model b2", "Model c2","Model d2","Model e2","Model a3","Model b3","Model c3","Model d3","Model e3"),
"AIC" = c(modela$aic,modelb$aic,modelc$aic,modeld$aic,modele$aic,modela2$aic,modelb2$aic,modelc2$aic,modeld2$aic,modele2$aic,modela3$aic,modelb3$aic,modelc3$aic,modeld3$aic,modele3$aic))
aic_model## Nama.Model AIC
## 1 Model a 1135.675
## 2 Model b 1135.929
## 3 Model c 1138.172
## 4 Model d 1140.463
## 5 Model e 1140.144
## 6 Model a2 1123.458
## 7 Model b2 1121.318
## 8 Model c2 1128.894
## 9 Model d2 1125.702
## 10 Model e2 1126.679
## 11 Model a3 1124.616
## 12 Model b3 1121.318
## 13 Model c3 1128.894
## 14 Model d3 1125.702
## 15 Model e3 1126.679dplyr::arrange(.data=aic_model, AIC)## Nama.Model AIC
## 1 Model b2 1121.318
## 2 Model b3 1121.318
## 3 Model a2 1123.458
## 4 Model a3 1124.616
## 5 Model d2 1125.702
## 6 Model d3 1125.702
## 7 Model e2 1126.679
## 8 Model e3 1126.679
## 9 Model c2 1128.894
## 10 Model c3 1128.894
## 11 Model a 1135.675
## 12 Model b 1135.929
## 13 Model c 1138.172
## 14 Model e 1140.144
## 15 Model d 1140.463Model tentatif dipilih berdasarkan nilai AIC minimum dengan pertimbangan bahwa seluruh parameter signifikan. Oleh karena itu, model A2 Arima(1,0,4)(1,0,1)12 yang menghasilkan hampir seluruh parameter signifikan serta nilai AIC yang cukup minimum akan dianalisis lebih lanjut sebagai model tentatif serta akan digunakan untuk melakukan peramalan terhadap data testing.
Analisis Sisaan
sisaan <- modela2$residuals
par(mfrow = c(2,2))
qqnorm(sisaan)
qqline(sisaan, col = "blue", lwd = 2)
plot(sisaan, type="o",
ylab = "Sisaan", xlab = "Order", main = "Sisaan M6d vs Order")
abline(h = 0, col='red')
acf(sisaan)
pacf(sisaan)
ks.test(sisaan,"pnorm") ##
## Exact one-sample Kolmogorov-Smirnov test
##
## data: sisaan
## D = 0.5, p-value < 2.2e-16
## alternative hypothesis: two-sidedBerdasarkan hasil uji di atas, didapat nilai P−Value<0.0001 yang berarti TOLAK H0. Artinya, pada taraf nyata 5%, ada bukti untuk menyatakan bahwa sisaan tidak mengikuti sebaran normal.
Box.test(sisaan, type = "Ljung") ##
## Box-Ljung test
##
## data: sisaan
## X-squared = 0.0029871, df = 1, p-value = 0.9564Berdasarkan hasil uji di atas, didapat nilai P−Value=0.9564 yang berarti TERIMA H0. Artinya, pada taraf nyata 5%, ada bukti untuk menyatakan bahwa tidak ada autokorelasi pada data.
t.test(sisaan, mu = 0, conf.level = 0.95) ##
## One Sample t-test
##
## data: sisaan
## t = 0.27044, df = 87, p-value = 0.7875
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -23.24717 30.56976
## sample estimates:
## mean of x
## 3.661292Berdasarkan hasil uji di atas, didapat nilai P−Value=0.7447 yang berarti TERIMA H0. Artinya, pada taraf nyata 5%, ada bukti untuk menyatakan bahwa nilai tengah sisaan bernilai nol
Diagnostik Model
#diagnostik model
checkresiduals(modela2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,4)(1,0,1)[12]
## Q* = 6.639, df = 4, p-value = 0.1562
##
## Model df: 6. Total lags used: 10tsdisplay(residuals(modela2), lag.max=50, main='ARIMA C Model Residuals')
library(tseries)
jarque.bera.test(residuals(modela2))##
## Jarque Bera Test
##
## data: residuals(modela2)
## X-squared = 7.0571, df = 2, p-value = 0.02935Selanjutnya, apabila dilihat melalui histogram sisaan, terlihat bahwa sisaaan menyebar normal, selain itu, berdasarkan pengujian menggunakan Jarque bera Test diperoleh p-value = 0.02935 > 0.05 yang berarti bahwa H0 ditolak sehingga dapat disimpulkan sisaan model tidak menyebar normal.
#validasi model
ramalan_arima3 = forecast(modela2, 120)
ramalan_arima3## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 89 10054.76 9878.190 10231.33 9784.721 10324.79
## 90 10082.62 9786.983 10378.26 9630.482 10534.76
## 91 10058.48 9677.557 10439.40 9475.910 10641.04
## 92 10048.70 9615.781 10481.63 9386.605 10710.80
## 93 10076.72 9617.894 10535.55 9375.004 10778.44
## 94 10124.36 9640.955 10607.76 9385.057 10863.66
## 95 10209.37 9702.546 10716.19 9434.251 10984.48
## 96 10343.94 9814.773 10873.11 9534.648 11153.24
## 97 10576.98 10026.276 11127.68 9734.753 11419.20
## 98 10592.23 10021.162 11163.30 9718.855 11465.61
## 99 10441.99 9851.254 11032.73 9538.535 11345.45
## 100 10215.49 9605.714 10825.26 9282.919 11148.06
## 101 10215.02 9581.987 10848.05 9246.881 11183.15
## 102 10275.90 9618.746 10933.06 9270.868 11280.94
## 103 10312.02 9631.505 10992.54 9271.260 11352.79
## 104 10333.44 9631.260 11035.63 9259.546 11407.34
## 105 10361.45 9639.604 11083.30 9257.481 11465.42
## 106 10409.07 9667.917 11150.22 9275.575 11542.56
## 107 10494.04 9733.976 11254.11 9331.620 11656.47
## 108 10628.57 9850.076 11407.06 9437.966 11819.17
## 109 10861.51 10065.084 11657.94 9643.479 12079.55
## 110 10876.76 10063.406 11690.12 9632.840 12120.69
## 111 10726.58 9896.641 11556.52 9457.296 11995.87
## 112 10500.16 9653.963 11346.36 9206.012 11794.31
## 113 10499.69 9633.683 11365.70 9175.245 11824.14
## 114 10560.56 9673.809 11447.30 9204.394 11916.72
## 115 10596.66 9689.577 11503.75 9209.395 11983.93
## 116 10618.07 9691.864 11544.29 9201.557 12034.59
## 117 10646.07 9702.258 11589.89 9202.633 12089.51
## 118 10693.67 9732.355 11654.99 9223.465 12163.87
## 119 10778.61 9799.966 11757.26 9281.902 12275.32
## 120 10913.09 9917.443 11908.73 9390.380 12435.79
## 121 11145.94 10133.733 12158.15 9597.901 12693.99
## 122 11161.19 10133.454 12188.92 9589.405 12732.97
## 123 11011.06 9968.037 12054.09 9415.893 12606.23
## 124 10784.73 9726.633 11842.83 9166.510 12402.95
## 125 10784.26 9707.914 11860.60 9138.132 12430.39
## 126 10845.10 9749.610 11940.59 9169.692 12520.51
## 127 10881.19 9766.818 11995.57 9176.904 12585.48
## 128 10902.60 9770.335 12034.86 9170.953 12634.24
## 129 10930.58 9781.730 12079.44 9173.564 12687.60
## 130 10978.16 9812.693 12143.63 9195.729 12760.60
## 131 11063.07 9881.049 12245.10 9255.322 12870.83
## 132 11197.50 9999.179 12395.82 9364.827 13030.17
## 133 11430.26 10216.075 12644.45 9573.321 13287.21
## 134 11445.50 10216.538 12674.47 9565.964 13325.04
## 135 11295.44 10051.871 12539.00 9393.568 13197.30
## 136 11069.19 9811.194 12327.18 9145.252 12993.12
## 137 11068.72 9793.358 12344.08 9118.223 13019.21
## 138 11129.54 9835.947 12423.12 9151.162 13107.91
## 139 11165.61 9853.987 12477.24 9159.654 13171.57
## 140 11187.01 9858.219 12515.80 9154.799 13219.22
## 141 11214.99 9870.189 12559.78 9158.296 13271.68
## 142 11262.55 9901.640 12623.46 9181.219 13343.88
## 143 11347.43 9970.406 12724.45 9241.455 13453.40
## 144 11481.80 10088.885 12874.71 9351.521 13612.08
## 145 11714.48 10306.095 13122.86 9560.542 13868.41
## 146 11729.71 10307.007 13152.41 9553.873 13905.55
## 147 11579.70 10142.818 13016.58 9382.179 13777.22
## 148 11353.54 9902.617 12804.46 9134.546 13572.53
## 149 11353.07 9885.340 12820.80 9108.371 13597.77
## 150 11413.86 9928.495 12899.23 9142.189 13685.54
## 151 11449.93 9947.068 12952.79 9151.502 13748.35
## 152 11471.32 9951.756 12990.88 9147.350 13795.28
## 153 11499.28 9964.089 13034.47 9151.407 13847.16
## 154 11546.83 9995.841 13097.81 9174.800 13918.85
## 155 11631.67 10064.848 13198.49 9235.422 14027.92
## 156 11765.99 10183.522 13348.46 9345.813 14186.17
## 157 11998.58 10400.896 13596.27 9555.132 14442.03
## 158 12013.81 10402.100 13625.52 9548.913 14478.71
## 159 11863.86 10238.244 13489.47 9377.697 14350.01
## 160 11637.78 9998.385 13277.18 9130.541 14145.02
## 161 11637.31 9981.488 13293.13 9104.949 14169.67
## 162 11698.08 10025.026 13371.14 9139.364 14256.80
## 163 11734.13 10043.963 13424.30 9149.240 14319.03
## 164 11755.51 10048.963 13462.06 9145.569 14365.46
## 165 11783.47 10061.537 13505.40 9150.002 14416.93
## 166 11830.99 10093.484 13568.50 9173.702 14488.29
## 167 11915.81 10162.637 13668.98 9234.565 14597.05
## 168 12050.08 10281.419 13818.73 9345.148 14755.00
## 169 12282.58 10498.870 14066.29 9554.633 15010.52
## 170 12297.80 10500.274 14095.32 9548.722 15046.88
## 171 12147.90 10336.666 13959.14 9377.855 14917.95
## 172 11921.91 10097.068 13746.76 9131.053 14712.78
## 173 11921.44 10080.444 13762.45 9105.877 14737.01
## 174 11982.19 10124.253 13840.13 9140.719 14823.67
## 175 12018.23 10143.448 13893.01 9150.998 14885.46
## 176 12039.60 10148.669 13930.53 9147.670 14931.53
## 177 12067.55 10161.412 13973.68 9152.366 14982.73
## 178 12115.05 10193.487 14036.62 9176.271 15053.84
## 179 12199.84 10262.729 14136.94 9237.287 15162.38
## 180 12334.05 10381.566 14286.54 9347.981 15320.13
## 181 12566.47 10599.042 14533.89 9557.551 15575.38
## 182 12581.68 10600.588 14562.78 9551.860 15611.50
## 183 12431.84 10437.172 14426.51 9381.257 15482.43
## 184 12205.94 10197.785 14214.10 9134.732 15277.15
## 185 12205.47 10181.363 14229.58 9109.865 15301.08
## 186 12266.20 10225.369 14307.02 9145.020 15387.37
## 187 12302.22 10244.754 14359.68 9155.598 15448.84
## 188 12323.58 10250.137 14397.03 9152.522 15494.64
## 189 12351.52 10263.000 14440.03 9157.406 15545.63
## 190 12399.01 10295.161 14502.85 9181.453 15616.56
## 191 12483.75 10364.454 14603.06 9242.564 15724.95
## 192 12617.92 10483.312 14752.53 9353.318 15882.53
## 193 12850.25 10700.779 14999.71 9562.919 16137.57
## 194 12865.46 10702.427 15028.49 9557.388 16173.52
## 195 12715.67 10539.167 14892.18 9386.994 16044.35
## 196 12489.86 10299.957 14679.76 9140.693 15839.02
## 197 12489.39 10283.689 14695.09 9116.062 15862.71
## 198 12550.09 10327.840 14772.34 9151.453 15948.73
## 199 12586.10 10347.369 14824.83 9162.257 16009.94
## 200 12607.46 10352.874 14862.04 9159.370 16055.54
## 201 12635.38 10365.822 14904.93 9164.392 16106.36
## 202 12682.85 10398.042 14967.66 9188.537 16177.16
## 203 12767.57 10467.360 15067.77 9249.704 16285.43
## 204 12901.68 10586.216 15217.15 9360.482 16442.88
## 205 13133.92 10803.649 15464.19 9570.079 16697.76
## 206 13149.12 10805.373 15492.87 9564.667 16733.58
## 207 12999.40 10642.244 15356.55 9394.443 16604.35
## 208 12773.67 10403.186 15144.15 9148.329 16399.00plot(ramalan_arima3)
Berdasarkan ukuran keakuratan ramalan di atas, dapat disimpulkan bahwa model ARIMA()x()12 dapat digunakan untuk meramalkan data harga beras premium dengan baik. Selain itu, apabila dilihat berdasarkan plot ramalan, nilai ramalan yang dihasilkan memiliki nilai yang .. dengan nilai data sebelumnya.
Peramalan (5 bulan ke depan)
Berikut merupakan hasil ramalan data harga beras premium untuk 5 bulan ke depan:
modelZt<-Arima(data$Premium,order=c(7,1,8),seasonal=c(1,0,2))
ramalan_arimaZt = forecast(modelZt, 5)
ramalan_arimaZt## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 116 9726.153 9554.557 9897.749 9463.719 9988.587
## 117 9809.721 9514.223 10105.218 9357.796 10261.645
## 118 9907.472 9524.601 10290.343 9321.921 10493.023
## 119 9930.126 9506.461 10353.792 9282.186 10578.067
## 120 9903.078 9463.719 10342.437 9231.137 10575.020plot(ramalan_arimaZt)
Daftar Pustaka
Fani E, Widjajati FA, Soehardjoepri. 2017. Perbandingan Metode Winter Eksponensial Smoothing dan Metode Event Based untuk Menentukan Penjualan Produk Terbaik di Perusahaan X. Jurnal Sains dan Seni ITS. 6(1):2337-3529.
H. R. Naufal and R Adrean, “R. naufal Hayâ and R. Adrean, ‘Sistem Informasi Inventory Berdasarkan Prediksi Data Penjualan Barang Menggunakan Metode Single Moving Average Pada CV. Agung Youanda,’ ProTekInfo (Pengembangan Ris. dan Obs.Tek. Inform., vol. 4, pp. 29??”33, 2017,” ProTekInfo (Pengembangan Ris. dan Obs. Tek.Inf., vol. 4, pp. 29??“33, 2017.
Hartono A, Dwijana D, Handiwidjojo. 2012. Perbandingan metode single exponential smoothing dan exponential smoothing adjusted for trend untuk meramalkan penjualan. Studi kasus: Toko onderdil mobil “Prodi, Purwodadi”. Jurnal EKSIS. 5(1):8-18.
Hapsari, V. 2013. Perbandingan Metode Dekomposisi Klasik dengan Metode Pemulusan Eksponensial Holt-Winter dalam meramalkan Tingkat Pencemaran Udara di Kota Bandung Periode 2003-2012. Universitas Pendidikan Indonesia.
Muzaki L. 2022. Contoh studi kasus metode double moving average, kelebihan kekurangan, dan cara menghitung [internet]. [diacu pada 2022 Oktober 25]. Tersedia dari: https://www.pengadaanbarang.co.id/2022/05/metode-double-moving-average.html.
M. Layakana, S. Iskandar. “Penerapan Metode Double Moving Average dan Double Eksponensial Smoothing Dalam Meramalkan Jumlah Produksi Crude Palm Oil(CPO) Pada PT.Perkebunan Nusantara IV Unit Dolok Sinumbah. 6(1): 47.
Purwanto A, Hanief S. 2017. Teknik peramalan dengan double exponential smoothing pada distributor gula. Jurnal Teknologi Informasi dan Komputer. 3(1):362-367