Analisis Model ARIMA Data Inflasi Indonesia Bulanan pada Periode Januari 2012 - Oktober 2022

Kelompok 1 Paralel :

  1. Elfina Tan (G5401201005)

  2. M. Daryl Fauzan (G5401201087)

  3. Denanda Aufadlan T (G5401201088)

  4. Annisa Permata S (G5401201091)

  5. Akmal Taufik (G54190002)

  6. Kreshna Bayu A. (G64190098)

2022-11-19

Pendahuluan

Indonesia termasuk ke dalam negara berkembang dengan jumlah penduduk yang sangat banyak. Hal tersebut membuat Indonesia menjadi negara konsumtif. Dalam pertumbuhan ekonomi di Indonesia, pemerintah harus mengaturnya agar tidak terjadi lonjakan harga yang tinggi (adanya inflasi). Pengendalian inflasi sangat penting untuk dilakukan. Jika kita bisa memperkirakan tingkat inflasi kedepannya, maka tindakan preventif bisa dilakukan untuk menanggulangi inflasi. Data inflasi merupakan data deret waktu, maka salah satu metode yang dapat digunakan dalam pemodelan data inflasi untuk melakukan peramalan adalah Autoregressive Integrated Moving Average (ARIMA). ARIMA adalah singkatan dari Autoregressive Integrated Moving Average yang merupakan bagian dari model peramalan data deret waktu yang umum digunakan untuk memodelkan data deret waktu yang stasioner. ARIMA sangat baik ketepatannya untuk peramalan jangka pendek, sedangkan untuk peramalan jangka panjang ketepatan peramalannya kurang baik karena biasanya plot akan cenderung flat (mendatar/konstan) untuk periode yang cukup panjang (Elvani et al, 2016).

Penelitian yang dilakukan ini memiliki tujuan untuk menganalisis model Autoregressive Integrated Moving Average (ARIMA) terbaik yang digunakan untuk menentukan target inflasi pada beberapa waktu ke depan. Penentuan model terbaik pada penelitian ini menggunakan beberapa indikator, antara lain signifikansi model, nilai AIC dan BIC.

Analisis Data

Package dan Library

#Import Library 
library(readxl)
library(TTR)
library(forecast)
library(tseries)
library(TSA)
library(dynlm)
library(lmtest)
library(imputeTS)
library(stats)
library(MASS)
library(kableExtra)
library(padr)
library(astsa)
library(tfarima)
library(ggplot2)

Import Data

Data analisis yang digunakan adalah data sekunder bulanan yang bersumber dari website https://www.bi.go.id/id/statistik/indikator/data-inflasi.aspx mengenai data inflasi periode Januari 2012 - Oktober 2022. Data ini terdiri atas satu peubah numerik dengan jumlah data 130 amatan. Data ini dipilih untuk dianalisis lebih lanjut mengenai forecasting dengan Double Moving Average, Double Exponential Smoothing, dan Forecasting menggunakan model ARIMA. 
data <- read_excel("C:/Users/denan/OneDrive/Documents/Kuliah/Kampus/Semester 5/MPDW/Data Inflasi.xlsx", sheet = "Sheet2")
dataMinyak <- read_excel("C:/Users/denan/OneDrive/Documents/Kuliah/Kampus/Semester 5/MPDW/Harga Minyak.xlsx", sheet = "Sheet2")
#Mengubah data ke dalam bentuk time series
data.ts<-ts(data,frequency=12, start=2012)
dataMinyak.ts <- ts(dataMinyak,frequency=12,start=2012)

Eksplorasi Data

Time Series Plot

#Plot semua data
plot(data.ts,xlab ="Waktu", ylab = "Data Inflasi (Persen)", col="black", main = "Plot Data Deret Waktu Data Inflasi")
points(data.ts)

plot(dataMinyak.ts,xlab ="Waktu", ylab = "Rata-Rata Harga Minyak", col="black", main = "Plot Data Deret Waktu Data Harga Minyak")
points(dataMinyak.ts)

Berdasarkan Time Series Plot, Pola data inflasi yang terbentuk dari data terlihat tidak stasioner, baik dalam ragam maupun nilai tengah yang cenderung memiliki pola siklik. Plot ini terlihat mengalami penurunan yang signifikan ketika kondisi perekonomian pada masa pandemi tahun 2020.

Standarisasi pada Data Banyaknya Ekspor Minyak

Standarisasi data diperlukan untuk menyamakan perbandingan data-data yang ingin dimodelkan. Data yang ingin distandarisasi adalah harga mintak dunia agar satuannya sama dengan satuan inflasi(persentase). Metode standarisasi yang digunakan adalah Min Max Normalization dengan interval [0,1]. 
#Min Max Normalization
min = min(dataMinyak.ts)
max = max(dataMinyak.ts)
stddataMinyak1 = as.data.frame(scale(dataMinyak.ts, center = min, scale = max-min)) 
stddataMinyak = ts(as.data.frame(scale(dataMinyak.ts, center = min, scale = max-min)), frequency = 12, start = 2012)
stddataMinyak
##             Jan        Feb        Mar        Apr        May        Jun
## 2012 0.88928837 0.94730837 1.00000000 0.95741752 0.85839643 0.72031962
## 2013 0.86887983 0.89510658 0.84223447 0.80429291 0.80961289 0.81349987
## 2014 0.83783272 0.86605220 0.85792058 0.86646567 0.87521750 0.90271341
## 2015 0.26941166 0.34883280 0.32853820 0.37729348 0.42859850 0.41619433
## 2016 0.09030924 0.10323025 0.16845551 0.20370402 0.25731760 0.27540701
## 2017 0.33642863 0.34431906 0.30865708 0.32168145 0.29821690 0.25969506
## 2018 0.46705143 0.43845292 0.44575760 0.49358257 0.54151090 0.52648807
## 2019 0.36737014 0.41440262 0.44182961 0.49134292 0.47332242 0.40020674
## 2020 0.41950211 0.33391334 0.11535877 0.00000000 0.09651133 0.19033508
## 2021 0.33656646 0.40747696 0.44227754 0.43318115 0.46884314 0.52469636
## 2022 0.64994401 0.74941855 0.94430183 0.85144285 0.92052718 0.98981824
##             Jul        Aug        Sep        Oct        Nov        Dec
## 2012 0.78260776 0.87067245 0.88112635 0.85138633 0.82829958 0.82849858
## 2013 0.87050983 0.90048379 0.90668685 0.87225912 0.84330941 0.87282265
## 2014 0.87018692 0.81667672 0.77326212 0.67227151 0.57834439 0.40995779
## 2015 0.34418124 0.25476785 0.26086657 0.26786114 0.22813334 0.16053062
## 2016 0.23860798 0.24636058 0.24808338 0.29201482 0.25032303 0.32640193
## 2017 0.27509691 0.29873374 0.32981308 0.35017659 0.40199845 0.41495392
## 2018 0.53362047 0.51725385 0.56149539 0.57558791 0.42663451 0.34025325
## 2019 0.41795159 0.37860281 0.40310104 0.37450254 0.40685675 0.43735033
## 2020 0.21731415 0.23154449 0.20211905 0.19491774 0.21976053 0.28615729
## 2021 0.53999483 0.49434060 0.53499871 0.63075200 0.60856232 0.53568783
## 2022 0.86870531 0.77453700 0.69439228 0.71616849

Splitting Data (Data Training dan Testing)

Data periode Januari 2012 sampai Desember 2019 digunakan sebagai data training untuk membangun model dan data Januari 2020 sampai Oktober 2022 digunakan sebagai data testing untuk memeriksa keakuratan model dalam memprediksi inflasi. 
#Split Data
data.train <- ts(data[1:96,1])
data.test <- ts(data[97:130,1])
dataMinyak.train <- ts(stddataMinyak[1:96,1])
dataMinyak.test <- ts(stddataMinyak[97:130,1])

#Time Series Data
training.ts<-ts(data.train,frequency=12, start=2012)
testing.ts<-ts(data.test,frequency=12, start=2020)

Minyaktraining.ts<-ts(dataMinyak.train,frequency=12, start=2012)
Minyaktesting.ts<-ts(dataMinyak.test,frequency=12, start=2020)


#Plot Data
plot(training.ts, xlab ="Waktu", ylab = "Data Inflasi (Persen)", col="red", main = "Plot Data Training")
points(training.ts)

plot(testing.ts, xlab ="Waktu", ylab = "Data Inflasi (Persen)", col="red", main = "Plot Data Testing")
points(testing.ts)

plot(Minyaktraining.ts,xlab ="Waktu", ylab = "Rata-Rata Harga Minyak", col="black", main = "Plot Data Deret Waktu Data Harga Minyak")
points(Minyaktraining.ts)

plot(Minyaktesting.ts,xlab ="Waktu", ylab = "Rata-Rata Harga Minyak", col="black", main = "Plot Data Deret Waktu Data Harga Minyak")
points(Minyaktesting.ts)

## Plot Data Training dan testing Inflasi
ts.plot(data.ts, xlab = "Periode", ylab ="Data Inflasi (Persen)", 
        main = "Plot Data Training dan Data Testing")
lines(training.ts, col = "blue")
lines(testing.ts, col="Red")
legend(2020,9,c("Data Training","Data Testing"), 
       lty=8, col=c("blue","red"), cex=0.8)
abline(v=2020, col=c("black"), lty=1, lwd=1)

#Plot Training - Testing Data Minyak
ts.plot(stddataMinyak, xlab = "Periode", ylab ="Data Inflasi (Persen)", 
        main = "Plot Data Training dan Data Testing")
lines(Minyaktraining.ts, col = "blue")
lines(Minyaktesting.ts, col="Red")
legend(2020,9,c("Data Training","Data Testing"), 
       lty=8, col=c("blue","red"), cex=0.8)
abline(v=2020, col=c("black"), lty=1, lwd=1)

Smoothing Data

Double Moving Average

iterasi_dma <- function(x,min,max,metode=c("SSE","MSE","RMSE","MAD","MAPE")){
  tabledma <- data.frame()
  z <- max-min+1
  temp <- sqrt(z)
  a <- round(temp) 
  b=a
  if(a*b<z){
    b=b+1
  }
  par(mfrow=c(a,b))
  for(i in min:max) {
    df_ts_sma <- SMA(x,  n=i)
    df_ts_dma <- SMA(df_ts_sma, n=i)
    At <- 2*df_ts_sma-df_ts_dma
    Bt <- df_ts_sma-df_ts_dma
    dmasmoothing <- At+Bt
    ramal_dma <- c(NA, dmasmoothing)
    df_dma <- cbind(df_aktual=c(x,NA), dmasmoothing=c(dmasmoothing,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)
    tabledma <- rbind(tabledma,ak)
  }
  opt <- tabledma$n[which.min(tabledma[,metode])]
  return(list(Akurasi=tabledma,n.optimum=paste("Nilai n optimum dengan menggunakan metode Double Moving Average berada pada n =",opt),smoothing=dmasmoothing, ramalan=ramal_dma))
}

dma_iter<-iterasi_dma(data.train,2,30, metode= "MAPE")
dma_iter[[2]]
## [1] "Nilai n optimum dengan menggunakan metode Double Moving Average berada pada n = 2"
#Single Moving  Average dengan N=2
df_ts_sma <- SMA(training.ts,  n=2)

#Double Moving  Average dengan N=2
df_ts_dma <- SMA(df_ts_sma, n=2)
At <- c(2*df_ts_sma-df_ts_dma)
Bt <- c(df_ts_sma-df_ts_dma)
pemulusan_dma <- c(At+Bt)
ramal_dma <- c(NA, pemulusan_dma)
df_dma <- cbind(df_aktual=c(training.ts,NA), pemulusan_dma=c(pemulusan_dma,NA), ramal_dma)
head(df_dma)
##      df_aktual pemulusan_dma ramal_dma
## [1,]      3.65            NA        NA
## [2,]      3.56            NA        NA
## [3,]      3.97         3.925        NA
## [4,]      4.50         4.705     3.925
## [5,]      4.45         4.715     4.705
## [6,]      4.53         4.505     4.715
#Plot
ts.plot(training.ts, xlab="Tahun", ylab="Data Inflasi", col="blue", lty=3, ylim=c(1, 12))
points(training.ts)

pemulusan_dma.ts<-ts(pemulusan_dma,frequency=12, start=2012)
ramal_dma.ts <- ts(ramal_dma, frequency=12, start=2012)

lines(pemulusan_dma.ts,col="red",lwd=2)
lines(ramal_dma.ts,col="black",lwd=2)
title("Double Moving Average N = 2", cex.main=1, font.main=4 ,col.main="black")
legend("topleft", c("Data Aktual","Pemulusan","Ramalan"),lty=3:1,col=c ("blue","red","black"))

### Ukuran Keakuratan Ramalan

# Ukuran Keakuratan
error.dma <- df_dma[, 1] - df_dma[, 3]

## SSE (Sum Square Error)
SSE.dma <- sum(error.dma^2, na.rm = T)

## MSE (Mean Squared Error)
MSE.dma <- mean(error.dma^2, na.rm = T)

## RMSE (Root Mean Square Error)
RMSE.dma <- sqrt(mean(error.dma^2, na.rm = T))

## MAD (Mean Absolute Deviation)
MAD.dma <- mean(abs(error.dma), na.rm = T)

## MAPE (Mean Absolute Percentage Error)
r.error.dma <- (error.dma/df_dma[, 1])*100 # Relative Error
MAPE.dma <- mean(abs(r.error.dma), na.rm = T)


akurasidma <- data.frame(
          "Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"), 
          "Double Moving Average N=2" = c(SSE.dma, MSE.dma, MAPE.dma, RMSE.dma, MAD.dma))
kable(akurasidma)
Ukuran.Keakuratan Double.Moving.Average.N.2
SSE 48.4388000
MSE 0.5208473
MAPE 10.4800038
RMSE 0.7216975
MAD 0.4926882
Diperoleh nilai MAPE pada data training sebesar 10.48%. Nilai MAPE yang diperoleh menghasilkan angka yang berada di rentang 10% - 20%, sehingga metode peramalan yang telah dilakukan pada data training dengan menggunakan metode Double Moving Average baik untuk digunakan. (Febrian et al. 2020)

Double Exponential Smoothing

Data yang digunakan memiliki pola yang tidak stasioner atau mengandung bentuk trend, maka digunakan metode Double Exponential Smoothing (DES) untuk pemulusan data. (Kristanti, Darsyah. 2017). Penerapan metode Double Exponential Smoothing menggunakan package Holt Winters 
Des <- HoltWinters(training.ts, gamma = FALSE)
## Warning in HoltWinters(training.ts, gamma = FALSE): optimization difficulties:
## ERROR: ABNORMAL_TERMINATION_IN_LNSRCH
Des
## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = training.ts, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 1
##  beta : 0.003543673
##  gamma: FALSE
## 
## Coefficients:
##          [,1]
## a  2.72000000
## b -0.06851632
tabledes <- data.frame("Alpha Optimum"=c(Des$alpha), "Beta Optimum"=c(Des$beta))
kable(tabledes)
Alpha.Optimum Beta.Optimum
alpha 1 0.0035437
Dengan menggunakan metode HoltWinters, diperoleh \(\alpha\) sebesar 1 dan \(\beta\) beta sebesar 0.0035437. Karena nilai \(\alpha\) pada metode Double Exponential Smoothing tidak dapat bernilai 1, maka digunakan angka terdekatnya, yaitu alpha = 0.999. (Habsari et al. 2020) 
Des1 <- HoltWinters(training.ts, gamma = FALSE, alpha = 0.999)
Des1
## Holt-Winters exponential smoothing with trend and without seasonal component.
## 
## Call:
## HoltWinters(x = training.ts, alpha = 0.999, gamma = FALSE)
## 
## Smoothing parameters:
##  alpha: 0.999
##  beta : 0.003504798
##  gamma: FALSE
## 
## Coefficients:
##          [,1]
## a  2.72021211
## b -0.06869918
tabledes1 <- data.frame("Alpha Optimum"=c(Des1$alpha), "Beta Optimum"=c(Des1$beta))
kable(tabledes1)
Alpha.Optimum Beta.Optimum
0.999 0.0035048

Dengan menggunakan metode HoltWinters, diperoleh \(\alpha\) sebesar 0.999 dan \(\beta\) beta sebesar 0.0035048

ramalanDES <- forecast(Des1, h=34)
ramalanDES
##          Point Forecast       Lo 80    Hi 80       Lo 95    Hi 95
## Jan 2020      2.6515129  1.87166493 3.431361  1.45883834 3.844188
## Feb 2020      2.5828137  1.47856197 3.687066  0.89400640 4.271621
## Mar 2020      2.5141146  1.15954531 3.868684  0.44247965 4.585749
## Apr 2020      2.4454154  0.87868774 4.012143  0.04931222 4.841519
## May 2020      2.3767162  0.62208810 4.131344 -0.30675591 5.060188
## Jun 2020      2.3080170  0.38262729 4.233407 -0.63661247 5.252647
## Jul 2020      2.2393179  0.15608288 4.322553 -0.94671510 5.425351
## Aug 2020      2.1706187 -0.06029182 4.401529 -1.24126450 5.582502
## Sep 2020      2.1019195 -0.26840149 4.472240 -1.52317363 5.727013
## Oct 2020      2.0332203 -0.46963061 4.536071 -1.79455986 5.861000
## Nov 2020      1.9645211 -0.66502255 4.594065 -2.05701890 5.986061
## Dec 2020      1.8958220 -0.85538656 4.647030 -2.31178838 6.103432
## Jan 2021      1.8271228 -1.04136507 4.695611 -2.55985081 6.214096
## Feb 2021      1.7584236 -1.22347807 4.740325 -2.80200147 6.318849
## Mar 2021      1.6897244 -1.40215339 4.781602 -3.03889464 6.418343
## Apr 2021      1.6210252 -1.57774794 4.819798 -3.27107619 6.513127
## May 2021      1.5523261 -1.75056312 4.855215 -3.49900704 6.603659
## Jun 2021      1.4836269 -1.92085612 4.888110 -3.72308056 6.690334
## Jul 2021      1.4149277 -2.08884846 4.918704 -3.94363553 6.773491
## Aug 2021      1.3462285 -2.25473255 4.947190 -4.16096620 6.853423
## Sep 2021      1.2775294 -2.41867670 4.973735 -4.37533000 6.930389
## Oct 2021      1.2088302 -2.58082915 4.998489 -4.58695362 7.004614
## Nov 2021      1.1401310 -2.74132124 5.021583 -4.79603794 7.076300
## Dec 2021      1.0714318 -2.90026996 5.043134 -5.00276189 7.145626
## Jan 2022      1.0027326 -3.05778008 5.063245 -5.20728567 7.212751
## Feb 2022      0.9340335 -3.21394578 5.082013 -5.40975334 7.277820
## Mar 2022      0.8653343 -3.36885214 5.099521 -5.61029504 7.340964
## Apr 2022      0.7966351 -3.52257630 5.115847 -5.80902871 7.402299
## May 2022      0.7279359 -3.67518844 5.131060 -6.00606169 7.461934
## Jun 2022      0.6592367 -3.82675263 5.145226 -6.20149195 7.519965
## Jul 2022      0.5905376 -3.97732752 5.158403 -6.39540922 7.576484
## Aug 2022      0.5218384 -4.12696699 5.170644 -6.58789590 7.631573
## Sep 2022      0.4531392 -4.27572065 5.181999 -6.77902785 7.685306
## Oct 2022      0.3844400 -4.42363430 5.192514 -6.96887511 7.737755
#Plot DES
plot(Des1)
points(training.ts)
legend("topright", c("Data aktual","Pemulusan DES"),lty=1:2,col=c ("black","red"))
Berdasarkan pada plot diatas, dapat dilihat bahwa data prediksi yang disimbolkan dengan garis merah mendekati data aktual yang disimbolkan dengan garis hitam

Keakuratan Ramalan Pada Data Training

SSEDES<-Des1$SSE 

MSEDES<-Des1$SSE/length(data.train) 

RMSEDES<-sqrt(MSEDES)

# Mencari MAPE dengan manual
fittedDES<-ramalanDES$fitted 
sisaanDES<-ramalanDES$residuals 
residDES<-training.ts-ramalanDES$fitted 
MAPEDES <-sum(abs(sisaanDES[3:length(data.train)]/data.train[3:length(data.train)])*100)/length(data.train)

akurasiDES <- data.frame(
          "Ukuran Keakuratan Data Training" = c("SSE", "MSE", "RMSE", "MAPE"), 
          "Double Exponential Smoothing" = c(SSEDES, MSEDES, RMSEDES, MAPEDES))
kable(akurasiDES)
Ukuran.Keakuratan.Data.Training Double.Exponential.Smoothing
SSE 34.8311655
MSE 0.3628246
RMSE 0.6023493
MAPE 8.0359336
Diperoleh nilai MAPE pada data training sebesar 8.04%. Nilai MAPE yang diperoleh menghasilkan angka yang lebih kecil dari 10%, sehingga metode peramalan yang telah dilakukan pada data training dengan menggunakan metode Double Exponential Smoothing sangat baik untuk digunakan. (Marlianah et al. 2019)

Perbandingan keakuratan pada data testing Double Moving Average dengan Double Exponential Smoothing

akurasiDESDMA <- data.frame(
          "Ukuran Keakuratan Data Training" = c("SSE", "MSE", "RMSE", "MAPE"), 
          "Double Exponential Smoothing" = c(SSEDES, MSEDES, RMSEDES, MAPEDES), "Double Moving Average N = 2" = c(SSE.dma, MSE.dma, RMSE.dma, MAPE.dma))
kable(akurasiDESDMA)
Ukuran.Keakuratan.Data.Training Double.Exponential.Smoothing Double.Moving.Average.N…2
SSE 34.8311655 48.4388000
MSE 0.3628246 0.5208473
RMSE 0.6023493 0.7216975
MAPE 8.0359336 10.4800038

Karena MAPE yang dihasilkan pada metode Double Exponential Smoothing lebih kecil dari MAPE yang dihasilkan oleh metode Double Moving Average, maka metode Double Exponential Smoothing lebih baik digunakan untuk dilakukan smoothing.

Keakuratan Ramalan Pada Data Testing

selisih.test <- c(ramalanDES$mean)-c(data.test)

SSEtesting<-sum(selisih.test^2) 
MSEtesting<-SSEtesting/length(data.test) 
RMSEtesting<-sqrt(MSEtesting)
MAPEtesting<-sum(abs(selisih.test))/length(data.test)

akurasi.test <- data.frame("Perbandingan Ukuran Keakuratan Data Testing" = c("SSE", "MSE", "RMSE", "MAPE"), "Double Exponential Smoothing" = c(SSEtesting, MSEtesting, RMSEtesting, MAPEtesting))
kable(akurasi.test)
Perbandingan.Ukuran.Keakuratan.Data.Testing Double.Exponential.Smoothing
SSE 133.974681
MSE 3.940432
RMSE 1.985052
MAPE 1.224675
Diperoleh nilai MAPE pada data testing sebesar 1.22%. Diperoleh nilai MAPE pada data testing sebesar 1.13%. Nilai MAPE yang diperoleh menghasilkan angka yang lebih kecil dari 10%, sehingga metode peramalan yang telah dilakukan pada data testing dengan menggunakan metode Dobule Exponential Smoothing sangat baik untuk digunakan. (Marlianah et al. 2019)

1. Mencari model ARIMA dari Variabel Independen (Harga Minyak)

Cek Kestasioneran Data

1. ACF & PACF Plot data inflasi dan data ekspor minyak rusia

acf(stddataMinyak, lag.max = 24, main="Plot ACF Data Inflasi")

pacf(stddataMinyak, lag.max = 24, main="Plot PACF Data Inflasi")

Berdasarkan Plot ACF Data minyak, nilai korelasi antar lag terlihat pada plot di atas menurun secara perlahan (tails off slowly). Hal tersebut mengindikasikan bahwa data minyak tidak stasioner.

2. Uji Formal (ADF-Test)

Secara formal, metode Augmented Dickey-Fuller (ADF) dapat memberikan hasil uji secara akurat untuk menentukan apakah sebuah data stasioner atau tidak. Namun, Uji ADF ini hanya mengukur tingkat stasioneritas berdasarkan nilai tengah saja. Dengan hipotesis yang diuji sebagai berikut :

\(H_{0}\): Data tidak stasioner

\(H_{1}\): Data stasioner

\(\alpha\) = 5% = 0.05

adf.test(dataMinyak.test)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dataMinyak.test
## Dickey-Fuller = -1.7649, Lag order = 3, p-value = 0.6641
## alternative hypothesis: stationary
kpss.test(dataMinyak.test)
## Warning in kpss.test(dataMinyak.test): p-value smaller than printed p-value
## 
##  KPSS Test for Level Stationarity
## 
## data:  dataMinyak.test
## KPSS Level = 0.82167, Truncation lag parameter = 3, p-value = 0.01
Berdasarkan hasil Augmented Dickey-Fuller Test (ADF Test) didapatkan p-value = 0.4858 > \(\alpha\), maka tak tolak \(H_{0}\). Artinya, cukup bukti untuk mengatakan bahwa data tidak stasioner pada taraf nyata 5%. Untuk mengatasi ketidakstasioneran tersebut, perlu dilakukan differencing.

Pemodelan Data dengan Data Training pada Data Inflasi

Karena data tidak stasioner, maka akan dilakukan differencing untuk menstasionerkan data.

Differencing 1

Berdasarkan plot diatas, setelah data dilakukan differencing sebanyak satu kali d=1, pola data inflasi periode Januari 2012 - Januari 2022 sudah stasioner dalam nilai tengah. 
data.dif1<-diff(data.train,differences = 1) 
plot.ts(data.dif1,lty=1,ylab= "Data Inflasi Pembedaan 1", main="Plot Differencing Data Inflasi")

Cek Kestasioneran Data Setelah Differencing

Pengujian menggunakan Augmented Dickey-Fuller Test

\(H_{0}\): Data tidak stasioner

\(H_{1}\): Data stasioner

\(\alpha\) = 5% = 0.05

adf.test(data.dif1)
## Warning in adf.test(data.dif1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.dif1
## Dickey-Fuller = -4.54, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
Berdasarkan hasil Augmented Dickey-Fuller Test (ADF Test) didapatkan p-value = 0.01 < \(\alpha\), maka tolak H_0. Artinya, cukup bukti untuk mengatakan bahwa data stasioner pada taraf nyata 5% setelah dilakukan differencing sebanyak 1 kali.

ACF & PACF Plot, EACF Matrix

acfdif <- acf(data.dif1, lag.max = 24)

pacfdif <- pacf(data.dif1, lag.max = 24)

plot(acfdif, main = "Plot ACF Data Inflasi Setelah Differencing satu kali")

plot(pacfdif, main = "Plot PACF Data Inflasi Setelah Differencing satu kali")

plot(acfdif, main = "Plot ACF Data Inflasi Setelah Differencing satu kali")

plot(pacfdif, main = "Plot PACF Data Inflasi Setelah Differencing satu kali")

Berdasarkan plot ACF dan PACF di atas,terlihat bahwa nilai korelasi antara data dengan lag seperti gambar di atas tidak turun secara perlahan, dimana pada plot ACF diperoleh cuts off pada lag ke-1 dan pada plot PACF diperoleh cut off di lag 1. Berdasarkan hasil eksplorasi di atas, model yang dapat dibentuk secara berurutan adalah ARIMA(0,1,1) dan ARIMA(1,1,0).
Berdasarkan plot ACF dan PACF di atas,terlihat bahwa nilai korelasi antara data dengan lag seperti gambar di atas tidak turun secara perlahan, dimana pada plot ACF diperoleh cuts off pada lag ke-1 dan pada plot PACF diperoleh cut off di lag 2. Berdasarkan hasil eksplorasi di atas, model yang dapat dibentuk secara berurutan adalah ARIMA(0,1,1) dan ARIMA(2,1,0). 
eacf(data.dif1)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o o o o o  x  o  o 
## 1 x x o o o o o o o o o  x  o  o 
## 2 o o o o o o o o o o o  x  o  o 
## 3 x o o o o o o o o o o  x  o  o 
## 4 x o o o o o o o o o o  x  o  o 
## 5 x x x o o o o o o o o  x  o  o 
## 6 x x x o o o o o o o o  x  o  o 
## 7 x o x o o o o o o o o  x  o  o

Dari Matriks EACF Data Inflasi dapat diduga model yang cocok adalah model :

  1. ARIMA(0,1,2)

  2. ARIMA (0,1,3)

  3. ARIMA (1,1,2)

  4. ARIMA (3,1,1)

  5. ARIMA (3,1,2)

Pemodelan Data dengan Data Training pada Data Harga Minyak

Karena data tidak stasioner, maka akan dilakukan differencing untuk menstasionerkan data.

dataMinyak.dif1<-diff(dataMinyak.train,differences = 1) 
plot.ts(dataMinyak.dif1,lty=1,ylab= "Data harga minyak Pembedaan 1", main="Plot Differencing Train Data Inflasi")

adf.test(dataMinyak.dif1)
## Warning in adf.test(dataMinyak.dif1): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  dataMinyak.dif1
## Dickey-Fuller = -4.207, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
acfMinyak <-acf(dataMinyak.dif1, lag.max = 24)

pacfMinyak <-pacf(dataMinyak.dif1, lag.max = 24)

plot(acfMinyak, main = "Plot ACF Harga Minyak Setelah Differencing satu kali")

plot(pacfMinyak,main = "Plot PACF Harga Minyak Setelah Differencing satu kali")

eacf(dataMinyak.dif1)
## AR/MA
##   0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x o o o o o o o o o o  o  o  o 
## 1 x o o o o o o o o o o  o  o  o 
## 2 x x o o o o o o o o o  o  o  o 
## 3 o x o o o o o o o o o  o  o  o 
## 4 x x o o o o o o o o o  o  o  o 
## 5 o x o 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 o o o o  o  o  o

Dari Matriks EACF Harga Minyak dapat diduga model yang cocok adalah model :

  1. ARIMA(0,1,2)

  2. ARIMA (1,1,2)

  3. ARIMA (2,1,2)

  4. ARIMA (3,1,2)

Pemodelan ARIMA Data Inflasi

Identifikasi Model Tentatif

Berdasarkan plot ACF, PACF, dan matriks EACF, diperoleh 5 model tentatif beserta orde parameternya, sebagai berikut:

ARIMA(0,1,1), dengan orde p = 0, d = 1, dan q = 1

ARIMA(2,1,0), dengan orde p = 2, d = 1, dan q = 0

ARIMA(0,1,2), dengan orde p = 0, d = 1, dan q = 2

ARIMA(0,1,3), dengan orde p = 0, d = 1, dan q = 3

ARIMA(1,1,2), dengan orde p = 1, d = 1, dan q = 2

ARIMA(3,1,1), dengan orde p = 3, d = 1, dan q = 1

ARIMA(3,1,2), dengan orde p = 3, d = 1, dan q = 2

model1 <- Arima(data.dif1, order=c(0,0,1), method="ML")   
model2 <- Arima(data.dif1, order=c(2,0,0), method="ML")  
model3 <- Arima(data.dif1, order=c(0,0,2), method="ML")
model4 <- Arima(data.dif1, order=c(0,0,3), method="ML")
model5 <- Arima(data.dif1, order=c(1,0,2), method="ML") 
model6 <- Arima(data.dif1, order=c(3,0,1), method="ML")
model7 <- Arima(data.dif1, order=c(3,0,2), method="ML")

Pendugaan Parameter Model

ARIMA(0,1,1)

summary(model1) 
## Series: data.dif1 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##          ma1     mean
##       0.3490  -0.0113
## s.e.  0.1007   0.0790
## 
## sigma^2 = 0.3345:  log likelihood = -81.83
## AIC=169.66   AICc=169.93   BIC=177.32
## 
## Training set error measures:
##                        ME      RMSE       MAE  MPE MAPE      MASE        ACF1
## Training set 0.0004262922 0.5722137 0.3513354 -Inf  Inf 0.7426473 -0.03000317
lmtest::coeftest((model1)) 
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ma1        0.348971   0.100691  3.4658 0.0005287 ***
## intercept -0.011280   0.078982 -0.1428 0.8864335    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diperoleh persamaan pada model:

\(\nabla Y_t\) = et - 0.348971 et-1 Untuk t = 1,2,3,…, 84
Berdasarkan data di atas diperoleh bahwa setiap parameter pada model ARIMA(0,1,1) signifikan. Hal ini dapat dilihat dari nilai Pr(>|z|) < 0.05 pada setiap parameter, sehingga model tersebut dapat dimasukan kedalam kandidat model terbaik.

ARIMA(2,1,0)

summary(model2) 
## Series: data.dif1 
## ARIMA(2,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2     mean
##       0.3020  -0.2419  -0.0108
## s.e.  0.0992   0.0987   0.0617
## 
## sigma^2 = 0.3288:  log likelihood = -80.53
## AIC=169.05   AICc=169.5   BIC=179.27
## 
## Training set error measures:
##                        ME      RMSE       MAE MPE MAPE     MASE          ACF1
## Training set 0.0009839349 0.5642494 0.3459638 Inf  Inf 0.731293 -0.0004607976
lmtest::coeftest((model2)) 
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value Pr(>|z|)   
## ar1        0.302001   0.099202  3.0443 0.002332 **
## ar2       -0.241863   0.098727 -2.4498 0.014293 * 
## intercept -0.010841   0.061723 -0.1756 0.860574   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diperoleh persamaan pada model:

\(\nabla Y_t\) = 0.302001 yt-1 - 0.241863 yt-2 + et Untuk t = 1,2,3,…, 84
Berdasarkan data di atas diperoleh bahwa setiap parameter pada model ARIMA(2,1,0) signifikan. Hal ini dapat dilihat dari nilai Pr(>|z|) < 0.05 pada setiap parameter, sehingga model tersebut dapat dimasukan kedalam kandidat model terbaik.

ARIMA(0,1,2)

summary(model3) 
## Series: data.dif1 
## ARIMA(0,0,2) with non-zero mean 
## 
## Coefficients:
##          ma1      ma2     mean
##       0.2843  -0.1095  -0.0106
## s.e.  0.1103   0.1086   0.0686
## 
## sigma^2 = 0.3346:  log likelihood = -81.33
## AIC=170.66   AICc=171.11   BIC=180.88
## 
## Training set error measures:
##                        ME      RMSE       MAE MPE MAPE      MASE       ACF1
## Training set 0.0003472967 0.5692085 0.3471446 Inf  Inf 0.7337888 0.01808154
lmtest::coeftest((model3)) 
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value Pr(>|z|)   
## ma1        0.284312   0.110346  2.5765  0.00998 **
## ma2       -0.109462   0.108619 -1.0078  0.31357   
## intercept -0.010646   0.068573 -0.1552  0.87663   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diperoleh persamaan pada model:

\(\nabla Y_t\) = et - 0.284312 et-1 + 0.109462 et-2 Untuk t = 1,2,3,…, 84
Berdasarkan data di atas diperoleh bahwa ada parameter pada model ARIMA(0,1,2) yang tidak signifikan. Hal ini dapat dilihat dari parameter yang memiliki nilai Pr(>|z|) > 0.05, sehingga model tersebut tidak dapat dimasukan kedalam kandidat model terbaik.

ARIMA(0,1,3)

summary(model4)
## Series: data.dif1 
## ARIMA(0,0,3) with non-zero mean 
## 
## Coefficients:
##          ma1      ma2      ma3     mean
##       0.3030  -0.1566  -0.1372  -0.0110
## s.e.  0.1023   0.1045   0.0994   0.0586
## 
## sigma^2 = 0.3315:  log likelihood = -80.41
## AIC=170.81   AICc=171.49   BIC=183.58
## 
## Training set error measures:
##                       ME      RMSE       MAE MPE MAPE      MASE         ACF1
## Training set 0.001378825 0.5634979 0.3481822 Inf  Inf 0.7359821 0.0006102943
lmtest::coeftest((model4))
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value Pr(>|z|)   
## ma1        0.302963   0.102343  2.9603 0.003074 **
## ma2       -0.156631   0.104480 -1.4991 0.133836   
## ma3       -0.137235   0.099393 -1.3807 0.167362   
## intercept -0.010982   0.058600 -0.1874 0.851346   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diperoleh persamaan pada model:

\(\nabla Y_t\) = et - 0.302963 et-1 + 0.156631 et-2 + 0.137235 et-3 Untuk t = 1,2,3,…, 84
Berdasarkan data di atas diperoleh bahwa ada parameter pada model ARIMA(0,1,3) yang tidak signifikan. Hal ini dapat dilihat dari parameter yang memiliki nilai Pr(>|z|) > 0.05, sehingga model tersebut tidak dapat dimasukan kedalam kandidat model terbaik.

ARIMA(1,1,2)

summary(model5)
## Series: data.dif1 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##          ar1      ma1      ma2     mean
##       0.8800  -0.6099  -0.3901  -0.0282
## s.e.  0.0599   0.1043   0.0995   0.0194
## 
## sigma^2 = 0.3221:  log likelihood = -79.98
## AIC=169.95   AICc=170.63   BIC=182.72
## 
## Training set error measures:
##                      ME      RMSE       MAE MPE MAPE      MASE        ACF1
## Training set 0.04241606 0.5554849 0.3527581 Inf  Inf 0.7456547 -0.01252043
lmtest::coeftest((model5))
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.879966   0.059945 14.6796 < 2.2e-16 ***
## ma1       -0.609942   0.104278 -5.8492 4.939e-09 ***
## ma2       -0.390054   0.099526 -3.9191 8.887e-05 ***
## intercept -0.028221   0.019389 -1.4555    0.1455    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Diperoleh persamaan pada model:

\(\nabla Y_t\) = 0.879966 yt-1 + et + 0.609942 et-1 + 0.390054 et-2 Untuk t = 1,2,3,…, 84
Berdasarkan data di atas diperoleh bahwa setiap parameter pada model ARIMA(1,1,2) signifikan. Hal ini dapat dilihat dari nilai Pr(>|z|) < 0.05 pada setiap parameter, sehingga model tersebut dapat dimasukan kedalam kandidat model terbaik.

ARIMA(3,1,1)

summary(model6) 
## Series: data.dif1 
## ARIMA(3,0,1) with non-zero mean 
## 
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): NaNs produced
##          ar1      ar2      ar3    ma1     mean
##       0.1284  -0.1871  -0.0517  0.171  -0.0109
## s.e.     NaN      NaN      NaN    NaN   0.0612
## 
## sigma^2 = 0.336:  log likelihood = -80.52
## AIC=173.04   AICc=173.99   BIC=188.36
## 
## Training set error measures:
##                       ME      RMSE       MAE MPE MAPE      MASE        ACF1
## Training set 0.001047293 0.5642167 0.3458751 Inf  Inf 0.7311055 0.001935606
lmtest::coeftest((model6)) 
## Warning in sqrt(diag(se)): NaNs produced
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value Pr(>|z|)
## ar1        0.128449        NaN     NaN      NaN
## ar2       -0.187099        NaN     NaN      NaN
## ar3       -0.051727        NaN     NaN      NaN
## ma1        0.170953        NaN     NaN      NaN
## intercept -0.010853   0.061186 -0.1774   0.8592

Diperoleh persamaan pada model:

\(\nabla Y_t\) = 0.128449 yt-1 - 0.187099 yt-2 - 0.051727 yt-3 + et - 0.170953 et-1 Untuk t = 1,2,3,…, 84
Berdasarkan data di atas diperoleh bahwa setiap parameter pada model ARIMA(3,1,1) tidak signifikan. Hal ini dapat dilihat dari parameter yang memiliki nilai Pr(>|z|) = NaN, sehingga model tersebut tidak dapat dimasukan kedalam kandidat model terbaik.

ARIMA(3,1,2)

summary(model7) 
## Series: data.dif1 
## ARIMA(3,0,2) with non-zero mean 
## 
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): NaNs produced
##          ar1      ar2      ar3     ma1      ma2     mean
##       0.1672  -0.1532  -0.0571  0.1311  -0.0464  -0.0108
## s.e.     NaN      NaN      NaN     NaN   0.1879   0.0604
## 
## sigma^2 = 0.3397:  log likelihood = -80.51
## AIC=175.02   AICc=176.31   BIC=192.9
## 
## Training set error measures:
##                       ME      RMSE       MAE MPE MAPE      MASE        ACF1
## Training set 0.001129378 0.5641619 0.3463942 Inf  Inf 0.7322027 0.002545374
lmtest::coeftest((model7)) 
## Warning in sqrt(diag(se)): NaNs produced
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value Pr(>|z|)
## ar1        0.167204        NaN     NaN      NaN
## ar2       -0.153203        NaN     NaN      NaN
## ar3       -0.057106        NaN     NaN      NaN
## ma1        0.131094        NaN     NaN      NaN
## ma2       -0.046401   0.187872 -0.2470   0.8049
## intercept -0.010843   0.060359 -0.1796   0.8574

Diperoleh persamaan pada model:

\(\nabla Y_t\) = 0.167204 yt-1 - 0.153203 yt-2 - 0.057106 yt-3 + et + 0.131094 et-1 + 0.046401 et-2 Untuk t = 1,2,3,…, 84
Berdasarkan data di atas diperoleh bahwa setiap parameter pada model ARIMA(3,1,2) tidak signifikan. Hal ini dapat dilihat dari parameter yang memiliki nilai Pr(>|z|) > 0.05, sehingga model tersebut tidak dapat dimasukan kedalam kandidat model terbaik.

Perbandingan Kebaikan Model Tentatif

#AIC ARIMA dan Signifikansi Parameter
modelaccuracy<-data.frame(
  "Model"=c("ARIMA(0,1,1)", "ARIMA(2,1,0)", "ARIMA(0,1,2)", "ARIMA(0,1,3)","ARIMA(1,1,2)", "ARIMA(3,1,1)", "ARIMA(3,1,2)"),
  "AIC"=c(model1$aic, model2$aic, model3$aic,model4$aic,model5$aic,model6$aic,model7$aic),
  "BIC"=c(model1$bic, model2$bic, model3$bic,model4$bic,model5$bic,model6$bic,model7$bic),
  "Signifikansi"=c("Signifikan","Signifikan","Tidak Signifikan","Tidak Signifikan","Signifikan","Tidak Signifikan","Tidak Signifikan"))

modelaccuracy
##          Model      AIC      BIC     Signifikansi
## 1 ARIMA(0,1,1) 169.6621 177.3237       Signifikan
## 2 ARIMA(2,1,0) 169.0506 179.2661       Signifikan
## 3 ARIMA(0,1,2) 170.6633 180.8788 Tidak Signifikan
## 4 ARIMA(0,1,3) 170.8122 183.5816 Tidak Signifikan
## 5 ARIMA(1,1,2) 169.9512 182.7206       Signifikan
## 6 ARIMA(3,1,1) 173.0399 188.3632 Tidak Signifikan
## 7 ARIMA(3,1,2) 175.0214 192.8985 Tidak Signifikan
kable(modelaccuracy)
Model AIC BIC Signifikansi
ARIMA(0,1,1) 169.6621 177.3237 Signifikan
ARIMA(2,1,0) 169.0506 179.2661 Signifikan
ARIMA(0,1,2) 170.6633 180.8788 Tidak Signifikan
ARIMA(0,1,3) 170.8122 183.5816 Tidak Signifikan
ARIMA(1,1,2) 169.9512 182.7206 Signifikan
ARIMA(3,1,1) 173.0399 188.3632 Tidak Signifikan
ARIMA(3,1,2) 175.0214 192.8985 Tidak Signifikan
Selain nilai AIC terkecil, perlu dilihat dari signifikasi pendugaan parameter model. Terdapat tiga model tentatif dengan semua parameter yang berpengaruh signifikan, yaitu model ARIMA (0,1,1), ARIMA (2,1,0), dan ARIMA (1,1,2). Model ARIMA(2,1,0) memiliki nilai AIC terendah dibandingkan ketiga kandidat model yang lainnya, sehingga model ARIMA(2,1,0) dipilih sebagai model terbaik

Pemodelan Arima Data Harga Minyak

Identifikasi Model Tentatif

modelm1 <- Arima(dataMinyak.dif1, order=c(1,0,0), method="ML")   
modelm2 <- Arima(dataMinyak.dif1, order=c(0,0,1), method="ML")  
modelm3 <- Arima(dataMinyak.dif1, order=c(0,0,2), method="ML")
modelm4 <- Arima(dataMinyak.dif1, order=c(0,0,3), method="ML")
modelm5 <- Arima(dataMinyak.dif1, order=c(1,0,2), method="ML") 
modelm6 <- Arima(dataMinyak.dif1, order=c(2,0,2), method="ML")
modelm7 <- Arima(dataMinyak.dif1, order=c(3,0,2), method="ML")
summary(modelm1) 
## Series: dataMinyak.dif1 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1     mean
##       0.3383  -0.0042
## s.e.  0.0969   0.0076
## 
## sigma^2 = 0.002473:  log likelihood = 151.27
## AIC=-296.54   AICc=-296.28   BIC=-288.88
## 
## Training set error measures:
##                         ME      RMSE        MAE      MPE     MAPE      MASE
## Training set -0.0002610692 0.0491984 0.04000873 169.6177 176.3294 0.8716311
##                    ACF1
## Training set 0.05174304
lmtest::coeftest((modelm1)) 
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value Pr(>|z|)    
## ar1        0.3382527  0.0968582  3.4922 0.000479 ***
## intercept -0.0042342  0.0075910 -0.5578 0.576985    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modelm2) 
## Series: dataMinyak.dif1 
## ARIMA(0,0,1) with non-zero mean 
## 
## Coefficients:
##          ma1     mean
##       0.3476  -0.0045
## s.e.  0.0876   0.0068
## 
## sigma^2 = 0.00245:  log likelihood = 151.71
## AIC=-297.42   AICc=-297.16   BIC=-289.76
## 
## Training set error measures:
##                         ME       RMSE        MAE      MPE     MAPE      MASE
## Training set -0.0001275517 0.04896879 0.03991593 155.6702 163.5104 0.8696092
##                    ACF1
## Training set 0.02560585
lmtest::coeftest((modelm2)) 
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## ma1        0.3476041  0.0875754  3.9692 7.211e-05 ***
## intercept -0.0045477  0.0067530 -0.6734    0.5007    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modelm3) 
## Series: dataMinyak.dif1 
## ARIMA(0,0,2) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2     mean
##       0.3999  0.1117  -0.0043
## s.e.  0.1080  0.1223   0.0075
## 
## sigma^2 = 0.002455:  log likelihood = 152.11
## AIC=-296.21   AICc=-295.77   BIC=-285.99
## 
## Training set error measures:
##                         ME       RMSE        MAE      MPE    MAPE      MASE
## Training set -0.0002182282 0.04875703 0.03944992 162.7797 172.368 0.8594567
##                     ACF1
## Training set -0.01320808
lmtest::coeftest((modelm3)) 
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## ma1        0.3999414  0.1079526  3.7048 0.0002116 ***
## ma2        0.1116879  0.1223251  0.9130 0.3612205    
## intercept -0.0043161  0.0075335 -0.5729 0.5667045    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modelm4)
## Series: dataMinyak.dif1 
## ARIMA(0,0,3) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2      ma3     mean
##       0.3584  0.0392  -0.1281  -0.0046
## s.e.  0.1050  0.1220   0.1054   0.0063
## 
## sigma^2 = 0.002445:  log likelihood = 152.81
## AIC=-295.62   AICc=-294.95   BIC=-282.85
## 
## Training set error measures:
##                         ME       RMSE       MAE      MPE     MAPE      MASE
## Training set -9.110006e-05 0.04839069 0.0388858 133.9833 146.9447 0.8471668
##                    ACF1
## Training set 0.01255828
lmtest::coeftest((modelm4))
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## ma1        0.3584091  0.1049821  3.4140 0.0006402 ***
## ma2        0.0391680  0.1220011  0.3210 0.7481755    
## ma3       -0.1281451  0.1053850 -1.2160 0.2239960    
## intercept -0.0045766  0.0063063 -0.7257 0.4680140    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modelm5)
## Series: dataMinyak.dif1 
## ARIMA(1,0,2) with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     ma2     mean
##       -0.3254  0.7193  0.2367  -0.0043
## s.e.   0.4794  0.4648  0.1808   0.0073
## 
## sigma^2 = 0.002471:  log likelihood = 152.31
## AIC=-294.63   AICc=-293.95   BIC=-281.86
## 
## Training set error measures:
##                         ME       RMSE        MAE      MPE     MAPE      MASE
## Training set -0.0001987893 0.04864788 0.03921654 155.5108 166.2217 0.8543724
##                     ACF1
## Training set -0.01014127
lmtest::coeftest((modelm5))
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value Pr(>|z|)
## ar1       -0.3254473  0.4793546 -0.6789   0.4972
## ma1        0.7192886  0.4647894  1.5476   0.1217
## ma2        0.2367229  0.1808414  1.3090   0.1905
## intercept -0.0043370  0.0073405 -0.5908   0.5546
summary(modelm6) 
## Series: dataMinyak.dif1 
## ARIMA(2,0,2) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2      ma1     ma2     mean
##       0.6013  -0.5547  -0.2406  0.3873  -0.0049
## s.e.  0.3117   0.2358   0.3316  0.2666   0.0059
## 
## sigma^2 = 0.002424:  log likelihood = 153.69
## AIC=-295.38   AICc=-294.42   BIC=-280.05
## 
## Training set error measures:
##                        ME       RMSE        MAE      MPE     MAPE      MASE
## Training set 0.0001048468 0.04792316 0.03825521 144.3942 163.2184 0.8334289
##                     ACF1
## Training set 0.007191712
lmtest::coeftest((modelm6)) 
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value Pr(>|z|)  
## ar1        0.6013382  0.3117021  1.9292  0.05371 .
## ar2       -0.5547321  0.2357832 -2.3527  0.01864 *
## ma1       -0.2405728  0.3316190 -0.7254  0.46818  
## ma2        0.3872516  0.2665591  1.4528  0.14628  
## intercept -0.0048618  0.0059245 -0.8206  0.41186  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(modelm7) 
## Series: dataMinyak.dif1 
## ARIMA(3,0,2) with non-zero mean 
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2     mean
##       0.9748  -1.1119  0.2345  -0.6169  0.8397  -0.0047
## s.e.  0.2279   0.2489  0.1625   0.2105  0.2062   0.0066
## 
## sigma^2 = 0.00242:  log likelihood = 154.09
## AIC=-294.18   AICc=-292.89   BIC=-276.3
## 
## Training set error measures:
##                        ME      RMSE        MAE      MPE    MAPE     MASE
## Training set 0.0001366945 0.0476176 0.03764542 102.7592 139.714 0.820144
##                     ACF1
## Training set 0.008919317
lmtest::coeftest((modelm7)) 
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.9747967  0.2278528  4.2782 1.884e-05 ***
## ar2       -1.1119163  0.2488990 -4.4673 7.920e-06 ***
## ar3        0.2344746  0.1624871  1.4430  0.149011    
## ma1       -0.6168930  0.2105051 -2.9305  0.003384 ** 
## ma2        0.8397036  0.2061955  4.0724 4.654e-05 ***
## intercept -0.0047314  0.0066123 -0.7155  0.474275    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#AIC ARIMA dan Signifikansi Parameter
modelmaccuracy<-data.frame(
  "Model"=c("ARIMA(1,1,0)", "ARIMA(0,1,1)", "ARIMA(0,1,2)", "ARIMA(0,1,3)","ARIMA(1,1,2)", "ARIMA(2,1,2)", "ARIMA(3,1,2)"),
  "AIC"=c(modelm1$aic, modelm2$aic, modelm3$aic,modelm4$aic,modelm5$aic,modelm6$aic,modelm7$aic),
  "BIC"=c(modelm1$bic, modelm2$bic, modelm3$bic,modelm4$bic,modelm5$bic,modelm6$bic,modelm7$bic),
  "Signifikansi"=c("Signifikan","Signifikan","Tidak Signifikan","Tidak Signifikan","Tidak Signifikan","Tidak Signifikan","Tidak Signifikan"))

modelmaccuracy
##          Model       AIC       BIC     Signifikansi
## 1 ARIMA(1,1,0) -296.5401 -288.8784       Signifikan
## 2 ARIMA(0,1,1) -297.4216 -289.7600       Signifikan
## 3 ARIMA(0,1,2) -296.2101 -285.9946 Tidak Signifikan
## 4 ARIMA(0,1,3) -295.6192 -282.8498 Tidak Signifikan
## 5 ARIMA(1,1,2) -294.6284 -281.8590 Tidak Signifikan
## 6 ARIMA(2,1,2) -295.3774 -280.0541 Tidak Signifikan
## 7 ARIMA(3,1,2) -294.1789 -276.3017 Tidak Signifikan
kable(modelmaccuracy)
Model AIC BIC Signifikansi
ARIMA(1,1,0) -296.5401 -288.8784 Signifikan
ARIMA(0,1,1) -297.4216 -289.7600 Signifikan
ARIMA(0,1,2) -296.2101 -285.9946 Tidak Signifikan
ARIMA(0,1,3) -295.6192 -282.8498 Tidak Signifikan
ARIMA(1,1,2) -294.6284 -281.8590 Tidak Signifikan
ARIMA(2,1,2) -295.3774 -280.0541 Tidak Signifikan
ARIMA(3,1,2) -294.1789 -276.3017 Tidak Signifikan
Selain nilai AIC terkecil, perlu dilihat dari signifikasi pendugaan parameter model. Terdapat tiga model tentatif dengan semua parameter yang berpengaruh signifikan, yaitu model ARIMA (0,1,1) dan ARIMA (0,1,2). Model ARIMA(1,1,0) memiliki nilai AIC terendah dibandingkan kedua kandidat model yang lainnya, sehingga model ARIMA(1,1,0) dipilih sebagai model terbaik

Diagnostik Model

Overfitting

Overfitting yang digunakan adalah model ARIMA(2,1,0). Model ARIMA(2,1,0) akan dibandingkan dengan model ARIMA(3,1,0) 
model2b <- Arima(data.dif1, order=c(3,0,0), method="ML")
modelm1b <- Arima(dataMinyak.dif1, order=c(2,0,0), method="ML")

ARIMA(3,1,0)

coeftest(model2)
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value Pr(>|z|)   
## ar1        0.302001   0.099202  3.0443 0.002332 **
## ar2       -0.241863   0.098727 -2.4498 0.014293 * 
## intercept -0.010841   0.061723 -0.1756 0.860574   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest((model2b))
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value Pr(>|z|)   
## ar1        0.298835   0.102398  2.9184 0.003519 **
## ar2       -0.238035   0.103396 -2.3022 0.021326 * 
## ar3       -0.012674   0.101811 -0.1245 0.900928   
## intercept -0.010852   0.060960 -0.1780 0.858708   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Karena parameter \(\phi _3\) atau ar3 tidak berbeda secara signifikan dari angka nol, dan

karena parameter \(\phi _1\) atau ar1 dan parameter \(\phi _2\) atau ar2 pada model ARIMA(3,1,0) tidak berubah secara signifikan dari parameter \(\phi _1\) atau ar1 dan parameter \(\phi _2\) atau ar2 pada model ARIMA(2,1,0),

Maka model ARIMA(2,1,0) tidak terjadi overfitting atau merupakan model yang baik untuk digunakan. (Cryer dan Chan. 2008) 
coeftest(modelm1)
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value Pr(>|z|)    
## ar1        0.3382527  0.0968582  3.4922 0.000479 ***
## intercept -0.0042342  0.0075910 -0.5578 0.576985    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest((modelm1b))
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.3895350  0.1013394  3.8439 0.0001211 ***
## ar2       -0.1574882  0.1017088 -1.5484 0.1215208    
## intercept -0.0046498  0.0064898 -0.7165 0.4736917    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Uji Diagnostik

Diagnostik Model: Eksploratif

sisaanm1 <- modelm1$residuals

par(mfrow=c(2,2))
qqnorm(sisaanm1)
box(col="black",lwd=2)
qqline(sisaanm1, col = "red", lwd =1, col.lab="black",
       col.axis="black",col.sub = "black")
box(col="black",lwd=2)
plot(c(1:length(sisaanm1)),sisaanm1,col="black",col.lab="black",col.axis="black")
box(col="black",lwd=2)
acf(sisaanm1,col="black",col.sub = "black",col.axis="black", col.lab="black")
box(col="black",lwd=2)
pacf(sisaanm1,col="black",col.sub = "black",col.axis="black", col.lab="black",col.main="black")
box(col="black",lwd=2)

sisaan2 <- model2$residuals

par(mfrow=c(2,2))
qqnorm(sisaan2)
box(col="black",lwd=2)
qqline(sisaan2, col = "red", lwd =1, col.lab="black",
       col.axis="black",col.sub = "black")
box(col="black",lwd=2)
plot(c(1:length(sisaan2)),sisaan2,col="black",col.lab="black",col.axis="black")
box(col="black",lwd=2)
acf(sisaan2,col="black",col.sub = "black",col.axis="black", col.lab="black")
box(col="black",lwd=2)
pacf(sisaan2,col="black",col.sub = "black",col.axis="black", col.lab="black",col.main="black")
box(col="black",lwd=2)
Berdasarkan hasil eksplorasi menggunakan Q-Q plot, terlihat bahwa sisaan berdistribusi tidak mengikuti garis normal, sehingga dapat dikatakan bahwa sisaan tidak menyebar normal. Kemudian, plot sisaan yang diperoleh menunjukkan bahwa sisaan memiliki pola acak dan tersebar di sekitar nilai nol. Sedangkan pada plot ACF dan PACF, nilai awal amatan tidak melewati garis signifikan, atau dapat dikatakan bahwa sisaan saling bebas.

Diagnostik Model: Uji Formal

1. Sisaan Menyebar Normal

Uji formal ini dilakukan dengan Shapiro-Wilk test.

Hipotesis yang diuji:

\(H_{0}\): Sisaan menyebar normal

\(H_{1}\): Sisaan tidak menyebar normal

\(\alpha\) = 5% = 0.05

shapiro.test(sisaanm1)
## 
##  Shapiro-Wilk normality test
## 
## data:  sisaanm1
## W = 0.98058, p-value = 0.171
shapiro.test(sisaan2)
## 
##  Shapiro-Wilk normality test
## 
## data:  sisaan2
## W = 0.85335, p-value = 2.943e-08
Berdasarkan hasil Shapiro-Wilk test, diperoleh p-value = 2.943 x 10-07 > α = 0.05, maka tolak \(H_{0}\). Artinya, tidak cukup bukti untuk menyatakan bahwa sisaan menyebar normal pada taraf nyata 5%.

2. Sisaan Saling Bebas

Uji formal ini dilakukan dengan LJung-Box test.

Hipotesis yang diuji:

\(H_0\): Sisaan antara lag saling bebas

\(H_1\): Sisaan antara lag tidak saling bebas

\(\alpha\) = 5% = 0.05

Box.test(sisaanm1, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  sisaanm1
## X-squared = 0.26246, df = 1, p-value = 0.6084
Box.test(sisaan2, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  sisaan2
## X-squared = 2.0816e-05, df = 1, p-value = 0.9964
Berdasarkan LJung-Box test, diperoleh p-value = 0.9964 > α = 0.05, maka tak tolak \(H_0\). Artinya, cukup bukti untuk menyatakan bahwa sisaan antara lag saling bebas atau dapat dikatakan tidak ada autokorelasi antara sisaan lag pada taraf nyata 5%.

3. Nilai Tengah Sisaan Sama dengan Nol

Uji formal ini dilakukan dengan t test.

Hipotesis yang diuji:

\(H_0\): Nilai tengah sisaan sama dengan nol

\(H_1\): Nilai tengah sisaan tidak sama dengan nol

\(\alpha\) = 5% = 0.05

t.test(sisaanm1, mu = 0, conf.level = 0.95) 
## 
##  One Sample t-test
## 
## data:  sisaanm1
## t = -0.051449, df = 94, p-value = 0.9591
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.010336323  0.009814184
## sample estimates:
##     mean of x 
## -0.0002610692
t.test(sisaan2, mu = 0, conf.level = 0.95) 
## 
##  One Sample t-test
## 
## data:  sisaan2
## t = 0.016907, df = 94, p-value = 0.9865
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1145691  0.1165370
## sample estimates:
##    mean of x 
## 0.0009839349

Berdasarkan t-test, diperoleh p-value = 0.9865 > α = 0.05, maka tak tolak \(H_0\). Artinya, cukup bukti untuk menyatakan bahwa nilai tengah sisaan sama dengan nol pada taraf nyata 5%.

Uji diagnostik yang dipenuhi mengindikasikan bahwa ARIMA(1,1,0) merupakan model terbaik untuk Variabel Independen (harga minyak)[

Prewhitening

Proses Prewhitening diperlukan untuk menduga input yang akan digunakan untuk membentuk Plot CCF 
umx <- as.um(modelm1)
umx
##          mu         ar1        sig2 
## 0.004234244 0.338252672 0.002472536
a <- residuals(umx, training.ts, method = "cond")
b <- residuals(umx, Minyaktraining.ts, method = "cond")

CCF Plot

Cross Corelation Function Plot digunakan untuk melihat keeratan dan arah hubungan antara variabel independen dan variabel respon. Selain itu, plot CCF juga dapat digunakan untuk menduga parameter (b,r,s) pada model ARIMAX 
#CCF Plot
ccf(b,a, ylab = "CCF", main = "Cross Corelation Function")

Output dari plot CCF itu menandakan bahwa dari order ke cenderung menurun secara perlahan.oleh karena itu, parameter yang diduga dari plot CCF tersebut adalah (0, 1, 0)

Model ARIMAX

modelx = arimax(data.train, order = c(0,1,0), xreg = data.frame(dataMinyak.train))
modelx
## 
## Call:
## arimax(x = data.train, order = c(0, 1, 0), xreg = data.frame(dataMinyak.train))
## 
## Coefficients:
##       dataMinyak.train
##                -0.6327
## s.e.            1.1716
## 
## sigma^2 estimated as 0.3594:  log likelihood = -86.19,  aic = 174.38

Plot ACF, PACF, dan EACF dari residual model arimax

acf(residuals(modelx))

pacf(residuals(modelx))

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

Identifikasi Model Tentatif

Berdasarkan plot ACF, PACF, dan matriks EACF, diperoleh 5 model tentatif beserta orde parameternya, sebagai berikut:

ARIMA(0,1,1), dengan orde p = 0, d = 0, dan q = 1

ARIMA(2,1,0), dengan orde p = 2, d = 0, dan q = 0

ARIMA(0,1,2), dengan orde p = 0, d = 0, dan q = 2

ARIMA(0,1,3), dengan orde p = 0, d = 0, dan q = 3

ARIMA(1,1,2), dengan orde p = 1, d = 0, dan q = 2

ARIMA(3,1,1), dengan orde p = 3, d = 0, dan q = 1

ARIMA(3,1,2), dengan orde p = 3, d = 0, dan q = 2

modelax1 <- Arima(data.train, order=c(0,0,1), method="ML", xreg = c(dataMinyak.train))   
modelax2 <- Arima(data.train, order=c(2,0,0), method="ML", xreg = c(dataMinyak.train))  
modelax3 <- Arima(data.train, order=c(0,0,2), method="ML", xreg = c(dataMinyak.train))
modelax4 <- Arima(data.train, order=c(0,0,3), method="ML", xreg = c(dataMinyak.train))
modelax5 <- Arima(data.train, order=c(1,0,2), method="ML", xreg = c(dataMinyak.train)) 
modelax6 <- Arima(data.train, order=c(3,0,1), method="ML", xreg = c(dataMinyak.train))
modelax7 <- Arima(data.train, order=c(3,0,2), method="ML", xreg = c(dataMinyak.train))

Signifikansi Parameter Model Tentatif

coeftest(modelax1)
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ma1       0.914591   0.084783 10.7874 < 2.2e-16 ***
## intercept 3.400731   0.431478  7.8816 3.233e-15 ***
## xreg      2.403292   0.727335  3.3042 0.0009523 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(modelax2)
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1        1.201994   0.097284 12.3555 < 2.2e-16 ***
## ar2       -0.278384   0.099107 -2.8089  0.004971 ** 
## intercept  4.391639   0.946698  4.6389 3.503e-06 ***
## xreg       0.139219   1.209881  0.1151  0.908391    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(modelax3)
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ma1       1.295620   0.086609 14.9595 < 2.2e-16 ***
## ma2       0.628170   0.060653 10.3569 < 2.2e-16 ***
## intercept 3.408556   0.463610  7.3522  1.95e-13 ***
## xreg      2.361064   0.766347  3.0809  0.002064 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(modelax4)
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ma1       1.377524   0.095936 14.3588 < 2.2e-16 ***
## ma2       1.011006   0.121304  8.3345 < 2.2e-16 ***
## ma3       0.368366   0.077624  4.7455  2.08e-06 ***
## intercept 3.531901   0.523775  6.7432  1.55e-11 ***
## xreg      2.111399   0.857052  2.4636   0.01376 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(modelax5)
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.912335   0.050826 17.9503 < 2.2e-16 ***
## ma1        0.345260   0.127398  2.7101  0.006726 ** 
## ma2       -0.051901   0.121111 -0.4285  0.668259    
## intercept  4.395149   0.987333  4.4515 8.526e-06 ***
## xreg       0.053733   1.191065  0.0451  0.964017    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(modelax6)
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.290158   0.099858  2.9057  0.003664 ** 
## ar2        0.762640   0.073549 10.3691 < 2.2e-16 ***
## ar3       -0.193106   0.101896 -1.8951  0.058076 .  
## ma1        0.999994   0.030027 33.3030 < 2.2e-16 ***
## intercept  4.463833   0.943423  4.7315 2.228e-06 ***
## xreg      -0.017592   1.134679 -0.0155  0.987630    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(modelax7)
## 
## z test of coefficients:
## 
##             Estimate Std. Error z value  Pr(>|z|)    
## ar1        1.2170781  0.6894693  1.7652   0.07752 .  
## ar2       -0.4477967  0.6071507 -0.7375   0.46080    
## ar3        0.1684925  0.2465278  0.6835   0.49431    
## ma1        0.0460946  0.6940931  0.0664   0.94705    
## ma2       -0.0288537  0.3973332 -0.0726   0.94211    
## intercept  4.3828220  1.0356297  4.2320 2.316e-05 ***
## xreg       0.0031994  1.1763546  0.0027   0.99783    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#AIC ARIMA dan Signifikansi Parameter
modelaccuracy<-data.frame(
  "Model"=c("ARIMAX(0,1,0)(0,0,1)", "ARIMAX(0,1,0)(2,0,0)", "ARIMAX(0,1,0)(0,0,2)", "ARIMAX(0,1,0)(0,0,3)","ARIMAX(0,1,0)(1,0,2)", "ARIMAX(0,1,0)(3,0,1)", "ARIMAX(0,1,0)(3,0,2)"),
  "AIC"=c(modelax1$aic, modelax2$aic, modelax3$aic,modelax4$aic,modelax5$aic,modelax6$aic,modelax7$aic),
  "BIC"=c(modelax1$bic, modelax2$bic, modelax3$bic,modelax4$bic,modelax5$bic,modelax6$bic,modelax7$bic),
  "Signifikansi"=c("Signifikan","Tidak Signifikan","Signifikan","Signifikan","Tidak Signifikan","Tidak Signifikan","Tidak Signifikan"))

modelaccuracy
##                  Model      AIC      BIC     Signifikansi
## 1 ARIMAX(0,1,0)(0,0,1) 273.4669 283.7243       Signifikan
## 2 ARIMAX(0,1,0)(2,0,0) 175.7913 188.6131 Tidak Signifikan
## 3 ARIMAX(0,1,0)(0,0,2) 223.3576 236.1794       Signifikan
## 4 ARIMAX(0,1,0)(0,0,3) 207.4901 222.8762       Signifikan
## 5 ARIMAX(0,1,0)(1,0,2) 174.5559 189.9420 Tidak Signifikan
## 6 ARIMAX(0,1,0)(3,0,1) 174.8558 192.8062 Tidak Signifikan
## 7 ARIMAX(0,1,0)(3,0,2) 177.6442 198.1590 Tidak Signifikan
kable(modelaccuracy)
Model AIC BIC Signifikansi
ARIMAX(0,1,0)(0,0,1) 273.4669 283.7243 Signifikan
ARIMAX(0,1,0)(2,0,0) 175.7913 188.6131 Tidak Signifikan
ARIMAX(0,1,0)(0,0,2) 223.3576 236.1794 Signifikan
ARIMAX(0,1,0)(0,0,3) 207.4901 222.8762 Signifikan
ARIMAX(0,1,0)(1,0,2) 174.5559 189.9420 Tidak Signifikan
ARIMAX(0,1,0)(3,0,1) 174.8558 192.8062 Tidak Signifikan
ARIMAX(0,1,0)(3,0,2) 177.6442 198.1590 Tidak Signifikan
Selain nilai AIC terkecil, pemilihan model terbaik perlu dilihat dari signifikasi pendugaan parameter model. Terdapat tiga model tentatif dengan semua parameter yang berpengaruh signifikan, yaitu model ARIMAX(0,1,0)(0,0,1), ARIMAX(0,1,0)(0,0,2), dan Model ARIMAX(0,1,0)(0,0,3). Model ARIMAX(0,1,0)(0,0,3) memiliki nilai AIC terendah dibandingkan kedua kandidat model yang lainnya, sehingga model ARIMAX(0,1,0)(0,0,3) dipilih sebagai model terbaik

Overfitting

modelax41 = Arima(data.train, order=c(0,0,4), method="ML", xreg = c(dataMinyak.train))
modelax42 = Arima(data.train, order=c(1,0,3), method="ML", xreg = c(dataMinyak.train))

coeftest(modelax41)
## 
## z test of coefficients:
## 
##           Estimate Std. Error z value  Pr(>|z|)    
## ma1        1.31053    0.13395  9.7839 < 2.2e-16 ***
## ma2        1.03918    0.15899  6.5363 6.306e-11 ***
## ma3        0.64945    0.14732  4.4084 1.041e-05 ***
## ma4        0.16970    0.12134  1.3985   0.16197    
## intercept  3.68036    0.55483  6.6333 3.284e-11 ***
## xreg       1.81889    0.90400  2.0120   0.04421 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
coeftest(modelax42)
## 
## z test of coefficients:
## 
##            Estimate Std. Error z value  Pr(>|z|)    
## ar1        0.931435   0.044591 20.8885 < 2.2e-16 ***
## ma1        0.338832   0.113873  2.9755  0.002925 ** 
## ma2       -0.109392   0.118707 -0.9215  0.356775    
## ma3       -0.107932   0.104000 -1.0378  0.299360    
## intercept  4.347407   1.036381  4.1948 2.731e-05 ***
## xreg       0.063067   1.176254  0.0536  0.957240    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
intinya ga ada overfitting. ARIMAX(0,1,0)(0,0,3) tetap menjadi model terbaik.

Diagnostik Model : Eksploratif

sisaanax4 <- modelax4$residuals

par(mfrow=c(2,2))
qqnorm(sisaanax4)
box(col="black",lwd=2)
qqline(sisaanax4, col = "red", lwd =1, col.lab="black",
       col.axis="black",col.sub = "black")
box(col="black",lwd=2)
plot(c(1:length(sisaanax4)),sisaanax4,col="black",col.lab="black",col.axis="black")
box(col="black",lwd=2)
acf(sisaanax4,col="black",col.sub = "black",col.axis="black", col.lab="black")
box(col="black",lwd=2)
pacf(sisaanax4,col="black",col.sub = "black",col.axis="black", col.lab="black",col.main="black")
box(col="black",lwd=2)

Diagnostik Model: Uji Formal

1. Sisaan Menyebar Normal

Uji formal ini dilakukan dengan Shapiro-Wilk test.

Hipotesis yang diuji:

\(H_{0}\): Sisaan menyebar normal

\(H_{1}\): Sisaan tidak menyebar normal

\(\alpha\) = 5% = 0.05

shapiro.test(sisaanax4)
## 
##  Shapiro-Wilk normality test
## 
## data:  sisaanax4
## W = 0.93775, p-value = 0.0001941
Berdasarkan hasil Shapiro-Wilk test, diperoleh p-value = 2.943 x 10-07 > α = 0.05, maka tolak \(H_{0}\). Artinya, tidak cukup bukti untuk menyatakan bahwa sisaan menyebar normal pada taraf nyata 5%.

2. Sisaan Saling Bebas

Uji formal ini dilakukan dengan LJung-Box test.

Hipotesis yang diuji:

\(H_0\): Sisaan antara lag saling bebas

\(H_1\): Sisaan antara lag tidak saling bebas

\(\alpha\) = 5% = 0.05

Box.test(sisaanax4, type = "Ljung")
## 
##  Box-Ljung test
## 
## data:  sisaanax4
## X-squared = 1.5489, df = 1, p-value = 0.2133
Berdasarkan LJung-Box test, diperoleh p-value = 0.9964 > α = 0.05, maka tak tolak \(H_0\). Artinya, cukup bukti untuk menyatakan bahwa sisaan antara lag saling bebas atau dapat dikatakan tidak ada autokorelasi antara sisaan lag pada taraf nyata 5%.

3. Nilai Tengah Sisaan Sama dengan Nol

Uji formal ini dilakukan dengan t test.

Hipotesis yang diuji:

\(H_0\): Nilai tengah sisaan sama dengan nol

\(H_1\): Nilai tengah sisaan tidak sama dengan nol

\(\alpha\) = 5% = 0.05

t.test(sisaanax4, mu = 0, conf.level = 0.95) 
## 
##  One Sample t-test
## 
## data:  sisaanax4
## t = 0.098091, df = 95, p-value = 0.9221
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1283083  0.1416468
## sample estimates:
##   mean of x 
## 0.006669249

Berdasarkan t-test, diperoleh p-value = 0.9865 > α = 0.05, maka tak tolak \(H_0\). Artinya, cukup bukti untuk menyatakan bahwa nilai tengah sisaan sama dengan nol pada taraf nyata 5%.

Uji diagnostik yang dipenuhi mengindikasikan bahwa ARIMAX(0,1,0)(0,0,3) merupakan model terbaik.[

Forecasting

ramalan <- forecast(Arima(data.train, order=c(2,1,0),method="ML",include.drift = TRUE),h = length(data.test))
data.ramalan <- ramalan$mean
data.ramalan.ts <- ts(data.ramalan, start = 2020, frequency = 12)
plot(ramalan,col="black",col.sub ="black",col.axis="black",
     col.lab="black",col.main="black",lwd=2)
box(col="black",lwd=2)

ts.plot(data.ts, ylab = "Data Inflasi(Persen)", col="black",lwd=2,main="Forecasting ARIMA(2,1,0)",gpars = list(col.main="black",col.axis="black",col.sub="black"))
lines(data.ramalan.ts, col = "blue",lwd=2)
lines(testing.ts, col = "red", lwd =2)
legend(2018,9,c("Data Training", "Data Testing","Data Forecast ARIMA"), 
       lwd=2, col=c("black","red","blue"), cex=0.8)
box(col="black",lwd=2)
Berdasarkan plot di atas, terlihat bahwa hasil forecasting memiliki garis yang cenderung linear dan tidak menglami trend 
perbandingan.temp<-matrix(data=c(data.test[1:6], data.ramalan[1:6]), nrow = 6, ncol = 2)
colnames(perbandingan.temp)<-c("Aktual","Hasil Forecast")
head(perbandingan.temp)
##      Aktual Hasil Forecast
## [1,]   2.68       2.656693
## [2,]   2.98       2.695106
## [3,]   2.96       2.711829
## [4,]   2.67       2.697400
## [5,]   2.19       2.678808
## [6,]   1.96       2.666494
error <- data.frame(data.test)-data.frame(data.ramalan[1:34]) 

## SSE (Sum Square Error)
SSE <- sum(error^2, na.rm = T)

## MSE (Mean Squared Error)
MSE<- sapply(error^2, mean, na.rm = T)

## RMSE (Root Mean Square Error)
RMSE <- sqrt(MSE)

## MAD (Mean Absolute Deviation)
MAD <- sapply(abs(error), mean, na.rm = T)

## MAPE (Mean Absolute Percentage Error)
r.error <- (error/data.frame(data.test))*100 # Relative Error
MAPE <- sapply(abs(r.error), mean, na.rm = T)

akurasifarima <- data.frame(
  "Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"), 
  "Forecasting" = c(SSE, MSE, MAPE, RMSE, MAD))
akurasifarima
##   Ukuran.Keakuratan Forecasting
## 1               SSE   62.598323
## 2               MSE    1.841127
## 3              MAPE   49.245009
## 4              RMSE    1.356881
## 5               MAD    1.088177
kable(akurasifarima)
Ukuran.Keakuratan Forecasting
SSE 62.598323
MSE 1.841127
MAPE 49.245009
RMSE 1.356881
MAD 1.088177

Forecasting With ARIMAX

ramalan2 <- forecast(Arima(data.train, order=c(3,1,0), include.drift = TRUE,method="ML",,xreg=dataMinyak.train),xreg=dataMinyak.train)
data.ramalan2 <- ramalan2$mean
data.ramalan.ts2 <- ts(data.ramalan2, start = 2020, frequency = 12)
plot(ramalan2,col="black",col.sub ="black",col.axis="black",
     col.lab="black",col.main="black",lwd=2)
box(col="black",lwd=2)

ts.plot(data.ts, ylab = "Data Inflasi(Persen)", col="black",lwd=2,main="Forecasting ARIMA(2,1,0)",gpars = list(col.main="black",col.axis="black",col.sub="black"))
lines(data.ramalan.ts2, col = "blue",lwd=2)
lines(testing.ts, col = "red", lwd =2)
legend(2018,9,c("Data Training", "Data Testing","Data Forecast ARIMA"), 
       lwd=2, col=c("black","red","blue"), cex=0.8)
box(col="black",lwd=2)

perbandingan.temp<-matrix(data=c(data.test[1:6], data.ramalan2[1:6]), nrow = 6, ncol = 2)
colnames(perbandingan.temp)<-c("Aktual","Hasil Forecast")
head(perbandingan.temp)
##      Aktual Hasil Forecast
## [1,]   2.68       2.491559
## [2,]   2.98       2.505609
## [3,]   2.96       2.502117
## [4,]   2.67       2.502972
## [5,]   2.19       2.519057
## [6,]   1.96       2.555693
error <- data.frame(data.test)-data.frame(data.ramalan2[1:34]) 

## SSE (Sum Square Error)
SSE <- sum(error^2, na.rm = T)

## MSE (Mean Squared Error)
MSE<- sapply(error^2, mean, na.rm = T)

## RMSE (Root Mean Square Error)
RMSE <- sqrt(MSE)

## MAD (Mean Absolute Deviation)
MAD <- sapply(abs(error), mean, na.rm = T)

## MAPE (Mean Absolute Percentage Error)
r.error <- (error/data.frame(data.test))*100 # Relative Error
MAPE <- sapply(abs(r.error), mean, na.rm = T)

akurasifarima <- data.frame(
  "Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"), 
  "Forecasting" = c(SSE, MSE, MAPE, RMSE, MAD))
akurasifarima
##   Ukuran.Keakuratan Forecasting
## 1               SSE   62.963601
## 2               MSE    1.851871
## 3              MAPE   44.044558
## 4              RMSE    1.360835
## 5               MAD    1.041699
kable(akurasifarima)
Ukuran.Keakuratan Forecasting
SSE 62.963601
MSE 1.851871
MAPE 44.044558
RMSE 1.360834
MAD 1.041699
Berdasarkan data keakuratan diatas, dengan menggunakan fungsi forecasting pada model ARIMA, diperoleh MAPE sebesar 49.24%. Karena nilai MAPE yang sangat besar, forecasting dengan menggunakan model ARIMA pada data yang digunakan diduga kurang tepat untuk dilakukan
Berdasarkan data perbandingan nilai data aktual dan forecasting, terjadi perbedaan nilai yang cukup signifikan antara nilai data aktual dengan nilai data forecast (yang diramalkan). Pada data aktual terdapat trend, sedangkan pada data forecasting menghasilkan garis yang linear. Dan berdasarkan hasil uji keakuratan yang telah dilakukan, diperoleh nilai MAPE yang sangat tinggi. Sehingga, dapat disimpulkan forecasting yang dilakukan dengan menggunakan metode ARIMA kurang baik untuk dilakukan.

Kesimpulan

Data yang digunakan pada penelitian ini adalah data inflasi Indonesia tahun 2012 sampai tahun 2022 yang dihitung perbulan. Data tersebut merupakan data yang tidak stasioner atau mengandung trend dalam. Karena data yang tidak stasioner dan mengandung trend, maka dipilihlah metode Double Exponential Smoothing untuk melakukan pemulusan dan peramalan pada data. Data perlu distasionerkan agar bisa dilakukan peramalan dengan menggunakan metode ARIMA. Berdasarkan analisis yang telah dilakukan, pendugaan model terbaik pada data inflasi di Indonesia pada periode Januari 2012 - Januari 2022 adalah ARIMA (2,1,0) dilihat dari signifikasi pada semua pendugaan parameter model dan nilai AIC terkecil. Persamaan Model ARIMA (2,1,0), sebagai berikut:
\(\nabla Y_t\) = 0.302001 yt-1 - 0.241863 yt-2 + et Untuk t = 1,2,3,…, 84

Daftar Pustaka

Habsari HDP, Purnamasari I, Yuniarti D. 2020. Peramalan menggunakan metode Double Exponential Smoothing dan verifikasi hasil peramalan menggunakan grafik pengendali tracking signal. Jurnal Ilmu Matematika dan Terapan. 14(1): 13-22. doi:10.30598/barekengvol14iss1pp013-022.

Cryer JD, Chan KS. 2008. Time Series Analysis With Applications in R. New York (NY): Springer New York.

Febrian D, Idrus SIA, Nainggolan DAJ. 2020. The comparison of Double Moving Average and Double Exponential Smoothing methods in forecasting the number of foreign tourists coming to North Sumatera, Journal of Physics: Conference Series. 1462(1): 012046. DOI:10.1088/1742-6596/1462/1/012046

Kristanti N, Darsyah MY. 2017. Perbandingan peramalan metode Single Exponential Smoothing dan Double Exsponential Smoothing pada karakteristik penduduk bekerja di Indonesia tahun 2017.Di dalam: Prosiding Seminar Nasional Mahasiswa Unimus(Vol.1)

Hartati H. 2017. Penggunaan Metode Arima Dalam Meramal Pergerakan Inflasi. Jurnal Matematika Sains Dan Teknologi. 18(1): 1–10.

Elvani SP, Utary AR, Yudaruddin R. 2016. Peramalan jumlah produksi tanaman kelapa sawit dengan menggunakan metode Arima (Autoregressive Integrated Moving Average). Jurnal Manajemen. 8(1): 95-112. doi

Hikmah A, Agoestanto A, Arifudin R. 2018. Peramalan deret waktu dengan menggunakan Autoregressive (AR), jaringan syaraf tiruan Radial basis function (Rbf) dan Hibrid Ar-Rbf pada Inflasi Indonesia. Unnes Journal of Mathematics. 7(2): 228-241. doi:10.15294/ujm.v7i2.1398

Hakimah M, Rahmawati WM, Afandi AY. 2020. PENGUKURAN KINERJA METODE PERAMALAN TIPE EXPONENTIAL SMOOTHING DALAM PARAMETER TERBAIKNYA. Jurnal Ilmiah NERO. 5(1): 44-50. - Tinjauan pustaka (SES & DES) - https://nero.trunojoyo.ac.id/index.php/nero/article/view/150/139