Data merupakan data sekunder yang diperoleh dari website resmi badan pusat statistika (BPS). Di sini diambil data inflasi bahan makanan pada tahun 2015 sampai 2019.

Langkah 1. install packages

library(readr)
## Warning: package 'readr' was built under R version 3.6.3
library(tsoutliers)
## Warning: package 'tsoutliers' was built under R version 3.6.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(forecast)
## Warning: package 'forecast' was built under R version 3.6.3
library(tseries)
## Warning: package 'tseries' was built under R version 3.6.3
library(lmtest)
## Warning: package 'lmtest' was built under R version 3.6.3
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(FitAR)
## Warning: package 'FitAR' was built under R version 3.6.3
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 3.6.3
## Loading required package: leaps
## Warning: package 'leaps' was built under R version 3.6.3
## Loading required package: ltsa
## Loading required package: bestglm
## Warning: package 'bestglm' was built under R version 3.6.3
## 
## Attaching package: 'FitAR'
## The following object is masked from 'package:forecast':
## 
##     BoxCox

Langkah 2. Input Data dan definisikan data menjadi data time series

##INPUT DATA
data=read.csv("bahanmakanan.csv")
head(data)
##    X0.6
## 1 -1.47
## 2 -0.73
## 3 -0.79
## 4  1.39
## 5  1.60
## 6  2.02
data=data
min(data)
## [1] -1.97
##UBAH DATA MENJADI FORMAT TIME SERIES (+ PLOT DATA)
Yt <- ts(data, start = c(2015,1), freq=12)
Yt
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov   Dec
## 2015 -1.47 -0.73 -0.79  1.39  1.60  2.02  0.91 -1.07 -1.06  0.33  3.20  2.20
## 2016 -0.58  0.69 -0.94  0.30  1.62  1.12 -0.68 -0.07 -0.21  1.66  0.50  0.66
## 2017 -0.31 -0.66 -1.13  0.86  0.69  0.21 -0.67 -0.53 -0.45  0.37  2.26  2.34
## 2018  0.13  0.14 -0.26  0.21  0.88  0.86 -1.10 -1.62  0.15  0.24  1.45  0.92
## 2019 -1.11 -0.01  1.45  2.02  1.63  0.80 -0.19 -1.97 -0.41  0.37  0.78

Langkah 3. Indentifikasi Data

plot(Yt)

Berdasarkan plot data diatas, plot time series dari data inflasi bahan makanan diatas dapat diketahui secara visual bahwa data inflasi bahan makanan dari bulan Januari 2015 hingga Desember 2019 belum stasioner. Namun, kestasioneran data tidak dapat dapat ditentukan hanya dengan menggunakan plot. Oleh karena itu, perlu dilakukan uji kestasioneran data yang meliputi uji kestasioneran dalam varians, dan juga uji kestasioneran dalam mean.Uji kestasioneran dalam varians dilakukan dengan melihat nilai lambda menggunakan metode Box-Cox transformation. Sedangkan uji kestasioneran data dalam mean dilakukan dengan menggunakan Augmented Dickey Fuller (ADF) test. PLot ACF dan PACF sebelum dilakukan differencing

##PLOT ACF PACF
par(mfrow=c(2,1))
acf(Yt, lag=100)
pacf(Yt, lag=100)

##UJI STASIONER TERHADAP VARIANS (BOX.COX)
lambda=BoxCox.lambda(Yt)
## Warning in guerrero(x, lower, upper): Guerrero's method for selecting a Box-Cox
## parameter (lambda) is given for strictly positive data.
lambda
## [1] 0.8943683

Pada pengujian Box-Cox transformation diperoleh nilai lambda sebesar 0.89, sehingga dapat disimpulkan bahwa data sudah stasioner secara varians karena nilai lambda sudah mendekati 1, sehingga tidak perlu dilakukan transformasi data. Tahap selanjutnya yaitu uji kestasioneran data dalam mean dan diperoleh hasil bahwa data juga belum stasioner dalam mean. Dapat dilihat pada ACF dan PACF nya, pada pola ACF masih berbentuk seperti musiman. Maka dari itu dilakukan differencing untuk menghilangkan pola musiman.

##DIFFERENCING
Yt1=diff(Yt, lag=12, differences = 1)
adf.test(Yt1)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  Yt1
## Dickey-Fuller = -3.5849, Lag order = 3, p-value = 0.04441
## alternative hypothesis: stationary
acf(Yt1, lag=100)

pacf(Yt1, lag=100)

Setelah dilakukan proses differencing orde ke-1 untuk data seasonal dan diperoleh bahwa pola ACF masih berpola musiman.

Yt2=diff(Yt1, differences = 1)
adf.test(Yt2)
## Warning in adf.test(Yt2): p-value smaller than printed p-value
## 
##  Augmented Dickey-Fuller Test
## 
## data:  Yt2
## Dickey-Fuller = -4.6227, Lag order = 3, p-value = 0.01
## alternative hypothesis: stationary

Karena ACF masih berpola musiman sehingga dilakukan differencing orde ke-1 untuk data yang sudah dilakukan differencing orde ke-1 untuk data seasonal.

Setelah dilakukan differencing orde ke-1 untuk data yang sudah dilakukan differencing orde ke-1 untuk data seasonal, dilakukan kembali adf.test dan p-value yang dihasilkan adalah 0.01 yang artinya data sudah stasioner dalam mean yang dapat dilihat pada plot berikut

plot(Yt2)

Setelah memperoleh data yang stasioner dalam varians dan mean, selanjutnya melakukan pendugaan model awal menggunakan plot ACF dan PACF. Plot ACF dan PACF inflasi bahan makanan setelah differencing orde ke-1 dan orde ke-2 dapat di lihat di gambar berikut

par(mfrow=c(2,1))
acf(Yt2, lag=100)
pacf(Yt2, lag=100)

Dari gambar diatas dapat dilihat bahwa plot ACF dan PACF stasioner. Pada Plot ACF mengalami cut off lag 2 pada data non musiman dan cut off lag 1 pada data musiman sedangkan dengan PACF mengalami cut off lag 2. Dari gambar di atas juga diketahui bahwa data tidak mengandung unsur musiman.

Langkah 4. Estimasi parameter untuk mengetahui apakah model- model sementara tersebut sudah signifikan. Sehingga dapat dilanjutkan untuk diagnostic checking. Dikatakan signifikan Ketika nilai p-value < 0.05.

##PENAKSIRAN PARAMETER
fitS1=Arima(Yt2, order=c(0,1,1),seasonal=list(order=c(0,1,0),period=12))
fitS2=Arima(Yt2, order=c(0,1,2),seasonal=list(order=c(0,1,0),period=12))
fitS3=Arima(Yt2, order=c(1,1,0),seasonal=list(order=c(0,1,0),period=12))
fitS4=Arima(Yt2, order=c(1,1,1),seasonal=list(order=c(0,1,0),period=12))
fitS5=Arima(Yt2, order=c(1,1,2),seasonal=list(order=c(0,1,0),period=12))
fitS6=Arima(Yt2, order=c(2,1,0),seasonal=list(order=c(0,1,0),period=12))
fitS7=Arima(Yt2, order=c(2,1,1),seasonal=list(order=c(0,1,0),period=12))
fitS8=Arima(Yt2, order=c(2,1,2),seasonal=list(order=c(0,1,0),period=12))
fitS9=Arima(Yt2, order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))
fitS10=Arima(Yt2, order=c(0,1,2),seasonal=list(order=c(0,1,1),period=12))
fitS11=Arima(Yt2, order=c(1,1,0),seasonal=list(order=c(0,1,1),period=12))
fitS12=Arima(Yt2, order=c(1,1,1),seasonal=list(order=c(0,1,1),period=12))
fitS13=Arima(Yt2, order=c(1,1,2),seasonal=list(order=c(0,1,1),period=12))
fitS14=Arima(Yt2, order=c(2,1,0),seasonal=list(order=c(0,1,1),period=12))
fitS15=Arima(Yt2, order=c(2,1,1),seasonal=list(order=c(0,1,1),period=12))
fitS16=Arima(Yt2, order=c(2,1,2),seasonal=list(order=c(0,1,1),period=12))
i=AIC(fitS1,fitS2,fitS3,fitS4,fitS5,fitS6,fitS7,fitS8,fitS9,fitS10,fitS11,fitS12,fitS13,fitS14,fitS15,fitS16)
min(i)
## [1] 2
##UJI SIGNIFIKANSI PARAMETER
printstatarima=function(Ytr,digits=4,se=T,...){
  if (length(Ytr$coef)>0){
    cat("\nCoefficients:\n")
    coef=round(Ytr$coef,digits=digits)
    if(se && nrow(Ytr$var.coef)){
      ses=rep(0, length(coef))
      ses[Ytr$mask] = round(sqrt(diag(Ytr$var.coef)),digits=digits)
      coef=matrix(coef,1,dimnames=list(NULL, names(coef)))
      coef=rbind(coef,s.e.=ses)
      statt=coef[1,]/ses
      pval=2*pt(abs(statt),df=length(Ytr$residuals)-1,lower.tail=F)
      coef=rbind(coef,t=round(statt,digits=digits),sign.=round(pval,digits=digits))
      coef=t(coef)
    }
    print.default(coef,print.gap=2)
  }
}
printstatarima(fitS1)#
## 
## Coefficients:
##            s.e.         t  sign.
## ma1  -1  0.0765  -13.0719      0
printstatarima(fitS2)#
## 
## Coefficients:
##                 s.e.        t   sign.
## ma1  -1.7834  0.2650  -6.7298  0.0000
## ma2   0.7834  0.2467   3.1755  0.0027
printstatarima(fitS3)
## 
## Coefficients:
##                s.e.       t   sign.
## ar1  -0.435  0.1587  -2.741  0.0088
printstatarima(fitS4)
## 
## Coefficients:
##                 s.e.         t   sign.
## ar1  -0.2331  0.1697   -1.3736  0.1764
## ma1  -1.0000  0.0795  -12.5786  0.0000
printstatarima(fitS5)
## 
## Coefficients:
##                 s.e.        t   sign.
## ar1  -0.7818  0.5627  -1.3894  0.1716
## ma1  -0.1539  0.5036  -0.3056  0.7613
## ma2  -0.8461  0.5020  -1.6855  0.0988
printstatarima(fitS6)#
## 
## Coefficients:
##                 s.e.        t  sign.
## ar1  -0.7093  0.1364  -5.2001      0
## ar2  -0.6297  0.1315  -4.7886      0
printstatarima(fitS7)
## 
## Coefficients:
##                 s.e.         t   sign.
## ar1  -0.3563  0.1508   -2.3627  0.0225
## ar2  -0.4859  0.1494   -3.2523  0.0022
## ma1  -1.0000  0.0912  -10.9649  0.0000
printstatarima(fitS8)
## 
## Coefficients:
##                 s.e.        t   sign.
## ar1   0.1840  0.3051   0.6031  0.5495
## ar2  -0.2825  0.2252  -1.2544  0.2162
## ma1  -1.7513  0.4196  -4.1737  0.0001
## ma2   0.7513  0.4086   1.8387  0.0726
printstatarima(fitS9)
## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -1.0000  0.1003  -9.9701  0.0000
## sma1  -0.9997  0.4953  -2.0184  0.0495
printstatarima(fitS10)
## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -1.8238  0.4171  -4.3726  0.0001
## ma2    0.8238  0.3873   2.1270  0.0389
## sma1  -0.9997  0.7404  -1.3502  0.1837
printstatarima(fitS11)
## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.4794  0.1535  -3.1231  0.0031
## sma1  -1.0000  0.4176  -2.3946  0.0209
printstatarima(fitS12)
## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.2454  0.1669  -1.4703  0.1484
## ma1   -1.0000  0.1071  -9.3371  0.0000
## sma1  -0.9999  0.5307  -1.8841  0.0660
printstatarima(fitS13)
## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.7727  0.7393  -1.0452  0.3015
## ma1   -0.2088  0.7255  -0.2878  0.7748
## ma2   -0.7912  0.7229  -1.0945  0.2796
## sma1  -1.0000  0.4990  -2.0040  0.0511
printstatarima(fitS14)
## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.7542  0.1382  -5.4573  0.0000
## ar2   -0.5947  0.1436  -4.1414  0.0001
## sma1  -1.0000  0.5083  -1.9673  0.0553
printstatarima(fitS15)
## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.3656  0.1562  -2.3406  0.0237
## ar2   -0.4281  0.1543  -2.7745  0.0080
## ma1   -1.0000  0.1269  -7.8802  0.0000
## sma1  -1.0000  0.7017  -1.4251  0.1610
printstatarima(fitS16)
## 
## Coefficients:
##                  s.e.        t   sign.
## ar1    0.3310  0.2183   1.5163  0.1364
## ar2   -0.1640  0.1946  -0.8428  0.4038
## ma1   -1.9894  0.3584  -5.5508  0.0000
## ma2    0.9998  0.3590   2.7850  0.0078
## sma1  -1.0000  0.5318  -1.8804  0.0665

diketahui bahwa hanya model SARIMA (2,1,0)〖(0,1,0)〗^12, SARIMA (0,1,1)〖(0,1,0)〗^12, SARIMA (1,1,1)〖(0,1,0)〗^12, SARIMA (0,1,2)〖(0,1,0)〗^12, dan SARIMA (0,1,2)〖(0,1,1)〗^12 yang parameternya signifikan (P-value berada dibawah level toleransi (α = 0,05). Dengan demikian model tersebut memenuhi syarat signifikansi parameter.

Langkah 5. Diagnostic Checking a. Uji Non Autokorelasi Setelah estimasi parameter, tahap selanjutnya adalah pemeriksaan diagnostik model. Pada tahap ini akan diuji apakah model sudah layak atau belum. Kelayakan tersebut dinilai dengan pengujian asumsi non autokorelasi. Ketika p-value > 0.05 maka asumsi white noise terpenuhi.

##UJI ASUMSI RESIDUAL
#Pengujian Non Autokorelasi
resY1=residuals(fitS1)
wnY=Box.test(resY1, type=c("Ljung-Box"))
wnY
## 
##  Box-Ljung test
## 
## data:  resY1
## X-squared = 2.2625, df = 1, p-value = 0.1325
resY2=residuals(fitS2)
wnY=Box.test(resY2, type=c("Ljung-Box"))
wnY
## 
##  Box-Ljung test
## 
## data:  resY2
## X-squared = 1.3264, df = 1, p-value = 0.2495
resY3=residuals(fitS6)
wnY=Box.test(resY3, type=c("Ljung-Box"))
wnY
## 
##  Box-Ljung test
## 
## data:  resY3
## X-squared = 1.7319, df = 1, p-value = 0.1882

Hasil Uji Ljung-Box, dapat diketahui bahwa model SARIMA (0,1,1)〖(0,1,0)〗^(12 )tidak memenuhi asumsi non autokorelasi sedangkan SARIMA (2,1,0)〖(0,1,0)〗^12, SARIMA (1,1,1)〖(0,1,0)〗^12, SARIMA (0,1,2)〖(0,1,0)〗^12, dan SARIMA (0,1,2)〖(0,1,1)〗^12 memenuhi asumsi non autokorelasi, dikarenakan p-value>5% sehingga H0 diterima.

  1. Uji Kenormalan Residual Uji kenormalan residual dilakukan dengan uji Kolmogorov-smirnov, dimana uji ini dilakukan untuk menguji kenormalan residual. Dikatakan normal ketika nilai p-value> α.
#Pengujian Normalitas
nY1=length(resY1)
sdY1=sd(resY1)
resnY1=rnorm(nY1,0,sdY1)
ks.test(resY1,resnY1)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  resY1 and resnY1
## D = 0.26087, p-value = 0.08726
## alternative hypothesis: two-sided
nY2=length(resY2)
sdY2=sd(resY2)
resnY2=rnorm(nY2,0,sdY2)
ks.test(resY2,resnY2)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  resY2 and resnY2
## D = 0.26087, p-value = 0.08726
## alternative hypothesis: two-sided
nY3=length(resY2)
sdY3=sd(resY3)
resnY3=rnorm(nY3,0,sdY3)
ks.test(resY3,resnY3)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  resY3 and resnY3
## D = 0.21739, p-value = 0.2287
## alternative hypothesis: two-sided

Dari tabel dapat diketahui hasil uji Kolmogorov-Smirnov model SARIMA (0,1,1)〖(0,1,0)〗^(12 )tidak memenuhi asumsi kenormalan residual sedangkan SARIMA (2,1,0)〖(0,1,0)〗^12, SARIMA (1,1,1)〖(0,1,0)〗^12, SARIMA (0,1,2)〖(0,1,0)〗^12, dan SARIMA (0,1,2)〖(0,1,1)〗^12 memenuhi asumsi kenormalan residual karena p-value lebih besar dari 0.05.

  1. Uji Homoskedastisitas Untuk menguji homoskdastisitas, maka dilakukan uji Q-Ljung (Q-statistics). Dapat disimpulkan bahwa residual bersifat homoskdastisitas ketika nilaip-value> α., maka asumsi homoskedastisitas terpenuhi.
#Pengujian Homoskedastisitas
Box.test(resY1^2, type=c("Ljung-Box"))
## 
##  Box-Ljung test
## 
## data:  resY1^2
## X-squared = 0.10337, df = 1, p-value = 0.7478
Box.test(resY2^2, type=c("Ljung-Box"))
## 
##  Box-Ljung test
## 
## data:  resY2^2
## X-squared = 1.2473, df = 1, p-value = 0.2641
Box.test(resY3^2, type=c("Ljung-Box"))
## 
##  Box-Ljung test
## 
## data:  resY3^2
## X-squared = 0.73235, df = 1, p-value = 0.3921

Dari tabel dapat diketahui hasil uji Homoskedastisitas, model SARIMA (0,1,1)〖(0,1,0)〗^(12 )tidak memenuhi asumsi homoskedastisitas sedangkan SARIMA (2,1,0)〖(0,1,0)〗^12, SARIMA (1,1,1)〖(0,1,0)〗^12, SARIMA (0,1,2)〖(0,1,0)〗^12, dan SARIMA (0,1,2)〖(0,1,1)〗^12 memenuhi memenuhi asumsi Homoskedastisitas karena p-value lebih besar dari 0.05.

Langkah 6. PEMILIHAN MODEL TERBAIK Setelah melalui proses estimasi parameter dan diagnostic checking dapat diketahui bahwa model SARIMA (2,1,0)〖(0,1,0)〗^12, SARIMA (1,1,1)〖(0,1,0)〗^12, SARIMA (0,1,2)〖(0,1,0)〗^12, dan SARIMA (0,1,2)〖(0,1,1)〗^12 menjadi model yang bisa dipakai dalam meramalkan jumlah data inflasi bahan makanan. Karena telah memenuhi asumsi non autokorelasi, residual normal, dan homoskedastisitas.

##HASIL (CEK MAPE)
summary(fitS1)
## Series: Yt2 
## ARIMA(0,1,1)(0,1,0)[12] 
## 
## Coefficients:
##           ma1
##       -1.0000
## s.e.   0.0765
## 
## sigma^2 estimated as 4.942:  log likelihood=-74.44
## AIC=152.89   AICc=153.29   BIC=155.88
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 0.02191867 1.854134 1.255512 289.7789 369.2069 0.7208273
##                    ACF1
## Training set -0.2147347
summary(fitS2)
## Series: Yt2 
## ARIMA(0,1,2)(0,1,0)[12] 
## 
## Coefficients:
##           ma1     ma2
##       -1.7834  0.7834
## s.e.   0.2650  0.2467
## 
## sigma^2 estimated as 3.613:  log likelihood=-70.63
## AIC=147.27   AICc=148.1   BIC=151.76
## 
## Training set error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set -0.1605507 1.560295 1.024553 13.00653 317.1335 0.5882268 0.1644129
summary(fitS6)
## Series: Yt2 
## ARIMA(2,1,0)(0,1,0)[12] 
## 
## Coefficients:
##           ar1      ar2
##       -0.7093  -0.6297
## s.e.   0.1364   0.1315
## 
## sigma^2 estimated as 6.091:  log likelihood=-76.21
## AIC=158.43   AICc=159.26   BIC=162.92
## 
## Training set error measures:
##                      ME     RMSE      MAE     MPE     MAPE      MASE       ACF1
## Training set 0.03621568 2.025959 1.200467 174.058 273.5792 0.6892244 -0.1878758

Model SARIMA (0,1,2)〖(0,1,1)〗^12dianggap model yang layak karena parameter-parameter yang ada di dalamnya telah signifikan serta residual-residualnya telah mengandung asumsi non autokorelasi dan berdistribusi normal.

Langkah 7. Plot Peramalan

Setelah dilakukan pengecekan diagnostik dan semua pengujian menunjukkan kesesuaian model, maka dari model umum ARIMA yang terbentuk tersebut dapat dilakukan peramalan atau forecasting. Secara matematis model SARIMA (0,1,2)〖(0,1,1)〗^12dapat dituliskan dalam bentuk seperti berikut ini:

##PLOT PERAMALAN
fc=forecast(fitS2,h = 12)
fc
##          Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## Dec 2019     -0.5511610 -3.031934 1.929612 -4.345177 3.242855
## Jan 2020      0.2000458 -2.895862 3.295953 -4.534736 4.934828
## Feb 2020      1.1100458 -1.985862 4.205953 -3.624736 5.844828
## Mar 2020      1.8800458 -1.215862 4.975953 -2.854736 6.614828
## Apr 2020      0.1200458 -2.975862 3.215953 -4.614736 4.854828
## May 2020     -1.0399542 -4.135862 2.055953 -5.774736 3.694828
## Jun 2020     -0.7899542 -3.885862 2.305953 -5.524736 3.944828
## Jul 2020      0.9900458 -2.105862 4.085953 -3.744736 5.724828
## Aug 2020     -1.2399542 -4.335862 1.855953 -5.974736 3.494828
## Sep 2020     -0.1899542 -3.285862 2.905953 -4.924736 4.544828
## Oct 2020      0.7100458 -2.385862 3.805953 -4.024736 5.444828
## Nov 2020     -0.7799542 -3.875862 2.315953 -5.514736 3.954828
par(mfrow=c(1,1))
plot(fc)
fc1=fitted(fitS2)
lines(fc1,col='blue')

Dari gambar diatas dapat dilihat bahwa plot peramalan inflasi bahan makanan mendekati data aktualnya. Sehingga, model SARIMA (0,1,2)〖(0,1,1)〗^12 sudah baik digunakan untuk meramalkan inflasi bahan makanan tahun 2020.