library(readxl)
library(forecast)
library(lmtest)
library(padr)
library(tseries)
library(TSA)
library(tidyverse)
library(MASS)
library(scales)Data
Data yang digunakan data harga minyak goreng (y), data harga CPO Internasional (x1) dan data harga CPO domestik (x2). Semua data berupa data bulanan periode Januari 2013 hingga Oktober 2022. Untuk keperluan analisis, data dibagi menjadi dua yaitu data pemodelan dan validasi. Data pemodelan diambil dari periode Januari 2013 hingga Januari 2022. Data validasi diambil dari periode Februari 2022 hingga Oktober 2022.
setwd("C:/Users/Nabil Izzany/Documents/Kuylah/Semester 5/MPDW/Kerkom/Buat dikumpulin")
data <- readxl::read_xlsx("Data Migor Rapi.xlsx")
data <- data[,c(-1,-6)]
head(data)## # A tibble: 6 × 5
## t Yt X1 X2 Keterangan
## <dbl> <dbl> <dbl> <dbl> <chr>
## 1 1 12664 8.53 7.66 Train
## 2 2 12607 8.90 7.92 Train
## 3 3 12554 8.65 7.74 Train
## 4 4 12531 8.50 7.59 Train
## 5 5 12441 8.57 7.73 Train
## 6 6 12461 8.66 7.93 Trainlog.data <- data[,2:4]
log.data <- log(log.data)
log.data <- as.ts(log.data)
log.data.ket <- data.frame(t = data$t,Yt = log.data[,1], X1 = log.data[,2], X2 = log.data[,3], Keterangan = data$Keterangan)Eksplorasi Data
q <- ggplot(data, aes(x=t, y=Yt, color=Keterangan ))+
geom_path(aes(group = 1))+
geom_point()+
xlab("Periode")+
ylab("Harga (Rp)")+
theme(legend.position = "bottom")
q+theme_classic()## Warning: Removed 1 row containing missing values (`geom_path()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).
v <- ggplot(log.data.ket, aes(x=t, y=Yt, color=Keterangan ))+
geom_path(aes(group = 1))+
geom_point()+
xlab("Periode")+
ylab("Harga (Rp)")+
theme(legend.position = "bottom")
v+theme_classic()## Don't know how to automatically pick scale for object of type <ts>. Defaulting
## to continuous.## Warning: Removed 1 row containing missing values (`geom_path()`).
## Removed 1 rows containing missing values (`geom_point()`).
cor(log.data[,1],log.data[,2])## [1] NAcor(log.data[,1], log.data[,3])## [1] NAtrain.log <- as.ts(log.data[1:109,])
test.log <- as.ts(log.data[110:118,])
train <- as.ts(data[1:109, 2:4])
test <- as.ts(data[110:118, 2:4])Transformasi
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## somey <- data$Yt
boxcox(y ~ 1)
ARIMAX
Metode merujuk pada (Umar 2018)
Secara singkat sebagai berikut: 1. Eksplorasi data 2. Pembentukan model regresi 3. Pemodelan sisaan wt menggunakan ARIMA 4. Pembentukan model ARIMAX 5. Evaluasi kebaikan model ARIMAX
Eksplorasi Data
Pembentukan model regresi
reg3 <- lm(train.log[,1]~train[,2]+train[,3], data = NULL)
summary(reg3)##
## Call:
## lm(formula = train.log[, 1] ~ train[, 2] + train[, 3], data = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.12842 -0.02866 0.01068 0.03332 0.20073
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.220391 0.039494 233.461 < 2e-16 ***
## train[, 2] 0.044126 0.011321 3.898 0.000171 ***
## train[, 3] -0.013959 0.008584 -1.626 0.106857
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04759 on 106 degrees of freedom
## Multiple R-squared: 0.6667, Adjusted R-squared: 0.6604
## F-statistic: 106 on 2 and 106 DF, p-value: < 2.2e-16reg4 <- lm(train.log[,1]~0+train[,2]+train[,3], data = NULL)
summary(reg4)##
## Call:
## lm(formula = train.log[, 1] ~ 0 + train[, 2] + train[, 3], data = NULL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5735 -0.3238 0.2183 0.8119 2.2164
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## train[, 2] 2.52525 0.08816 28.64 <2e-16 ***
## train[, 3] -1.76827 0.09373 -18.87 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.075 on 107 degrees of freedom
## Multiple R-squared: 0.9875, Adjusted R-squared: 0.9873
## F-statistic: 4242 on 2 and 107 DF, p-value: < 2.2e-16Pengujian sisaan white nose
H0 : 𝜌1=𝜌2=⋯=𝜌𝐾=0 (sisaan saling bebas atau sisaan white noise)
H1 : minimal ada satu 𝜌𝑘≠0, untuk k = 1, 2, …, K
sisaan2 <- reg4$residuals
plot(sisaan2, type="o",ylab="Sisaan", xlab="Order")
abline(h=0,col="red")
Box.test(sisaan2)##
## Box-Pierce test
##
## data: sisaan2
## X-squared = 54.437, df = 1, p-value = 1.605e-13Nilai p lebih kecil dari taraf nyata 5%, maka dapat disimpulkan bahwa sisaan belum memenuhi asumsi white noise. Selanjutnya dilakukan pemodelan sisaan menggunakan ARIMA.
Pemodelan sisaan model regresi menggunakan ARIMA
acf(sisaan2)
pacf(sisaan2)
adf.test(sisaan2, k=12)##
## Augmented Dickey-Fuller Test
##
## data: sisaan2
## Dickey-Fuller = -2.3857, Lag order = 12, p-value = 0.4166
## alternative hypothesis: stationaryPlot ACF menunjukkan pola yang menurun secara perlahan secara eksponensial yang mengindikasikan bahwa data tidak stasioner, sehingga penentuan ordo untuk model sementara tidak dapat dilakukan.
sisaan.dif2 <- diff(sisaan2, differences = 1)
adf.test(sisaan.dif2, k=12)## Warning in adf.test(sisaan.dif2, k = 12): p-value smaller than printed p-value##
## Augmented Dickey-Fuller Test
##
## data: sisaan.dif2
## Dickey-Fuller = -4.5213, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationaryUji ADF menghasilkan nilai-p sebesar 0.01 yang membuktikan bahwa data tidak mengandung akar unit atau sudah stasioner. Selanjutnya dilakukan identifikasi orde model ARIMAX dengan melihat plot ACF dan PACF serta EACF.
acf(sisaan.dif2, lag.max = 48)
pacf(sisaan.dif2, lag.max = 48)
eacf(sisaan.dif2)## 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 x o o o o o
## 1 x o o o o o o o x x o o o o
## 2 x o o o o o o o x x o o o o
## 3 x o o o o o o o x o o o o o
## 4 x o x x o o o o x o o o o o
## 5 x x x x o o o o x o o o o o
## 6 x x o x o o o o x o o o o o
## 7 x x o x o o o o x o o o o oDilihat dari plot ACF dan PACF belum terlihat jelas ordo AR dan MAnya oleh karena itu, perlu plot eacf untuk mempermudah. Maka dugaan model sementaara berdasarkan plot EACF adalah ARIMAX(0,1,1) dan ARIMAX(1,1,1).
Pembentukan model ARIMAX
#Yg bener
model1 <- Arima(train.log[,1], order = c(0,1,1), xreg = train[,2:3], method = "ML")
model3 <- Arima(train.log[,1], order = c(1,1,1), xreg = train[,2:3], method = "ML")#Y doang di LN #Yg bener
round(coeftest(model1), digits = 3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ma1 -0.019 0.085 -0.219 0.826
## X1 -0.001 0.009 -0.061 0.951
## X2 0.004 0.007 0.521 0.602round(coeftest(model3), digits = 3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.789 0.181 -4.352 <2e-16 ***
## ma1 0.588 0.183 3.215 0.001 ***
## X1 -0.006 0.009 -0.628 0.530
## X2 0.008 0.007 1.070 0.285
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Genap semua di LN
round(model1$aic, 3)## [1] -467.438round(model3$aic, 3)## [1] -469.603Hasil uji signifikansi hanya ARIMAX (1,1,1) karena terdapat dua parameternya lebih kecil dari taraf nyata 5%. Selain itu ARIMAX (1,1,1) memiliki AIC yang paling kecil dibanding model lainnya. Selanjutnya model ARIMAX (1,1,1) akan di-overfitting.
#Yg bener
over3 <- Arima(train.log[,1], order = c(1,1,2), xreg = train[,2:3], method = "ML")
over4 <- Arima(train.log[,1], order = c(2,1,1), xreg = train[,2:3], method = "ML")#Hanya Y di ln
round(coeftest(over3), digits = 3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 -0.012 0.145 -0.082 0.935
## ma1 -0.024 0.097 -0.249 0.804
## ma2 0.638 0.113 5.650 <2e-16 ***
## X1 -0.016 0.006 -2.654 0.008 **
## X2 0.011 0.005 2.008 0.045 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1round(coeftest(over4), digits = 3)##
## z test of coefficients:
##
## Estimate Std. Error z value Pr(>|z|)
## ar1 0.164 0.184 0.890 0.373
## ar2 0.647 0.105 6.166 <2e-16 ***
## ma1 -0.212 0.182 -1.166 0.244
## X1 -0.020 0.007 -2.858 0.004 **
## X2 0.013 0.006 2.337 0.019 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#yg bener
round(over3$aic, 3)## [1] -491.377round(over4$aic, 3)## [1] -492.602Hasil uji signifikansi hanya ARIMAX (2,1,1) lebih banyak parameternya lebih kecil dari taraf nyata 5% serta AIC lebih kecil. Oleh karena itu dipilih ARIMAX (5,1,1).
sisaan.arimax <- over4$residuals
plot(sisaan.arimax, type="o",ylab="Sisaan", xlab="Order")
abline(h=0,col="red")
Box.test(sisaan.arimax)##
## Box-Pierce test
##
## data: sisaan.arimax
## X-squared = 0.034136, df = 1, p-value = 0.8534acf(sisaan.arimax, lag.max = 48)
pacf(sisaan.arimax,lag.max = 48)
Evaluasi kebaikan model ARIMAX
ramalan <- forecast(over4, h = length(test.log[,1]), xreg = test[,2:3])
data.ramalan <- ramalan$mean
data.ramalan <- exp(data.ramalan)
data.ramalan2 <-data.frame("Y" = data.ramalan)
fit <- exp(ramalan$fitted)
fit <- data.frame("Y" = fit)
banding <- matrix(c(test[,1], data.ramalan), ncol = 2, byrow = FALSE)
colnames(banding) <- c("aktual", "ramalan")
knitr::kable(banding)| aktual | ramalan |
|---|---|
| 17685 | 21709.09 |
| 21433 | 24397.79 |
| 26464 | 24863.38 |
| 26101 | 26632.78 |
| 25808 | 28508.27 |
| 24300 | 31848.91 |
| 22500 | 31692.45 |
| 21700 | 34167.03 |
| 21000 | 35280.03 |
Peramalan
peramalan <- Arima(log.data[1:118,1], order = c(2,1,1), xreg = as.matrix(data[1:118, 3:4]), method = "ML")
nov <- forecast(peramalan, h = 1, xreg = as.matrix(data[119, 3:4]))
round(exp(9.980272), digits = 0)## [1] 21596Simpulan
Model ARIMAX lebih baik dibandingkan model ARIMA. Harga CPO domestik dan harga CPO internasional cukup mampu menjelaskan harga minyak goreng pada periode Januari 2013 hingga Januari 2022. Akan tetapi, model cukup buruk jika digunakan untuk menjelaskan harga minyak goreng pada periode Februari 2022 hingga Oktober 2022.
Daftar Pustaka
Umar MA, Afendi FM, Rizki A, Waryanto B. 2018. Penerapan Model ARIMAX pada Data Harga Cabai Merah Keriting di Indonesia [skripsi] . Bogor Institut Pertanian Bogor. http://repository.ipb.ac.id/handle/123456789/93668