Penggunaan Metode Arima Box-Jenkins Untuk Meramal Harga Minyak Goreng Kemasan di Indonesia Tahun 2022

Anggota Kelompok 4 - Paralel 2:

1. Dicky Lihardo Girsang

2. Syifa Maulany

3. Alifa Shakila

4. Nabil Izzany

5. M. Ficky Haris Ardiansyah

Pendahuluan

Latar Belakang

Minyak goreng adalah bahan pangan dengan komposisi utama dari trigliserida dengan atau tanpa perubahan kimiawi. Pada umumnya berbentuk cair pada suhu ruang dan digunakan untuk menggoreng makanan (Chairunisa 2013). Minyak goreng mengandung zat penting untuk menjaga kesehatan tubuh manusia serta berperan memberi nilai kalori paling besar di antara zat gizi lainnya. Minyak goreng umumnya diperoleh dari tumbuhan seperti kelapa, kelapa sawit, kacang-kacangan, jagung dan kanola.

Minyak goreng di Indonesia dianggap sebagai salah satu bahan pokok rakyat Indonesia. Awal tahun 2008, harga minyak goreng di Indonesia mengalami kenaikan secara signifikan. Salah satu penyebabnya adalah kenaikan harga Crude Palm Oil (CPO) di pasar internasional. Kemudian kembali terjadi lonjakan harga pada akhir tahun 2021 hingga awal tahun 2022 karena terbatasnya produksi kelapa sawit. (USDA 2021) Direktur KPBN, seperti yang dikutip di beberapa situs media, dan laporan Outlook CPOPC 2022 mengaitkan turunnya produktivitas sawit dengan tiga faktor utama, yakni kurangnya tenaga kerja, cuaca buruk, serta meningkatnya harga pupuk 50-80% di pertengahan tahun 2021.

Tujuan Penelitian

Tujuan dari penelitian ini adalah sebagai berikut: 1. Membangun model ARIMA yang dapat menjelaskan hubungan antara harga minyak goreng dunia saat ini dengan harga minyak goreng dunia pada periode sebelumnya 2. Melakukan peramalan harga minyak goreng dunia untuk periode 2021 - 2022 berdasarkan model terbaik

Pre-processing Data

Data yang digunakan adalah data harga minyak goreng kemasan di Indonesia. Data ini dipilih untuk dianalisis lebih lanjut mengenai forecasting dengan model ARIMA.

Library

library(tseries)
library(forecast)
library(TTR)
library(TSA)
library(imputeTS)
library(ggplot2)
library(lmtest)
library(readxl)

Input data

setwd("C:/Users/Dicky Girsang/Desktop/Semester 5/STA1341 Metode Peramalan Deret Waktu")
data <- read_excel("Data Migor Rapi.xlsx")

data.ts <- ts(data$Yt)
head(data.ts)

## Time Series:
## Start = 1 
## End = 6 
## Frequency = 1 
## [1] 12664 12607 12554 12531 12441 12461

Splitting Data ( Training dan Testing )

training <- data[1:106,]
testing <- data[107:117,]

training.ts <- ts(training[,2])
testing.ts <- ts(testing[,2])

Plot Data

#Plot time series
ggplot(data= data, aes(x=t, y=Yt, group=1)) +
  geom_line()+
  geom_point()+
  labs(title= "Time Series Plot Harga Minyak Goreng Kemasan di Indonesia Tahun 2013 - 2022 (September)",
       x= "Periode", y = "Harga (Rp)")

par(mfrow=c(1,2))

ggplot(data= training, aes(x=t, y=Yt, group=1)) +
  geom_line()+
  geom_point()+
  labs(title= "Plot data training (1-106)",
       x= "Periode", y = "Harga (Rp)")

ggplot(data= testing, aes(x=t, y=Yt, group=1)) +
  geom_line()+
  geom_point()+
  labs(title= "Plot data testing (107-117)",
       x= "Periode", y = "Harga (Rp)")

Berdasarkan plot di atas, secara observatif bahwa data tidak stasioner dalam rata-rata dan ragam. Terdapat inkonsistensi lebar pita hingga tahun 2022 dan indikasi pola trend naik.

Pemulusan

Diperlukan pencarian nilai n optimum untuk melakukan pemulusan. Metode Single Moving Average (SMA) dan Double Moving Average (DMA) akan diterapkan untuk pemulusan data harga minyak goreng. Proses pemulusan akan dilakukan terhadap data training, selanjutnya dilakukan perhitungan akurasi forecasting pada data testing. Pada fungsi di bawah ini didapatkan nilai optimisasi nilai n = 2.

Single Moving Average (SMA)

data.sma <- SMA(data.ts, n=2)
data.ramal <- c(NA, data.sma)

data.gab<-cbind(n = 1:122,aktual=c(data.ts,rep(NA,5)),pemulusan=c(data.sma,rep(NA,5)),ramalan=c(data.ramal,rep(data.ramal[length(data.ramal)],4)))

ts.plot(data.gab[,2], xlab="Periode Waktu", ylab="Harga", main= "SMA Migor")
points(data.gab[,2])
lines(data.gab[,3],col="green",lwd=2)
lines(data.gab[,4],col="red",lwd=2)
legend("topleft",c("data aktual","data pemulusan","data peramalan"), lty=8, col=c("black","green","red"), cex=0.8)

Plot di atas menunjukkan hasil pemulusan dan peramalan metode SMA pada n=2. Terlihat bahwa hasil garis putus hijau yang menunjukkan hasil pemulusan cukup menggambarkan pola data harga minyak goreng, namun pada garis putus merah yang menunjukkan hasil peramalan kurang menggambarkan pola data aktual (data testing).

# Akurasi SMA
data.ts <- as.numeric(data.ts)
error.sma = data.ts-data.ramal[1:length(data.ts)]
SSE.sma = sum(error.sma[3:length(data.ts)]^2)
MSE.sma = mean(error.sma[3:length(data.ts)]^2)
MAPE.sma = mean(abs((error.sma[3:length(data.ts)]/data.ts[3:length(data.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 = 2")
akurasi.sma

##      Akurasi m = 2
## SSE   1.009510e+08
## MSE   8.778346e+05
## MAPE  1.811912e+00

Single Exponential Smoothing

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.9802239
##  beta : FALSE
##  gamma: FALSE
## 
## Coefficients:
##       [,1]
## a 15998.24

plot(sesopt)

ramalanopt<- forecast(sesopt, h=10)
ramalanopt

##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 107       15998.24 15692.24 16304.23 15530.26 16466.22
## 108       15998.24 15569.75 16426.72 15342.92 16653.55
## 109       15998.24 15475.20 16521.27 15198.32 16798.15
## 110       15998.24 15395.30 16601.17 15076.12 16920.35
## 111       15998.24 15324.81 16671.66 14968.32 17028.15
## 112       15998.24 15261.04 16735.44 14870.79 17125.69
## 113       15998.24 15202.35 16794.12 14781.04 17215.43
## 114       15998.24 15147.71 16848.76 14697.47 17299.01
## 115       15998.24 15096.37 16900.10 14618.95 17377.52
## 116       15998.24 15047.80 16948.67 14544.67 17451.80

Selain itu, dilakukan juga penghitungan parameter alpha secara manual. Setelah dilakukan looping, diperoleh hasil sebagai berikut:

alpha	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	Optimum (0,98)
SES	3.485657e+07	2.019821e+07	1.369971e+07	10397124.000	8532202.0000	7408135.0000	6714562.0000	6300378.0000	6088830.0000	6039313.0000
MSE	3.288355e+05	1.905492e+05	1.292426e+05	98086.070	80492.4700	69888.0700	63344.9300	59437.5300	57441.7900	56974.6500
RMSE	5.734418e+02	4.365194e+02	3.595032e+02	313.187	283.7119	264.3635	251.6842	243.7981	239.6702	238.6936

Parameter optimum pada data training yang didapat dari metode SES berdasarkan nilai RMSE terkecil adalah α = 0.98

Double Exponential Smoothing

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: 0.9664833
##  beta : 0.04055307
##  gamma: FALSE
## 
## Coefficients:
##          [,1]
## a 15999.24005
## b    63.76021

plot(desopt)

ramalandesopt<- forecast(desopt, h=10)
ramalandesopt

##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 107       16063.00 15751.65 16374.35 15586.84 16539.16
## 108       16126.76 15685.20 16568.32 15451.45 16802.07
## 109       16190.52 15642.06 16738.98 15351.72 17029.32
## 110       16254.28 15610.29 16898.28 15269.38 17239.19
## 111       16318.04 15585.21 17050.87 15197.28 17438.80
## 112       16381.80 15564.47 17199.13 15131.80 17631.80
## 113       16445.56 15546.67 17344.45 15070.83 17820.29
## 114       16509.32 15530.93 17487.71 15013.00 18005.64
## 115       16573.08 15516.63 17629.53 14957.38 18188.78
## 116       16636.84 15503.36 17770.33 14903.33 18370.36

Parameter optimum pada data training yang didapat dari metode DES berdasarkan nilai RMSE terkecil adalah α = 0.5957 dan β = 0.4887.

Perbandingan akurasi pemulusan SES dan DES

Akurasi <- data.frame("DES"=akurasi.des1[,1], "SES"=akurasi.ses1[,1])
rownames(Akurasi) <- c("SSE","MSE", "RMSE")
Akurasi

##               DES          SES
## SSE  6170601.9457 6039313.1539
## MSE    58213.2259   56974.6524
## RMSE     241.2742     238.6936

Berdasarkan perbandingan nilai RMSE kedua metode dengan parameter optimum, metode SES dinyatakan lebih baik dibandingkan metode SES karena memiliki nilai RMSE yang lebih kecil.

Kestasioneran Data

Plot ACF dan PACF

Plot Autocorrelation Function

Plot ACF dapat digunakan secara observatif untuk melihat apakah data sudah stasioner dalam mean dan varians atau sebaliknya. Apabila plot berada dalam nilai mean, maka dikatakan bahwa data telah stasioner dalam mean. Konsep yang serupa juga berlaku pada kestasioneran dalam ragam

#Memeriksa data
acf(data.ts, lag.max = 48, plot = F) %>% 
  autoplot() + 
  ggtitle("Plot ACF Harga Minyak Goreng Kemasan di Indonesia tahun 2013 - 2022\n(sebelum dilakukan diferensiasi)") + 
  theme_classic()

Dari Plot ACF di atas, pola menurun pada grafik menunjukkan bahwa terdapat indikasi ketidakstasioneran dalam rataan pada data yang digunakan.

Plot Partial Autocorrelation Function

pacf(data.ts)

Pada plot PACF terlihat bahwa pola cut-off di lag 1, 7, 8, 9, dan 15. Oleh karena itu, diperlukan proses Diferensiasi untuk mentransformasi data yang dianalisis menjadi bentuk stasioner. Sebelum itu, sebagai bentuk validasi, dilakukan uji Stasioneritas Augmented Dickey-Fuller

Uji ADF

Salah satu uji akar unit yang digunakan untuk melihat kestasioneran data time-series dalam mean (rataan) adalah uji ADF (Box et al.). Uji ini digunakan karena telah mempertimbangkan kemungkinan adanya autokorelasi pada error term jika pola data yang digunakan non-stasioner Hipotesis yang digunakan dalam uji formal ADF H0 : data tidak stasioner H1 : data stasioner

adf.test(data.ts, k=12)   #lag order default adalah k=4

## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.ts
## Dickey-Fuller = -0.70882, Lag order = 12, p-value = 0.9673
## alternative hypothesis: stationary

Keputusan: Diperoleh nilai p-value dari hasil uji di atas adalah 0.9673 > 0.05 (pada tingkat kepercayaan 95%), sehingga disimpulkan bahwa tak tolak H0. Artinya, terdapat ketidaksioneran pada data tersebut atau masih terdapat unit root pada data dan perlu dilakukan diferensiasi terlebih dahulu.

Diferensiasi

Perilaku data yang stasioner antara lain tidak mempunyai variasi yang terlalu besar dan mempunyai kecenderungan untuk mendekati nilai rata-ratanya, dan sebaliknya untuk data yang tidak stasioner. Diferensiasi atau pembedaan merupakan jawaban untuk mengatasi ketidakstasioneran pada data.

Diferensiasi 1

#Differensiasi
data.dif1 <- diff(training.ts, differences = 1)
data.dif1 <- ts(data.dif1, frequency = 12)

Identifikasi kestasioneran data setelah dilakukan differencing 1

plot.ts(data.dif1,lty=1,xlab="Periode", ylab="Harga Minyak Goreng Kemasan\nDifferencing 1")

acf(data.dif1, lag.max = 25, plot = F) %>% 
  autoplot() + 
  ggtitle("Plot ACF Diff 1") + 
  theme_classic()

adf.test(data.dif1)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.dif1
## Dickey-Fuller = -3.8467, Lag order = 4, p-value = 0.01929
## alternative hypothesis: stationary

Berdasarkan uji ADF di atas, diperoleh nilai p-value (0.01929) < 0.05, artinya data sudah stasioner dan tidak perlu kembali dilakukan diferensi kedua kalinya.

Plot EACF

eacf(data.dif1)

## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 o o o o o o o o o o o  x  o  o 
## 1 x o o o o o o o o o o  x  o  o 
## 2 x x o o o o o o o o o  x  o  o 
## 3 x o x o o o o o o o o  x  o  o 
## 4 x o x o o o o o o o o  x  o  o 
## 5 o x x o o o o o o o o  x  o  o 
## 6 x o x o o o o o o o o  x  o  o 
## 7 x x o o o o o o o o o  x  o  o

Secara teoritis model ARMA(p,q) mempunyai pola segitiga nol dimana nilai pada pojok kiri atas bersesuaian dengan ordo ARMA.

Pemodelan ARIMA

model1 <- Arima(training.ts, order = c(0,1,0), method = "ML")
model2 <- Arima(training.ts, order = c(3,1,1), method = "ML")
model3 <- Arima(training.ts, order = c(5,1,0), method = "ML")
model4 <- Arima(training.ts, order = c(2,1,0), method = "ML")
model5 <- auto.arima(training.ts)
summary(model5)

## Series: training.ts 
## ARIMA(0,1,0) 
## 
## sigma^2 = 57547:  log likelihood = -724.41
## AIC=1450.81   AICc=1450.85   BIC=1453.47
## 
## Training set error measures:
##                    ME     RMSE      MAE       MPE      MAPE      MASE
## Training set 31.59117 238.7548 95.59117 0.2084013 0.6763793 0.9918056
##                     ACF1
## Training set -0.04253585

Dari model auto-arima di atas, diperoleh nilai parameter p,i,q yang sama dengan model1.

lmtest::coeftest(model2)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)    
## ar1 -0.932617   0.097574 -9.5580  < 2e-16 ***
## ar2  0.099353   0.130649  0.7605  0.44698    
## ar3  0.194601   0.096274  2.0213  0.04325 *  
## ma1  0.986600   0.053130 18.5695  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lmtest::coeftest(model3)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 -0.019162   0.096443 -0.1987   0.8425
## ar2  0.115042   0.096114  1.1969   0.2313
## ar3  0.075616   0.095926  0.7883   0.4305
## ar4 -0.056984   0.095180 -0.5987   0.5494
## ar5  0.112875   0.095565  1.1811   0.2375

lmtest::coeftest(model4)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 -0.021558   0.096411 -0.2236   0.8231
## ar2  0.116869   0.096128  1.2158   0.2241

Pemilihan model terbaik

banding <- matrix(c(model1$aic, model2$aic,
                    model3$aic, model4$aic
                   ), nrow = 4, byrow = T)
colnames(banding) <- c("AIC")
row.names(banding) <- c("(0,1,0)", "(3,1,1)", "(5,1,0)", "(2,1,0)")
banding

##              AIC
## (0,1,0) 1450.811
## (3,1,1) 1449.265
## (5,1,0) 1456.586
## (2,1,0) 1453.282

Overfitting

Overfitting adalah suatu keadaan data yang digunakan pada data training memiliki prediksi yang terlalu baik, tetapi prediksinya buruk pada data testing. Ketika sebuah model Overfit, model tidak dapat melakukan generalisasi dengan baik sehingga apabila dilakukan tes dengan menggunakan data yang berbeda dapat mengurangi akurasi (hasil yang dibuat tidak sesuai yang diharapkan).

Model overfitting adalah dengan menambahkan nilai P dan Q pada model terbaik sementara, yakni model2 ARIMA(1,1,3). Adapun pendefinisian model tersebut adalah sebagai berikut:

over1 <- Arima(training.ts, order = c(4,1,1), method = "ML")
over2 <- Arima(training.ts, order = c(3,1,2), method = "ML")

banding.over <- matrix(c(model2$aic, over1$aic, over2$aic), nrow = 3, byrow = T)
colnames(banding.over) <- "AIC"
row.names(banding.over) <- c("(3,1,1)","(4,1,1)", "(3,1,2)")
banding.over

##              AIC
## (3,1,1) 1449.265
## (4,1,1) 1451.264
## (3,1,2) 1451.264

Berdasarkan nilai AIC pada matriks di atas, diketahui bahwa model ARIMA(3,1,1) adalah model dengan AIC terkecil. Oleh karena itu, disimpulkan bahwa ARIMA(3,1,1) adalah model terbaik dari data harga minyak goreng.

Analisis Sisaan Model Terbaik

residual <- model2$residuals

plot(residual, type="o",ylab="Sisaan", xlab="Order",main="Plot Sisaan vs Order")
abline(h=0,col="red")

Box.test(residual)

## 
##  Box-Pierce test
## 
## data:  residual
## X-squared = 0.019575, df = 1, p-value = 0.8887

qqnorm(residual)
qqline(residual, col="red", lwd=2, lty=1)

acf(residual, main = "Plot ACF Residual")

pacf(residual, main = "Plot PACF Residual")

Forecasting

ramalan <- forecast(model2, h = 11)
ramalan

##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 107       15994.74 15697.30 16292.18 15539.85 16449.63
## 108       15990.98 15559.12 16422.83 15330.51 16651.44
## 109       16011.67 15469.35 16553.99 15182.26 16841.08
## 110       15990.97 15332.46 16649.49 14983.86 16998.09
## 111       16011.60 15272.48 16750.72 14881.21 17141.99
## 112       15994.34 15163.97 16824.70 14724.40 17264.27
## 113       16008.46 15109.81 16907.11 14634.09 17382.83
## 114       15997.59 15024.62 16970.55 14509.56 17485.61
## 115       16005.77 14971.61 17039.93 14424.16 17587.38
## 116       15999.81 14902.27 17097.35 14321.26 17678.35
## 117       16004.07 14850.51 17157.62 14239.86 17768.27

Berdasarkan hasil forecasting di atas, harga ramalan minyak goreng di Indonesia, berfluktuasi di harga 16000 rupiah.

plot(ramalan)

Pemulusan dan Pemodelan

coming soon

Daftar Pustaka

Box GEP, Jenkins GM, Reinsel GC, Ljung GM. 2016. Time Series Analysis - Forecasting & Control, 5th ed. New Jersey: John Wiley & Sons, Inc.