Anggota Kelompok :
Hanung Safrizal (G1401201050)
Maulana Ahsan Fadillah (G1401201062)
Muhammad Irsyad Robbani (G1401201064)
Muhammad raiv Zulfikar (G1401201083)
Dhea Adinda Puspita (G1401201090)Pendahuluan
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). Maka dari itu, pada penugasan kali ini kelompok kami bertujuan untuk membahas mengenai kualitas serta harga dari produk beras dari waktu ke waktu.
Tujuan
Pada penugasan kali ini kelompok kami bertujuan untuk membahas mengenai kualitas serta harga dari produk beras dari waktu ke waktu.
Tinjauan Pustaka
Tinjauan Pustaka dari penelitian ini merupakan hasil dari penelitian-penelitian terdahulu.
Sukiyono, K & Rosdiana (2018) dalam penelitiannya mengenai Pendugaan Model Peramalan Harga Beras pada Tingkat Grosir. Berdasarkan hasil dari peneilitan tersebut, disimpulkan bahwa model MA(2) adalah model terbaik jika digunakan model moving average. Akan tetapi, jika menggunakan model single exponential, model dengan α = 0.9 merupakan model terbaik. Pemilihan model tersebut didasarkan pada kriteria MSE, MAPE, dan MAD terkecil. Model Single Exponential α = 0.9 merupakan model terbaik jika dibandingkan dengan model lainnya. Menurut peneliti, Dibutuhkan kehati-hatian dalam pemilihan model peramalan terbaik, khususnya jika jumlah data yang digunakan juga berubah.
Listiowarni et al (2020) dalam penelitiannya mengenai Perbandingan Metode Double Exponential Smoothing dan Doube Moving Average untuk Peramalan Harga Beras Eceran di Kabupaten Pamekasan. Berdasarkan hasil penelitian tersebut yaitu menggunakan metode Double Moving Average dan Double Exponential Smoothing, didapatkan bahwa metode Double Moving Average lebih baik dibandingkan dengan metode Double Exponential Smoothing. Hal ini terlihat pada perolehan nilai MSE dan MAPE dari metode Double Moving Average yang lebih kecil. Semakin tinggi nilai α yang digunakan pada metode Double Exponential Smoothing kasus ini, semakin besar nilai error yang diperoleh, yang berakibat pada nilai MAPE dan MSE yang semakin besar juga. Sedangkan Pengujuan periode (orde waktu) pada metode Double Moving Average mengalami peningkatan seiring bertambah besarnya nilai periode.
Nugraheni et al (2022) dalam penelitiannya mengenai Penerapan Meotde Exponential Smoothin Winters pada Prediksi Harga Beras. Berdasarkan hasil penelitian tersebut didapat bahwa perhitungan prediksi harga beras di Kabupaten Sukoharjo menggunakan Exponential Smoothing Winters sudah memberikan hasil sesuai yang diharapkan dengan nilai MAPE sebesar 3.91% yang termasuk ke dalam kategori < 10 dimana hasil ini memiliki peramalan yang baik dengan α = 0.4, β = 0.1, dan μ = 0.3 untuk harga beras premium dan MAPE sebesar 4.24% yang termasuk ke dalam kategori < 10 dimana hasil ini juga memiliki peramalan yang baik dengan α = 0.4, β = 0.1, dan μ = 0.3.
Kestasioneran Deret Waktu
Stasioneritas berarti bahwa tidak terdapat perubahan yang drastis pada data. Fluktuasi data berada disekitar suatu nilai rata-rata yang konstan, tidak tergantung pada waktu dan variansi dari fluktuasi tersebut. Data time series dikatakan stasioner jika rata-rata dan variansinya konstan, tidak ada unsur trend dalam data, dan tidak ada unsur musiman (Deviana et.al, 2021)
Uji Stasioneritas Data Deret Waktu
Apabila data tidak stasioner, maka sebelum mencari model terbaik, data yang ada perlu distasionerkan terlebih dahulu. Apabila data yang digunakan dalam model ada yang tidak stasioner, maka ada kemungkinan terjadinya spurious regression. Spurious regression adalah regresi yang memiliki \(R^2\) yang tinggi, tetapi tidak mempunyai hubungan yang berarti. Secara umum, uji untuk mengetahui kestasioneran data time series dapat dikategorikan menjadi tiga, yakni melalui grafik, korelogram, dan uji unit root (ADF Test & Phillips-Perron Test).
1. Uji Stasioneritas Menggunakan Grafik
Uji stasioneritas data time series menggunakan grafik dapat dilakukan dengan membuat plot antara data observasi dengan variabel waktu (t). Jika dari plot tersebut, terlihat rata-rata dan variansnya konstan, maka data time series dikatakan stasioner. Sebaliknya, jika grafik tidak menunjukkan rata-rata dan varians konstan, maka data time series dikatakan tidak stasioner.
2. Uji Stasioneritas Menggunakan Korelogram (Correlogram)
Pada dasarnya, korelogram merupakan metode pengujian stasioneritas data time series berdasarkan fungsi autokorelasi (ACF) yang diperoleh dengan memplotkan antara \(\rho_k\) dan \(k\)(lag). Untuk data yang stasioner, korelogram menurun dengan cepat seiring dengan meningkatnya \(k\). Sedangkan untuk data yang tidak stasioner, korelogram cenderung tidak menuju nol (turun lambat).
3. Uji Stasioneritas Menggunakan Uji Unit Root
Uji unit root (uji akar unit) merupakan uji untuk mengetahui stasioneritas data time series yang sering digunakan. Uji ini dikembangkan oleh David Dickey dan Wayne Fuller sehingga dikenal dengan sebutan Augmented Dickey-Fuller Test (ADF Test). Uji unit root lain yang juga sering digunakan yaitu Uji Phillips-Perron.
Model Time Series Stasioner
Model Autoregressive (AR) adalah model regresi time series yang menghubungkan nilai pengamatan aktual dengan nilai pengamatan sebelumnya. Ini dapat dilakukan ketika sebuah pengamatan tidak lepas dari pengamatan yang terjadi sebelumnya. Konsep dasar model ini yaitu meregresikan pengamatan aktual dengan nilai pengamatan sebelumnya untuk melakukan peramalan nilai ke depan. Plot PACF dapat digunakan untuk mengidentifikasi orde dari model AR(p). Adapun persamaan model autoregressive dengan order p atau AR(p) dapat dituliskan sebagai berikut: \[ yt=\phi +\phi_1Y_{t-1}+\phi_2Y_{t-2}+...+\phi_pY_{t-p}+\epsilon_t \] dimana \(\epsilon_t \sim WN(0,\sigma^2)\)
Pada model moving average (MA) digunakan nilai error dalam model regresi. Adapun persamaan model moving average dengan orde q atau MA(q) dapat dituliskan sebagai: \[ y_t=\mu +wt+\theta_1w_{t-1}+\theta_2w_{t-2}+...+\theta_pw_{t-p}+\epsilon_t \]
Model ARMA merupakan campuran antara model autoregressive (AR) yang mengasumsikan bahwa data sekarang dipengaruhi oleh data sebelumnya dan model moving average (MA) yang mengasumsikan bahwa data sekarang dipengaruhi oleh nilai residual data sebelumnya. Secara umum, model ARMA dengan orde p dan q atau ARMA(p,q) diberikan sebagai berikut: \[ yt=\gamma +\phi_1Y_{t-1}+\phi_2Y_{t-2}+...+\phi_pY_{t-p}+\epsilon_t-\theta_1e_{t-1}-\theta_2e_{t-2}-...-\theta_qe_{t-q} \]
Model Time Series Tak Stasioner
Differencing (pembedaan) digunakan untuk menstasionerkan data yang tidak stasioner. Jika pembedaan pertama (first difference) berhasil membuat data menjadi stasioner, berarti kita peroleh orde d = 1 untuk ARIMA. Selanjutkan kita menentukan orde p dan q berdasarkan plot ACF dan PACF sampel dari data yang telah distasionerkan tersebut. Namun, ada kalanya pembedaan pertama belum berhasil membuat data menjadi stasioner. Jika kondisi demikian terjadi maka kita dapat lakukan pembedaan kedua (second difference), dan jika datanya telah stasioner berarti kita peroleh d = 2, dan selanjutnya menentukan p dan q berdasarkan plot ACF dan PACF sampel dari data yang telah stasioner. Perlu diketahui bahwa untuk diferensiasi atau pembedaan yang lebih tinggi tidaklah disarankan. Umumnya, paling banyak kita hanya melakukan pembedaan sebanyak dua kali. Ini dikarenakan semakin tinggi pembedaan yang dilakukan, maka model akan semakin kompleks. Jika pembedaan kedua tidak berhasil membuat data menjadi stasioner, mungkin ada baiknya mempertimbangkan untuk dilakukan transformasi data tersebut dahulu.
Model ARIMA atau dikenal juga dengan model Box-Jenkins menghadirkan kelas model yang sangat kuat dan fleksibel untuk analisis dan peramalan data deret waktu. Selama bertahun-tahun, model tersebut telah dengan sangat berhasil diterapkan pada banyak masalah dalam penelitian dan praktik. Namun, mungkin ada situasi tertentu di mana model tersebut mungkin gagal memberikan jawaban yang “benar”. Misalnya, dalam model ARIMA, meramalkan pengamatan di masa depan sangat bergantung pada data masa lalu dan secara implisit mengasumsikan bahwa kondisi di mana data dikumpulkan akan tetap sama di masa mendatang. Dalam banyak situasi, asumsi ini mungkin (dan kemungkinan besar akan) tidak sesuai. Untuk kasus tersebut, fungsi transfer-model noise, di mana sekumpulan variabel input yang mungkin memiliki efek pada deret waktu ditambahkan ke model, menyediakan opsi yang sesuai. Model ARIMA(p,d,q) dapat dituliskan sebagai: \[ yt =\gamma_0+\theta_1Y_{t-1}+\theta_2Y_{t-2}+\theta_nY_{t-q}-\lambda_1e_{t-1}-\lambda_2e_{t-2}-...-\lambda_ne_{t-q} \]
Input 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_excel("C:/R/Data MPDW.xlsx", sheet = "Sheet1")
data <- data[,3]
data2 <- read_excel("C:/R/Harga Gabah.xlsx")
data2 <- data2[,3]
data <- cbind(data,data2)
data$Premium <- as.numeric(data$Premium)
data$`Harga Gabah` <- as.numeric(data$`Harga Gabah`)
y <- data$Premium
x <- data$`Harga Gabah`Train-Test Split
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. Data dibagi dengan perbandingan 80% : 20%. Jumlah data training sebanyak 92 amatan (amatan 1 sampai 92) dan data testing sebanyak 23 amatan (amatan 93 sampai 115).
train <- data[1:92, 1]
test <- data[93:115, 1]
# Time-Series Object
data.ts <- ts(data, start = 2013, frequency = 12)
train.ts <- ts(train, start = 2013, frequency = 12)
test.ts <- ts(test)Eksplorasi Data
Plot time-series
plot(data.ts, col = "blue", type = "o", lwd = 1.5, main = "Time Series Plot Harga Beras dan Gabah 2013-2022",
cex = 0.7)
plot(data$`Harga Gabah`, data$Premium, pch = 20, col = "blue", main = "Harga Beras Premium vs Harga Gabah Kering")
Berdasarkan scatter plot di atas, dapat dilihat bahwa hubungan peubah harga beras premium (X) dan harga gabah kering (Y) memiliki hubungan yang positif cukup kuat yaitu dengan nilai korelasi sebesar 0.817009, artinya dengan kenaikan harga beras premium, harga gabah juga naik.
Pemulusan
Moving Average
# Parameter Optimum
#Fungsi Mencari Parameter Optimum
pemulusan <- function(x,min,max,method=c("SMA","DMA"),
metrik=c("SSE","MSE","RMSE","MAD","MAPE")){
df.master <- data.frame()
if(method=="SMA"){
for(i in min:max) {
sma <- TTR::SMA(x,i)
ramal <- c(NA,sma)
df <- cbind(aktual=c(x,NA),mulus=c(sma,NA),ramal)
error <- df[,1]-df[,3]
sse <- sum(error^2,na.rm=T)
mse <- mean(error^2,na.rm=T)
rmse <- sqrt(mse)
mad <- mean(abs(error),na.rm=T)
rerror <- error/df[,1]*100
mape <- mean(abs(rerror),na.rm=T)
ak <- data.frame(n=i,SSE=sse,MSE=mse,RMSE=rmse,MAD=mad,MAPE=mape)
df.master <- rbind(df.master,ak)
}
}
else if(method=="DMA") {
for(i in min:max) {
df_ts_sma <- TTR::SMA(x, n=i)
df_ts_dma <- TTR::SMA(df_ts_sma, n=i)
At <- 2*df_ts_sma-df_ts_dma
Bt <- df_ts_sma-df_ts_dma
pemulusan_dma <- At+Bt
ramal_dma <- c(NA, pemulusan_dma)
df_dma <- cbind(df_aktual=c(x,NA), pemulusan_dma=c(pemulusan_dma,NA), ramal_dma)
error.dma <- df_dma[, 1] - df_dma[, 3]
SSE.dma <- sum(error.dma^2, na.rm = T)
MSE.dma <- mean(error.dma^2, na.rm = T)
RMSE.dma <- sqrt(mean(error.dma^2, na.rm = T))
MAD.dma <- mean(abs(error.dma), na.rm = T)
r.error.dma <- (error.dma/df_dma[, 1])*100
MAPE.dma <- mean(abs(r.error.dma), na.rm = T)
ak <- data.frame(n=i,SSE=SSE.dma,MSE=MSE.dma,RMSE=RMSE.dma,
MAD=MAD.dma,MAPE=MAPE.dma)
df.master <- rbind(df.master,ak)
}
}
opt <- df.master$n[which.min(df.master[,metrik])]
return(list(Akurasi=df.master,n.optimum=paste("n optimum adalah",opt)))
}pemulusan(train.ts, min = 2, max = 10, method = "DMA", metrik = "RMSE")## $Akurasi
## n SSE MSE RMSE MAD MAPE
## 1 2 4346898 48841.55 221.0012 159.7996 1.740322
## 2 3 8846711 101686.34 318.8830 223.7614 2.429434
## 3 4 12622862 148504.26 385.3625 272.8248 2.938180
## 4 5 14936523 179958.11 424.2147 306.9679 3.286957
## 5 6 16599084 204926.97 452.6886 340.5738 3.636755
## 6 7 17606585 222868.16 472.0891 369.2973 3.936377
## 7 8 17268000 224259.74 473.5607 379.6653 4.041152
## 8 9 15802163 210695.51 459.0158 371.7431 3.936761
## 9 10 14283373 195662.64 442.3377 355.5422 3.759233
##
## $n.optimum
## [1] "n optimum adalah 2"Dari metode Moving Average n optimum adalah 2
# SMA
df_ts_sma <- SMA(train.ts, n=2)
# DMA
df_ts_dma <- SMA(df_ts_sma, n=2)
At <- 2*df_ts_sma - df_ts_dma
Bt <- df_ts_sma - df_ts_dma
pemulusan_dma <- At + Bt
ramal_dma <- c(NA, pemulusan_dma)
df_dma <- cbind(df_aktual = c(train.ts,NA), pemulusan_dma = c(pemulusan_dma,NA), ramal_dma)
head(df_dma)## df_aktual pemulusan_dma ramal_dma
## [1,] 7797.63 NA NA
## [2,] 7773.26 NA NA
## [3,] 7576.27 7564.085 NA
## [4,] 7420.72 7322.225 7564.085
## [5,] 7545.40 7467.625 7322.225
## [6,] 7548.22 7610.560 7467.625# Akurasi DMA
error.dma <- df_dma[, 1] - df_dma[, 3]
sse.dma <- sum(error.dma^2, na.rm = T)
mse.dma <- mean(error.dma^2, na.rm = T)
rmse.dma <- sqrt(mean(error.dma^2, na.rm = T))
mad.dma <- mean(abs(error.dma), na.rm = T)
r.error.dma <- (error.dma/df_dma[, 1])*100
mape.dma <- mean(abs(r.error.dma), na.rm = T)#Plot
ts.plot(train.ts, xlab = "periode waktu", ylab="Yt", col="blue", lty=3)
points(train.ts)
lines(pemulusan_dma , col = "red",lwd = 1.5)
lines(ramal_dma,col = "black",lwd = 2)
title("Rataan Bergerak Berganda N=2", cex.main = 1, font.main=4 ,col.main = "black")
legend("topleft", c("Data Aktual","Pemulusan","Ramalan"),lty=1:3,col=c ("blue","red","black"))
Exponential Smoothing
# Mencari Parameter Optimum
pemulusan2 <- function(x,alfa=NULL,beta=NULL,method=c("SES","DES"),
metrik=c("SSE","MSE","RMSE")){
df.master <- data.frame()
if(method=="SES"){
for(i in alfa){
df_ses <- HoltWinters(x, alpha = i, beta=F, gamma=F)
sse <- df_ses$SSE
mse <- sse/length(x)
rmse <-sqrt(mse)
ak <- data.frame(Alpha=i,SSE=sse,MSE=mse,RMSE=rmse)
df.master <- rbind(df.master,ak)
datases <- data.frame(x, c(NA, df_ses$fitted[,1]))
colnames(datases) = c("y","yhat")
}
}
else if(method=="DES"){
for (i in alfa) {
for (j in beta) {
df_des <- HoltWinters(x, alpha = i, beta=j, gamma=F)
sse <- df_des$SSE
mse <- sse/length(x)
rmse <-sqrt(mse)
ak <- data.frame(Alpha=i,Beta=j,SSE=sse,MSE=mse,RMSE=rmse)
df.master <- rbind(df.master,ak)
datades <- data.frame(x, c(NA, NA, df_des$fitted[,1]))
colnames(datades) = c("y","yhat")
}
}
}
if(method=="SES"){
opt <- df.master$Alpha[which.min(df.master[,metrik])]
return(list(Akurasi=df.master,n.optimum=paste("Alpha optimum adalah",opt)))
}
else if(method=="DES"){
opt <- df.master$Alpha[which.min(df.master[,metrik])]
opt2 <- df.master$Beta[which.min(df.master[,metrik])]
return(list(Akurasi=df.master,n.optimum=paste("Alpha optimum adalah",opt,
",Beta optimum adalah",opt2)))
}
}# Double Exponential Smoothing
pemulusan2(train.ts,alfa=c(seq(0.1, 0.9, by = 0.1)),
beta=c(seq(0.1, 0.9, by = 0.1)),method = "DES",metrik = "RMSE")## $Akurasi
## Alpha Beta SSE MSE RMSE
## 1 0.1 0.1 14725530 160060.10 400.0751
## 2 0.1 0.2 13869511 150755.56 388.2725
## 3 0.1 0.3 13396774 145617.11 381.5981
## 4 0.1 0.4 12441012 135228.39 367.7341
## 5 0.1 0.5 11565893 125716.23 354.5648
## 6 0.1 0.6 11229872 122063.82 349.3763
## 7 0.1 0.7 11364004 123521.78 351.4567
## 8 0.1 0.8 11860092 128914.04 359.0460
## 9 0.1 0.9 12704474 138092.11 371.6075
## 10 0.2 0.1 9173314 99709.93 315.7688
## 11 0.2 0.2 9207195 100078.20 316.3514
## 12 0.2 0.3 9563410 103950.11 322.4130
## 13 0.2 0.4 10391176 112947.56 336.0767
## 14 0.2 0.5 11679513 126951.22 356.3022
## 15 0.2 0.6 13355266 145165.94 381.0065
## 16 0.2 0.7 15309467 166407.25 407.9305
## 17 0.2 0.8 17426170 189414.90 435.2182
## 18 0.2 0.9 19584436 212874.31 461.3830
## 19 0.3 0.1 7650006 83152.23 288.3613
## 20 0.3 0.2 8217346 89318.98 298.8628
## 21 0.3 0.3 9150939 99466.73 315.3835
## 22 0.3 0.4 10356913 112575.14 335.5222
## 23 0.3 0.5 11659922 126738.28 356.0032
## 24 0.3 0.6 12855795 139736.90 373.8140
## 25 0.3 0.7 13751032 149467.74 386.6106
## 26 0.3 0.8 14231123 154686.12 393.3016
## 27 0.3 0.9 14327482 155733.50 394.6308
## 28 0.4 0.1 6738970 73249.67 270.6468
## 29 0.4 0.2 7414360 80590.87 283.8853
## 30 0.4 0.3 8257797 89758.66 299.5975
## 31 0.4 0.4 9087879 98781.29 314.2949
## 32 0.4 0.5 9738273 105850.79 325.3472
## 33 0.4 0.6 10126925 110075.28 331.7759
## 34 0.4 0.7 10279690 111735.76 334.2690
## 35 0.4 0.8 10286904 111814.17 334.3863
## 36 0.4 0.9 10236856 111270.17 333.5718
## 37 0.5 0.1 5973417 64928.45 254.8106
## 38 0.5 0.2 6557447 71276.59 266.9768
## 39 0.5 0.3 7156285 77785.71 278.9009
## 40 0.5 0.4 7631852 82954.91 288.0189
## 41 0.5 0.5 7928171 86175.77 293.5571
## 42 0.5 0.6 8073858 87759.32 296.2420
## 43 0.5 0.7 8135275 88426.90 297.3666
## 44 0.5 0.8 8171712 88822.96 298.0318
## 45 0.5 0.9 8219837 89346.05 298.9081
## 46 0.6 0.1 5304187 57654.21 240.1129
## 47 0.6 0.2 5762651 62637.51 250.2749
## 48 0.6 0.3 6169939 67064.55 258.9682
## 49 0.6 0.4 6453286 70144.42 264.8479
## 50 0.6 0.5 6616639 71919.99 268.1790
## 51 0.6 0.6 6703743 72866.77 269.9385
## 52 0.6 0.7 6756874 73444.28 271.0061
## 53 0.6 0.8 6798782 73899.81 271.8452
## 54 0.6 0.9 6832463 74265.90 272.5177
## 55 0.7 0.1 4728148 51392.91 226.7000
## 56 0.7 0.2 5079103 55207.64 234.9631
## 57 0.7 0.3 5360119 58262.16 241.3756
## 58 0.7 0.4 5539526 60212.24 245.3818
## 59 0.7 0.5 5637022 61271.97 247.5318
## 60 0.7 0.6 5684147 61784.20 248.5643
## 61 0.7 0.7 5700830 61965.55 248.9288
## 62 0.7 0.8 5692439 61874.34 248.7455
## 63 0.7 0.9 5656418 61482.81 247.9573
## 64 0.8 0.1 4238608 46071.82 214.6435
## 65 0.8 0.2 4506361 48982.19 221.3192
## 66 0.8 0.3 4702920 51118.69 226.0944
## 67 0.8 0.4 4816920 52357.83 228.8183
## 68 0.8 0.5 4868282 52916.10 230.0350
## 69 0.8 0.6 4877996 53021.69 230.2644
## 70 0.8 0.7 4857726 52801.37 229.7855
## 71 0.8 0.8 4812666 52311.59 228.7173
## 72 0.8 0.9 4747042 51598.28 227.1526
## 73 0.9 0.1 3827248 41600.52 203.9621
## 74 0.9 0.2 4033070 43837.71 209.3746
## 75 0.9 0.3 4172580 45354.13 212.9651
## 76 0.9 0.4 4244502 46135.90 214.7927
## 77 0.9 0.5 4267077 46381.27 215.3631
## 78 0.9 0.6 4257145 46273.32 215.1123
## 79 0.9 0.7 4225968 45934.43 214.3232
## 80 0.9 0.8 4182088 45457.47 213.2076
## 81 0.9 0.9 4133918 44933.90 211.9762
##
## $n.optimum
## [1] "Alpha optimum adalah 0.9 ,Beta optimum adalah 0.1"df_des <- HoltWinters(train.ts, alpha = 0.9, beta = 0.1, gamma = F)
df_des## Holt-Winters exponential smoothing with trend and without seasonal component.
##
## Call:
## HoltWinters(x = train.ts, alpha = 0.9, beta = 0.1, gamma = F)
##
## Smoothing parameters:
## alpha: 0.9
## beta : 0.1
## gamma: FALSE
##
## Coefficients:
## [,1]
## a 9960.6782
## b 13.7826datades <- data.frame(train.ts, c(NA,NA, df_des$fitted[,1]))
colnames(datades) <- c("y","yhat")
head(datades)## y yhat
## 1 7797.63 NA
## 2 7773.26 NA
## 3 7576.27 7748.890
## 4 7420.72 7553.626
## 5 7545.40 7382.143
## 6 7548.22 7491.900fore.des <- forecast(df_des, h = 23)# Akurasi DES
selisih.des <- fore.des$mean - test
sse.des <- sum(selisih.des^2)
mse.des <- sse.des / length(test)
rmse.des <- sqrt(mse.des)# Plot
plot(train.ts, xlab="periode waktu", ylab="Yt", col="blue", lty=3)
points(train.ts)
lines(datades[,2], col = "red",lwd = 1.5, lty = 1)
title("Double Exponential Smoothing (Alpha = 0.9 dan Beta = 0.1)", cex.main=1, font.main=4
,col.main="black")
legend("topleft", c("Data aktual","Fitted DES"), lty=1:3,col=c ("blue","red"))
Holt Winters Method
multi <- HoltWinters(train.ts, seasonal = "multiplicative")
fore.multi <- forecast(multi, h = 23)# Akurasi
selisih<-as.numeric(fore.multi$mean)-as.numeric(test.ts)
selisih## [1] 71.71791 180.80142 380.64213 411.52909 692.72212 740.53418
## [7] 746.70011 477.30260 401.41512 662.63677 983.77328 941.05463
## [13] 962.19453 1019.68502 1034.23521 1008.27105 1141.07786 1178.08853
## [19] 1050.04472 916.26935 980.09073 1172.73768 1233.87886sse.holt <- sum(selisih^2)
mse.holt <- sse.holt / length(selisih)
rmse.holt <- sqrt(mse.holt)Perbandingan Akurasi Smoothing
tabel1 <- matrix(c(sse.dma, mse.dma, rmse.dma, sse.des, mse.des, rmse.des,
sse.holt, mse.holt, rmse.holt), ncol = 3, byrow = T)
colnames(tabel1) <- c("SSE","MSE","RMSE")
rownames(tabel1) <- c("DMA","DES","Holt")
kable(tabel1)| SSE | MSE | RMSE | |
|---|---|---|---|
| DMA | 4346898 | 48841.55 | 221.0012 |
| DES | 6286996 | 273347.66 | 522.8266 |
| Holt | 17161598 | 746156.45 | 863.8035 |
Kestasioneran 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)
Dilihat dari grafik time-series, data harga beras 2013-2022 tidak stasioner dalam nilai tengah maupun ragam
ACF & PACF
acf(data.ts[,1], lag.max = 35)
pacf(data.ts[,1], lag.max = 35)
Plot ACF menunjukkan pola Tailing off slowly, yang mengindikasikan bahwa data tidak stasioner, selanjutnya akan dilakukan uji statistik untuk memastikan.
Uji Statistik
adf.test(data.ts[,1])##
## Augmented Dickey-Fuller Test
##
## data: data.ts[, 1]
## Dickey-Fuller = -2.1214, Lag order = 4, p-value = 0.5262
## alternative hypothesis: stationaryUji Statistik menghasilkan p-value yang lebih besar dari 0.05 yang berarti pada taraf nyata 5% dapat dikatakan data tidak stasioner sehingga diperlukan penanganan ketidakstasioneran data sebelum meentukan model.
Penanganan Stasioneritas
Differencing
data.diff <- diff(data.ts[,1], differences = 1)Cek Stasioneritas Data setelah Differencing
# cek Stasioner Eksploratif
plot(data.diff, type = "l")
# Uji Statistik
adf.test(data.diff)## Warning in adf.test(data.diff): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: data.diff
## Dickey-Fuller = -5.7563, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationaryDari grafik time-series terlihat setelah dilakukan differencing, data sudah cenderung stasioner dalam nilai tengah, dan hasil test ADF memberikan p-value < 0.05 sehingga dapat dipastikan bahwa data sudah stasioner.
ACF & PACF plot
acf(data.diff, lag.max = 20)
pacf(data.diff, lag.max = 20)
Secara eksploratif, dianggap plot ACF diatas memiliki pola “cut-off” pada lag ke-2 dan plot PACF nya berpola “tail-off”. Sehingga didapat model tentatif ARIMA(0,1,2)
EACF plot
eacf(data.diff)## 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 x o x 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 x x o o o o o o o o o o o
## 7 o o x o o o o x o o o o o oDari plot EACF kita bisa mendapat model tentatif campuran dengan melihat notasi “o” pada pojok kiri atas yang membentuk pola segitiga. Dari plot diatas kita mengambil ARIMA(0,1,4) dan ARIMA(2,1,1) sebagai model tentatif.
Model Tentatif
ARIMA(0,1,2)
model1 <- Arima(data.diff, order = c(0,0,2))
summary(model1)## Series: data.diff
## ARIMA(0,0,2) with non-zero mean
##
## Coefficients:
## ma1 ma2 mean
## 0.5122 0.1085 16.8327
## s.e. 0.0950 0.1187 23.5814
##
## sigma^2 = 24981: log likelihood = -737.55
## AIC=1483.09 AICc=1483.46 BIC=1494.04
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.04198917 155.9613 118.0654 -139.877 328.9805 0.8443722
## ACF1
## Training set -0.005396536Model diatas memiliki persamaan : \[ \nabla Y=e_t-0.51e_{t-1}+0.11e_{t-2} \]
ARIMA(0,1,4)
model2 <- Arima(data.diff, order = c(0,0,4))
summary(model2)## Series: data.diff
## ARIMA(0,0,4) with non-zero mean
##
## Coefficients:
## ma1 ma2 ma3 ma4 mean
## 0.4244 -0.0019 -0.3581 -0.3532 17.5829
## s.e. 0.0868 0.0901 0.0872 0.0829 9.9341
##
## sigma^2 = 22077: log likelihood = -729.75
## AIC=1471.51 AICc=1472.29 BIC=1487.92
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.4629717 145.2892 108.0816 -76.17541 312.1493 0.7729707
## ACF1
## Training set 0.004736279Didapat model dengan persamaan : \[ \nabla Y=e_t-0.42e_{t-1}+0.0019e_{t-2}+0.39e_{t-3}+0.35e_{t-4} \] ### ARIMA(2,1,1)
model3 <- Arima(data.diff, order = c(2,0,1))
summary(model3)## Series: data.diff
## ARIMA(2,0,1) with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 mean
## 1.1005 -0.5284 -0.6904 18.0037
## s.e. 0.1216 0.0789 0.1272 10.1927
##
## sigma^2 = 22552: log likelihood = -731.38
## AIC=1472.75 AICc=1473.31 BIC=1486.43
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.9697895 147.5164 111.3264 -119.9052 369.9621 0.7961764
## ACF1
## Training set 0.007336485Didapat model dengan persamaan : \[ \nabla Y=1.1Y_{t-1}-0.53Y_{t-2}+e_t+0.69e_{t-1} \]
Model dengan parameter optimum
model4 <- auto.arima(data.diff)
summary(model4)## Series: data.diff
## ARIMA(2,0,1)(1,0,0)[12] with zero mean
##
## Coefficients:
## ar1 ar2 ma1 sar1
## 1.0271 -0.4514 -0.6601 0.2615
## s.e. 0.1454 0.0869 0.1473 0.0952
##
## sigma^2 = 21560: log likelihood = -729.13
## AIC=1468.26 AICc=1468.81 BIC=1481.94
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 15.58471 144.2338 108.9301 -35.74622 311.8069 0.7790387
## ACF1
## Training set -0.02334064Dengan menggunakan fingsi auto.arima() kita bisa mendapatkan secara langsung model dengan parameter optimum dan didapat model dengan persamaan : \[
\nabla Y=1.02Y_{t-1}-0.45Y_{t-2}+0.66e_{t-1}+e_t
\] ## Uji Signifikansi Model
coeftest(model1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.51215 0.09504 5.3888 7.093e-08 ***
## ma2 0.10846 0.11868 0.9139 0.3608
## intercept 16.83272 23.58143 0.7138 0.4753
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(model2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 0.4243901 0.0867964 4.8895 1.011e-06 ***
## ma2 -0.0019207 0.0900930 -0.0213 0.98299
## ma3 -0.3581145 0.0871844 -4.1076 3.999e-05 ***
## ma4 -0.3532249 0.0828568 -4.2631 2.016e-05 ***
## intercept 17.5828644 9.9341410 1.7699 0.07674 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(model3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.100492 0.121642 9.0469 < 2.2e-16 ***
## ar2 -0.528377 0.078917 -6.6954 2.151e-11 ***
## ma1 -0.690367 0.127188 -5.4279 5.701e-08 ***
## intercept 18.003720 10.192698 1.7663 0.07734 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1coeftest(model4)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.027098 0.145387 7.0646 1.611e-12 ***
## ar2 -0.451426 0.086864 -5.1969 2.026e-07 ***
## ma1 -0.660050 0.147299 -4.4810 7.428e-06 ***
## sar1 0.261518 0.095187 2.7474 0.006007 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Perbandingan Akurasi
AIC dan BIC
aic1 <- model1$aic
aic2 <- model2$aic
aic3 <- model3$aic
aic4 <- model4$aic
bic1 <- model1$bic
bic2 <- model2$bic
bic3 <- model3$bic
bic4 <- model4$bic
tabel <- matrix(c(aic1,aic2,aic3,aic4,bic1,bic2,bic3,bic4), ncol = 4, byrow = T)
colnames(tabel) <- c("Model 1","Model 2","Model 3","Model 4")
rownames(tabel) <- c("AIC","BIC")
kable(tabel)| Model 1 | Model 2 | Model 3 | Model 4 | |
|---|---|---|---|---|
| AIC | 1483.092 | 1471.506 | 1472.751 | 1468.257 |
| BIC | 1494.037 | 1487.924 | 1486.432 | 1481.938 |
MAPE
mape1 <- 328.9805
mape2 <- 312.1493
mape3 <- 369.9621
mape4 <- 311.8069
tabel2 <- matrix(c(mape1, mape2, mape3, mape4), ncol = 4)
rownames(tabel2) <- "MAPE"
colnames(tabel2) <- c("ARIMA(0,1,2)","ARIMA(0,1,4)","ARIMA(2,1,1)","SARIMA(2,1,0)(1,1,0)")
kable(tabel2)| ARIMA(0,1,2) | ARIMA(0,1,4) | ARIMA(2,1,1) | SARIMA(2,1,0)(1,1,0) | |
|---|---|---|---|---|
| MAPE | 328.9805 | 312.1493 | 369.9621 | 311.8069 |
Overfitting
over1 <- Arima(train, order = c(3,1,1))
summary(over1)## Series: train
## ARIMA(3,1,1)
##
## Coefficients:
## ar1 ar2 ar3 ma1
## 1.0749 -0.5314 0.0123 -0.6218
## s.e. 0.2730 0.2071 0.1542 0.2522
##
## sigma^2 = 26898: log likelihood = -591.44
## AIC=1192.88 AICc=1193.59 BIC=1205.43
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 28.23258 159.4859 120.3887 0.3033602 1.320557 0.9171439 -0.0282721coeftest(over1)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.074909 0.272987 3.9376 8.23e-05 ***
## ar2 -0.531378 0.207093 -2.5659 0.01029 *
## ar3 0.012254 0.154160 0.0795 0.93664
## ma1 -0.621848 0.252205 -2.4656 0.01368 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1over2 <- Arima(train, order = c(2,1,2))
summary(over2)## Series: train
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## 1.0523 -0.5081 -0.5998 -0.0119
## s.e. 0.1722 0.1498 0.1992 0.1675
##
## sigma^2 = 26898: log likelihood = -591.44
## AIC=1192.88 AICc=1193.59 BIC=1205.43
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 28.20955 159.4869 120.4067 0.3031096 1.320743 0.9172805
## ACF1
## Training set -0.02761302coeftest(over2)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 1.052312 0.172203 6.1109 9.909e-10 ***
## ar2 -0.508062 0.149811 -3.3913 0.0006955 ***
## ma1 -0.599785 0.199159 -3.0116 0.0025989 **
## ma2 -0.011916 0.167498 -0.0711 0.9432851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Perbandingan AIC
tabel3 <- data.frame("Model" = c("ARIMA(2,1,1)","ARIMA(3,1,1)","ARIMA(2,1,2)"),
"AIC" = c(model3$aic, over1$aic, over2$aic))
kable(tabel3)| Model | AIC |
|---|---|
| ARIMA(2,1,1) | 1472.751 |
| ARIMA(3,1,1) | 1192.879 |
| ARIMA(2,1,2) | 1192.881 |
Analisis Sisaan
residual <- over1$residuals
qqnorm(residual)
qqline(residual, col="red", lwd=2, lty=1)
plot(residual, type="o",ylab="Sisaan", xlab="Order",main="Plot Sisaan vs Order")
abline(h=0,col="red")
acf(residual, main="Plot ACF Residual")
pacf(residual, main="Plot PACF Residual")
Forecasting
ramalan <- forecast(over1, h = 23)
ramalan## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## 93 9975.676 9765.495 10185.86 9654.232 10297.12
## 94 9973.393 9602.653 10344.13 9406.395 10540.39
## 95 9964.397 9489.938 10438.86 9238.775 10690.02
## 96 9956.099 9425.281 10486.92 9144.283 10767.91
## 97 9951.932 9390.861 10513.00 9093.848 10810.02
## 98 9951.753 9370.459 10533.05 9062.740 10840.76
## 99 9953.672 9353.537 10553.81 9035.845 10871.50
## 100 9955.779 9334.342 10577.22 9005.373 10906.19
## 101 9957.022 9311.289 10602.76 8969.459 10944.59
## 102 9957.262 9285.772 10628.75 8930.307 10984.22
## 103 9956.886 9260.054 10653.72 8891.173 11022.60
## 104 9956.368 9235.675 10677.06 8854.162 11058.57
## 105 9956.015 9213.067 10698.96 8819.774 11092.26
## 106 9955.906 9191.924 10719.89 8787.497 11124.32
## 107 9955.970 9171.707 10740.23 8756.543 11155.40
## 108 9956.093 9151.980 10760.21 8726.308 11185.88
## 109 9956.189 9132.531 10779.85 8696.513 11215.86
## 110 9956.228 9113.339 10799.12 8667.141 11245.32
## 111 9956.220 9094.469 10817.97 8638.286 11274.15
## 112 9956.193 9075.990 10836.40 8610.039 11302.35
## 113 9956.167 9057.928 10854.41 8582.429 11329.91
## 114 9956.155 9040.268 10872.04 8555.426 11356.88
## 115 9956.154 9022.970 10889.34 8528.972 11383.34data.ramalan <- ramalan$mean
plot(ramalan)
accuracy(data.ramalan, test)## ME RMSE MAE MPE MAPE
## Test set -317.2596 346.3558 317.2596 -3.312673 3.312673