MPDW BARU BANGET ARIMAX

2022-11-29

Package dan Library

library(readxl)
library(TTR)
library(forecast)

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

library(tseries)
library(TSA)

## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast

## 
## Attaching package: 'TSA'

## The following objects are masked from 'package:stats':
## 
##     acf, arima

## The following object is masked from 'package:utils':
## 
##     tar

library(dynlm)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(lmtest)
library(imputeTS)

## 
## Attaching package: 'imputeTS'

## The following object is masked from 'package:zoo':
## 
##     na.locf

## The following object is masked from 'package:tseries':
## 
##     na.remove

library(stats)
library(MASS)
library(kableExtra)

## Warning in !is.null(rmarkdown::metadata$output) && rmarkdown::metadata$output
## %in% : 'length(x) = 3 > 1' in coercion to 'logical(1)'

library(padr)
library(astsa)

## 
## Attaching package: 'astsa'

## The following object is masked from 'package:forecast':
## 
##     gas

library(tfarima)

## 
## Attaching package: 'tfarima'

## The following object is masked from 'package:TSA':
## 
##     spec

## The following objects are masked from 'package:forecast':
## 
##     easter, seasadj

library(ggplot2)

Import Data

dataset <- read_excel("~/Downloads/DATA KURS FIX MPDW.xlsx")
datagbp <- dataset$GBP

Mengubah Data Menjadi Time Series

datagbp.ts <- ts(datagbp)

Eksplorasi Data

Plot Data Kurs Jual

plot(datagbp.ts,xlab ="Waktu", ylab = "Data Kurs GBP (Ribu)", col="black", main = "Plot Data Deret Waktu Data urs GBP Kurs")
points(datagbp.ts)

Splitting Data (Training dan Testing)

Splitting Data Kurs Jual

datagbp.train <- ts(datagbp[1:145], )
datagbp.test <- ts(datagbp[146:206], start = 146)

Plot Data Kurs Jual

par(mfrow=c(2,1))
plot(datagbp.train, xlab ="Waktu", ylab = "Data Kurs (Ribu)", col="red", main = "Plot Data Training")
points(datagbp.train)

plot(datagbp.test, xlab ="Waktu", ylab = "Data Kurs (Ribu)", col="red", main = "Plot Data Testing")
points(datagbp.test)

Plot Pembagian Training Testing

## Plot Data Training dan testing kurs
ts.plot(datagbp.ts, xlab = "Periode", ylab ="Data Inflasi (Persen)", 
        main = "Plot Data Training dan Data Testing")
lines(datagbp.train, col = "blue")
lines(datagbp.test, col="Red")
legend(0,19600,c("Data Training","Data Testing"), 
       lty=8, col=c("blue","red"), cex=0.8)
abline(v=146, col=c("black"), lty=1, lwd=1)

Smoothing Data

DMA (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(datagbp.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(datagbp.train,  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(datagbp.train,NA), pemulusan_dma=c(pemulusan_dma,NA), ramal_dma)
head(df_dma)

##      df_aktual pemulusan_dma ramal_dma
## [1,]  20177.45            NA        NA
## [2,]  20117.72            NA        NA
## [3,]  20182.52      20152.66        NA
## [4,]  20208.16      20240.56  20152.66
## [5,]  20239.71      20252.53  20240.56
## [6,]  20101.33      20117.11  20252.53

Plot DMA

ts.plot(datagbp.train, xlab="Tahun", ylab="Data Kurs", col="blue", lty=3)
points(datagbp.train)

pemulusan_dma.ts<-ts(pemulusan_dma)
ramal_dma.ts <- ts(ramal_dma)

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"))

Uji Keakuratan dan Error

# 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	9.656670e+05
MSE	6.800472e+03
MAPE	3.326418e-01
RMSE	8.246498e+01
MAD	6.677553e+01

DES (Double Exponential Smoothing)

DES <- HoltWinters(datagbp.train, gamma = FALSE)
tabledes <- data.frame("Alpha Optimum"=c(DES$alpha), "Beta Optimum"=c(DES$beta))
kable(tabledes)

	Alpha.Optimum	Beta.Optimum
alpha	1	0.0556507

ramalanDES <- forecast(DES)
ramalanDES

##     Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 146       19506.81 19410.08 19603.55 19358.87 19654.75
## 147       19495.83 19355.17 19636.49 19280.71 19710.95
## 148       19484.84 19307.81 19661.88 19214.09 19755.59
## 149       19473.86 19263.90 19683.82 19152.75 19794.97
## 150       19462.87 19221.89 19703.85 19094.33 19831.42
## 151       19451.89 19181.04 19722.74 19037.66 19866.12
## 152       19440.90 19140.89 19740.92 18982.07 19899.73
## 153       19429.92 19101.17 19758.67 18927.14 19932.70
## 154       19418.93 19061.69 19776.17 18872.58 19965.28
## 155       19407.95 19022.33 19793.56 18818.20 19997.69

Plot DES

plot(DES)
points(datagbp.train)
legend("topright", c("Data aktual","Pemulusan DES"),lty=1:2,col=c ("black","red"))

Uji Keakuratan dan Error

# ukuran keakuratan
SSEDES<-DES$SSE 

MSEDES<-DES$SSE/length(datagbp.train) 

RMSEDES<-sqrt(MSEDES)

# Mencari MAPE dengan manual
fittedDES<-ramalanDES$fitted 
sisaanDES<-ramalanDES$residuals 
residDES<-datagbp.train-ramalanDES$fitted 
MAPEDES <-sum(abs(sisaanDES[3:length(datagbp.train)]/datagbp.train[3:length(datagbp.train)])*100)/length(datagbp.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	8.143943e+05
MSE	5.616512e+03
RMSE	7.494339e+01
MAPE	2.956961e-01

Perbandingan DMA dan DES

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	8.143943e+05	9.656670e+05
MSE	5.616512e+03	6.800472e+03
RMSE	7.494339e+01	8.246498e+01
MAPE	2.956961e-01	3.326418e-01

Uji Keakuratan dan Error

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

## Warning in c(ramalanDES$mean) - c(datagbp.test): longer object length is not a
## multiple of shorter object length

selisih.test <- na.omit(selisih.test)
SSEtesting<-sum(selisih.test^2) 
MSEtesting<-SSEtesting/length(datagbp.test) 
RMSEtesting<-sqrt(MSEtesting)
MAPEtesting<-sum(abs(selisih.test))/length(datagbp.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	3343118.0341
MSE	54805.2137
RMSE	234.1051
MAPE	204.9760

Pemodelan Arima

Differencing 1

data.dif1<-diff(datagbp.train,differences = 1) 
plot.ts(data.dif1,lty=1,ylab= "Data Inflasi Pembedaan 1", main="Plot Differencing Data Kurs")

Uji Kestasioneran Data

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 = -6.0691, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Plot Data ACF

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

Plot Data PACF

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

Plot Data EACF

eacf(data.dif1)

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

Differencing 2

data.dif2<-diff(datagbp.train,differences = 2) 
plot.ts(data.dif2,lty=1,ylab= "Data Inflasi Pembedaan 1", main="Plot Differencing Data Kurs")

Cek Kestasioneran Data

adf.test(data.dif2)

## Warning in adf.test(data.dif2): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  data.dif2
## Dickey-Fuller = -8.0178, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Plot Data ACF

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

cut off di lag 1

Plot Data PACF

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

cut off di lag off

Plot Data EACF

eacf(data.dif2)

## 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 x 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 x x x o o o o o o o o  o  o  o 
## 4 x o x o o o o o o o o  o  o  o 
## 5 x o x o o o o o o o o  o  o  o 
## 6 x x x o x o x o o o o  o  o  o 
## 7 x x o x x o o o o o o  o  o  o

diperoleh model : 1. Arima(0,2,1) 2. Arima(1,2,2) 3. Arima(2,2,2) 5. Arima(3,2,3) 6. Arima(1,2,1)

Pemodelan

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

Penduga Parameter Model

Arima(0,2,1)

summary(model1)

## Series: data.dif2 
## ARIMA(0,2,1) 
## 
## Coefficients:
##           ma1
##       -1.0000
## s.e.   0.0175
## 
## sigma^2 = 31702:  log likelihood = -932.72
## AIC=1869.43   AICc=1869.52   BIC=1875.33
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set 4.779887 176.1719 138.2666 Inf  Inf 0.9895546 -0.6509523

lmtest::coeftest((model1))

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -1.000000   0.017453 -57.296 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Arima(1,2,2)

summary(model2)

## Series: data.dif2 
## ARIMA(1,2,2) 
## 
## Coefficients:
##           ar1      ma1     ma2
##       -0.4698  -1.9957  0.9999
## s.e.   0.0743   0.0263  0.0256
## 
## sigma^2 = 8462:  log likelihood = -844
## AIC=1695.99   AICc=1696.29   BIC=1707.79
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE      MASE       ACF1
## Training set 3.254189 90.36892 72.42125 NaN  Inf 0.5183085 -0.1423265

lmtest::coeftest((model2))

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.469768   0.074255  -6.3264  2.51e-10 ***
## ma1 -1.995653   0.026329 -75.7981 < 2.2e-16 ***
## ma2  0.999904   0.025594  39.0676 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Arima(2,2,2)

summary(model3)

## Series: data.dif2 
## ARIMA(2,2,2) 
## 
## Coefficients:
##           ar1      ar2      ma1     ma2
##       -0.6191  -0.3044  -1.9967  1.0000
## s.e.   0.0805   0.0801   0.0264  0.0261
## 
## sigma^2 = 7657:  log likelihood = -837.19
## AIC=1684.38   AICc=1684.83   BIC=1699.13
## 
## Training set error measures:
##                   ME     RMSE      MAE  MPE MAPE      MASE        ACF1
## Training set 2.76951 85.65117 68.13987 -Inf  Inf 0.4876674 -0.04692459

lmtest::coeftest((model3))

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.619134   0.080461  -7.6949 1.416e-14 ***
## ar2 -0.304352   0.080059  -3.8016 0.0001438 ***
## ma1 -1.996673   0.026425 -75.5589 < 2.2e-16 ***
## ma2  0.999994   0.026128  38.2728 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Arima(3,2,3)

summary(model4)

## Series: data.dif2 
## ARIMA(3,2,3) 
## 
## Coefficients:
##           ar1      ar2      ar3      ma1     ma2      ma3
##       -0.0537  -0.0942  -0.0828  -2.8189  2.6413  -0.8220
## s.e.   0.0898   0.0875   0.0864   0.0378  0.0760   0.0388
## 
## sigma^2 = 6248:  log likelihood = -824.69
## AIC=1663.38   AICc=1664.23   BIC=1684.03
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE      MASE        ACF1
## Training set 2.682303 76.80059 60.74851 Inf  Inf 0.4347684 -0.01194438

lmtest::coeftest((model4))

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value Pr(>|z|)    
## ar1 -0.053742   0.089773  -0.5986   0.5494    
## ar2 -0.094194   0.087494  -1.0766   0.2817    
## ar3 -0.082791   0.086399  -0.9582   0.3379    
## ma1 -2.818851   0.037815 -74.5429   <2e-16 ***
## ma2  2.641340   0.075974  34.7665   <2e-16 ***
## ma3 -0.822023   0.038834 -21.1676   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Arima(1,2,1)

summary(model5)

## Series: data.dif2 
## ARIMA(1,2,1) 
## 
## Coefficients:
##           ar1      ma1
##       -0.6512  -1.0000
## s.e.   0.0633   0.0176
## 
## sigma^2 = 18183:  log likelihood = -893.8
## AIC=1793.6   AICc=1793.77   BIC=1802.44
## 
## Training set error measures:
##                    ME     RMSE      MAE MPE MAPE     MASE       ACF1
## Training set 4.514104 132.9469 106.5713 NaN  Inf 0.762716 -0.3346775

lmtest::coeftest((model5))

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.651248   0.063301 -10.288 < 2.2e-16 ***
## ma1 -0.999996   0.017622 -56.748 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Perbandingan Kebaikan Model

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

modelaccuracy

##          Model      AIC      BIC     Signifikansi
## 1 ARIMA(0,2,1) 1869.434 1875.332       Signifikan
## 2 ARIMA(1,2,2) 1695.993 1707.788       Signifikan
## 3 ARIMA(2,2,2) 1684.384 1699.128       Signifikan
## 4 ARIMA(3,2,3) 1663.385 1684.026 Tidak Signifikan
## 5 ARIMA(1,2,1) 1793.597 1802.443       Signifikan

Dengan model terbaik Arima(2,2,2)

Diagnostik Model

Overfitting

menggunakan model Arima(2,2,2) dengan Arima(3,2,2) dan Arima(2,2,3)

model5a <- Arima(data.dif2, order=c(3,2,2), method="ML")
model5b <- Arima(data.dif2, order=c(2,2,3), method="ML")

Model Arima(2,2,2)

coeftest(model5)

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.651248   0.063301 -10.288 < 2.2e-16 ***
## ma1 -0.999996   0.017622 -56.748 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model Arima(3,2,2)

coeftest(model5a)

## 
## z test of coefficients:
## 
##      Estimate Std. Error  z value  Pr(>|z|)    
## ar1 -0.675460   0.083681  -8.0718 6.926e-16 ***
## ar2 -0.415024   0.095306  -4.3546 1.333e-05 ***
## ar3 -0.171327   0.083550  -2.0506   0.04031 *  
## ma1 -1.996895   0.026545 -75.2279 < 2.2e-16 ***
## ma2  0.999940   0.026324  37.9857 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model Arima(2,2,3)

coeftest(model5b)

## Warning in sqrt(diag(se)): NaNs produced

## 
## z test of coefficients:
## 
##      Estimate Std. Error z value Pr(>|z|)
## ar1 -0.057046   0.086750 -0.6576   0.5108
## ar2 -0.071034   0.086805 -0.8183   0.4132
## ma1 -2.876170        NaN     NaN      NaN
## ma2  2.760864        NaN     NaN      NaN
## ma3 -0.884392        NaN     NaN      NaN

Uji Diagnostik

Diagnostik Model : Eksploratif

sisaan5 <- model5a$residuals

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

Diagnostik Model : Uji Formal

Sisaan Menyebar Normal

shapiro.test(sisaan5)

## 
##  Shapiro-Wilk normality test
## 
## data:  sisaan5
## W = 0.99351, p-value = 0.7664

Sisaan Saling Bebas

Box.test(sisaan5, type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sisaan5
## X-squared = 0.20539, df = 1, p-value = 0.6504

Nilai Tengah Sisaan Sama dengan Nol

t.test(sisaan5, mu = 0, conf.level = 0.95)

## 
##  One Sample t-test
## 
## data:  sisaan5
## t = 0.27576, df = 142, p-value = 0.7831
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -12.01121  15.90555
## sample estimates:
## mean of x 
##  1.947169

Forcasting

ramalan <- forecast(Arima(datagbp.train, order=c(3,2,2),method="ML",include.drift = TRUE),h = length(datagbp.test))

## Warning in Arima(datagbp.train, order = c(3, 2, 2), method = "ML", include.drift
## = TRUE): No drift term fitted as the order of difference is 2 or more.

data.ramalan <- ramalan$mean
data.ramalan.ts <- ts(data.ramalan, start = 146)
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(datagbp.train, ylab = "Data Inflasi(Persen)", col="black",lwd=2,main="Forecasting ARIMA(3,2,2)",gpars = list(col.main="black",col.axis="black",col.sub="black"), xlim = c(0,211), ylim = c(19000,21000))
lines(data.ramalan.ts, col = "blue",lwd=2)
lines(datagbp.test, col = "red", lwd =2)
legend(61,16400,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(datagbp.test[1:10], data.ramalan[1:10]), nrow = 10, ncol = 2)
colnames(perbandingan.temp)<-c("Aktual","Hasil Forecast")
head(perbandingan.temp)

##        Aktual Hasil Forecast
## [1,] 19409.94       19514.97
## [2,] 19439.03       19502.44
## [3,] 19445.60       19484.68
## [4,] 19490.78       19466.04
## [5,] 19428.94       19447.40
## [6,] 19478.72       19429.07

error <- data.frame(datagbp.test)-data.frame(data.ramalan[1:61]) 

## 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(datagbp.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 7.928205e+06
## 2               MSE 1.299706e+05
## 3              MAPE 1.635303e+00
## 4              RMSE 3.605143e+02
## 5               MAD 3.153171e+02

kable(akurasifarima)

Ukuran.Keakuratan	Forecasting
SSE	7.928205e+06
MSE	1.299706e+05
MAPE	1.635303e+00
RMSE	3.605143e+02
MAD	3.153171e+02

ARIMAX

Data

datausd <- dataset$USD
dataemas <- dataset$EMAS
datasaham <- dataset$TESLA
datausd.ts <- ts(datausd)
dataemas.ts <- ts(dataemas)
datasaham.ts <- ts(datasaham)

Plot Time Series

plot(datausd.ts,xlab ="Waktu", ylab = "Data Inflasi Global(Persen)", col="black", main = "Plot Data Inflasi Global (Persen)")
points(datasaham.ts)

plot(dataemas.ts,xlab ="Waktu", ylab = "Data Inflasi Global(Persen)", col="black", main = "Plot Data Inflasi Global (Persen)")
points(dataemas.ts)

plot(datasaham.ts,xlab ="Waktu", ylab = "Data Inflasi Global(Persen)", col="black", main = "Plot Data Inflasi Global (Persen)")
points(datasaham.ts)

Splitting Data

datausd.train <- ts(datausd[1:145], )
datausd.test <- ts(datausd[146:206], start = 146)

dataemas.train <- ts(dataemas[1:145], )
dataemas.test <- ts(dataemas[146:206], start = 146)

datasaham.train <- ts(datasaham[1:145], )
datasaham.test <- ts(datasaham[146:206], start = 146)

Plot Data

par(mfrow=c(2,1))
plot(datausd.train,xlab ="Waktu", ylab = "Data Inflasi (Persen)", col="red", main = "Data Training Inflasi")
points(datausd.train)

plot(datausd.test,xlab ="Waktu", ylab = "Data Inflasi (Persen)", col="red", main = "Data Testing Inflasi")
points(datausd.test)

par(mfrow=c(2,1))
plot(dataemas.train,xlab ="Waktu", ylab = "Data Inflasi (Persen)", col="red", main = "Data Training Inflasi")
points(dataemas.train)

plot(dataemas.test,xlab ="Waktu", ylab = "Data Inflasi (Persen)", col="red", main = "Data Testing Inflasi")
points(dataemas.test)

par(mfrow=c(2,1))
plot(datasaham.train,xlab ="Waktu", ylab = "Data Inflasi (Persen)", col="red", main = "Data Training Inflasi")
points(datasaham.train)

plot(datasaham.test,xlab ="Waktu", ylab = "Data Inflasi (Persen)", col="red", main = "Data Testing Inflasi")
points(datasaham.test)

### Korelasi Data

korelasix <- c("USD vs EMAS" = cor(datausd, dataemas),"USD vs SAHAM" =cor(datausd, datasaham),"EMAS vs SAHAM" =cor(dataemas, datasaham))
korx <- data.frame("Nilai Korelasi" = korelasix)
kable(korx)

	Nilai.Korelasi
USD vs EMAS	-0.3209930
USD vs SAHAM	-0.4054646
EMAS vs SAHAM	-0.0504935

antar variabel x gaada yang besar nilai kor (tidak ada yg multikol)

korelasitest <- c("GBP vs USD" = cor(datagbp.test, dataemas.test), "GBP vs SAHAM" =cor(datagbp.test, datasaham.test),"GBP vs EMAS" =cor(datagbp.test, dataemas.test))
kort <- data.frame("Nilai Korelasi" = korelasitest)
kable(kort)

	Nilai.Korelasi
GBP vs USD	-0.0743528
GBP vs SAHAM	-0.0349986
GBP vs EMAS	-0.0743528

korelasitest <- c("GBP vs USD" = cor(datagbp, dataemas), "GBP vs SAHAM" =cor(datagbp, datasaham),"GBP vs EMAS" =cor(datagbp, dataemas))
kort <- data.frame("Nilai Korelasi" = korelasitest)
kable(kort)

	Nilai.Korelasi
GBP vs USD	0.1418983
GBP vs SAHAM	-0.8383975
GBP vs EMAS	0.1418983

korelasitest <- c("USD vs EMAS" = cor(datausd.test, dataemas.test), "USD vs SAHAM" =cor(datausd.test, datasaham.test),"SAHAM vs EMAS" =cor(datasaham.test, dataemas.test))
kort <- data.frame("Nilai Korelasi" = korelasitest)
kable(kort)

	Nilai.Korelasi
USD vs EMAS	-0.0636976
USD vs SAHAM	0.2130446
SAHAM vs EMAS	0.3652165

Pembentukan Model Regresi

X1 = datausd
X2 = dataemas
X3 = datasaham
y <- datagbp
reg2  <- lm(y~X1+X2+X3)
reg2

## 
## Call:
## lm(formula = y ~ X1 + X2 + X3)
## 
## Coefficients:
## (Intercept)           X1           X2           X3  
##   4126.3609       0.9533       1.8235      -5.0303

summary(reg2)

## 
## Call:
## lm(formula = y ~ X1 + X2 + X3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -505.82 -113.63    8.76  111.76  663.81 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 4126.3609  2216.6529   1.862   0.0641 .  
## X1             0.9533     0.1310   7.280 7.30e-12 ***
## X2             1.8235     0.3329   5.478 1.27e-07 ***
## X3            -5.0303     0.2654 -18.956  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 199.9 on 202 degrees of freedom
## Multiple R-squared:  0.7725, Adjusted R-squared:  0.7691 
## F-statistic: 228.7 on 3 and 202 DF,  p-value: < 2.2e-16

Metode Best Subset Regression

library(olsrr)

## 
## Attaching package: 'olsrr'

## The following object is masked from 'package:MASS':
## 
##     cement

## The following object is masked from 'package:datasets':
## 
##     rivers

bs.hc <- ols_step_best_subset(reg2)
bs.hc

## Best Subsets Regression
## -----------------------
## Model Index    Predictors
## -----------------------
##      1         X3       
##      2         X1 X3    
##      3         X1 X2 X3 
## -----------------------
## 
##                                                            Subsets Regression Summary                                                           
## ------------------------------------------------------------------------------------------------------------------------------------------------
##                        Adj.        Pred                                                                                                          
## Model    R-Square    R-Square    R-Square     C(p)         AIC         SBIC          SBC           MSEP            FPE          HSP        APC  
## ------------------------------------------------------------------------------------------------------------------------------------------------
##   1        0.7029      0.7015      0.6954    61.8180    2824.2891    2238.7137    2834.2727    10646807.9911    52185.2916    254.5986    0.3029 
##   2        0.7387      0.7362      0.7297    32.0075    2799.8202    2214.5198    2813.1317     9409394.2346    46340.7192    226.1165    0.2690 
##   3        0.7725      0.7691      0.7631     4.0000    2773.2888    2188.8437    2789.9282     8233156.7512    40740.8283    198.8299    0.2365 
## ------------------------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria 
##  SBIC: Sawa's Bayesian Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
##  MSEP: Estimated error of prediction, assuming multivariate normality 
##  FPE: Final Prediction Error 
##  HSP: Hocking's Sp 
##  APC: Amemiya Prediction Criteria

model 3 terbaik dengan nilai aic terendah

sisaanx <- reg2$residuals
par(mfrow=c(2,2))
qqnorm(sisaanx)
box(col="black",lwd=2)
qqline(sisaanx, col = "red", lwd =1, col.lab="black",
       col.axis="black",col.sub = "black")
box(col="black",lwd=2)
plot(c(1:length(sisaanx)),sisaanx,col="black",col.lab="black",col.axis="black")
box(col="black",lwd=2)
acf(sisaanx,col="black",col.sub = "black",col.axis="black", col.lab="black")
box(col="black",lwd=2)
pacf(sisaanx,col="black",col.sub = "black",col.axis="black", col.lab="black",col.main="black")
box(col="black",lwd=2)

## White Noise Regressi

Box.test(sisaanx, lag = 12)

## 
##  Box-Pierce test
## 
## data:  sisaanx
## X-squared = 537.7, df = 12, p-value < 2.2e-16

Box.test(sisaanx, lag = 24)

## 
##  Box-Pierce test
## 
## data:  sisaanx
## X-squared = 553.84, df = 24, p-value < 2.2e-16

Box.test(sisaanx, lag = 36)

## 
##  Box-Pierce test
## 
## data:  sisaanx
## X-squared = 599.66, df = 36, p-value < 2.2e-16

Box.test(sisaanx, lag = 48)

## 
##  Box-Pierce test
## 
## data:  sisaanx
## X-squared = 639.27, df = 48, p-value < 2.2e-16

Pemodelan sisaan regresi

ACF dan PACF

acf(sisaanx)

pacf(sisaanx)

eacf(sisaanx)

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

UJI FORMAL ADF TEST

adf.test(sisaanx)

## Warning in adf.test(sisaanx): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  sisaanx
## Dickey-Fuller = -4.3893, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Dari Matriks EACF dapat diduga model yang cocok adalah model :

ARIMAX(2,0,5)

ARIMAX(3,0,5)

ARIMAX(2,0,6)

ARIMAX(4,0,5)

ARIMAX(1,0,1)

modelx1 <- Arima(datagbp.train, order = c(2,0,5), xreg = cbind(datausd.train,dataemas.train,datasaham.train), method = "ML")
modelx2 <- Arima(datagbp.train, order = c(3,0,5), xreg = cbind(datausd.train,dataemas.train,datasaham.train), method = "ML")
modelx3 <- Arima(datagbp.train, order = c(2,0,6), xreg = cbind(datausd.train,dataemas.train,datasaham.train), method = "ML")
modelx4 <- Arima(datagbp.train, order = c(4,0,5), xreg = cbind(datausd.train,dataemas.train,datasaham.train), method = "ML")
modelx5 <- Arima(datagbp.train, order = c(1,0,1), xreg = cbind(datausd.train,dataemas.train,datasaham.train), method = "ML")

Arimax(2,0,5)

summary(modelx1)

## Series: datagbp.train 
## Regression with ARIMA(2,0,5) errors 
## 
## Coefficients:
##          ar1     ar2     ma1      ma2      ma3      ma4      ma5  intercept
##       0.1808  0.8025  0.7922  -0.1076  -0.1588  -0.2061  -0.1254  15932.268
## s.e.  0.3640  0.3604  0.3700   0.1159   0.1108   0.1162   0.0942   2500.195
##       datausd.train  dataemas.train  datasaham.train
##              0.1876          0.3676           2.6951
## s.e.         0.1654          0.4004           1.0596
## 
## sigma^2 = 5069:  log likelihood = -820.19
## AIC=1664.37   AICc=1666.74   BIC=1700.09
## 
## Training set error measures:
##                     ME     RMSE      MAE        MPE      MAPE      MASE
## Training set -4.847586 68.44612 52.47789 -0.0255555 0.2614566 0.9218812
##                     ACF1
## Training set -0.01325761

coeftest(modelx1)

## 
## z test of coefficients:
## 
##                    Estimate  Std. Error z value  Pr(>|z|)    
## ar1              1.8084e-01  3.6400e-01  0.4968   0.61932    
## ar2              8.0247e-01  3.6038e-01  2.2268   0.02596 *  
## ma1              7.9224e-01  3.6995e-01  2.1415   0.03224 *  
## ma2             -1.0758e-01  1.1591e-01 -0.9282   0.35331    
## ma3             -1.5880e-01  1.1076e-01 -1.4337   0.15165    
## ma4             -2.0605e-01  1.1618e-01 -1.7736   0.07613 .  
## ma5             -1.2537e-01  9.4237e-02 -1.3303   0.18341    
## intercept        1.5932e+04  2.5002e+03  6.3724 1.861e-10 ***
## datausd.train    1.8763e-01  1.6537e-01  1.1346   0.25655    
## dataemas.train   3.6761e-01  4.0038e-01  0.9181   0.35854    
## datasaham.train  2.6951e+00  1.0596e+00  2.5434   0.01098 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Arimax(3,0,5)

summary(modelx2)

## Series: datagbp.train 
## Regression with ARIMA(3,0,5) errors 
## 
## Coefficients:
##           ar1     ar2     ar3     ma1     ma2      ma3      ma4      ma5
##       -0.2920  0.8057  0.4601  1.2607  0.3346  -0.2267  -0.3120  -0.1985
## s.e.   0.5334  0.1913  0.4031  0.5248  0.3964   0.1594   0.1477   0.0952
##       intercept  datausd.train  dataemas.train  datasaham.train
##       15489.033         0.2096          0.4370           2.7051
## s.e.   2499.783         0.1651          0.4035           1.0424
## 
## sigma^2 = 5065:  log likelihood = -819.63
## AIC=1665.26   AICc=1668.04   BIC=1703.96
## 
## Training set error measures:
##                   ME     RMSE      MAE         MPE      MAPE     MASE
## Training set -4.7773 68.15924 52.36688 -0.02519341 0.2609256 0.919931
##                      ACF1
## Training set -0.006930038

coeftest(modelx2)

## 
## z test of coefficients:
## 
##                    Estimate  Std. Error z value  Pr(>|z|)    
## ar1                -0.29200     0.53336 -0.5475  0.584057    
## ar2                 0.80567     0.19129  4.2118 2.533e-05 ***
## ar3                 0.46014     0.40307  1.1416  0.253624    
## ma1                 1.26072     0.52480  2.4023  0.016294 *  
## ma2                 0.33463     0.39642  0.8441  0.398607    
## ma3                -0.22666     0.15940 -1.4220  0.155024    
## ma4                -0.31199     0.14774 -2.1118  0.034701 *  
## ma5                -0.19851     0.09523 -2.0845  0.037112 *  
## intercept       15489.03298  2499.78281  6.1962 5.786e-10 ***
## datausd.train       0.20960     0.16510  1.2695  0.204257    
## dataemas.train      0.43703     0.40346  1.0832  0.278713    
## datasaham.train     2.70508     1.04235  2.5952  0.009455 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Arimax(2,0,6)

summary(modelx3)

## Series: datagbp.train 
## Regression with ARIMA(2,0,6) errors 
## 
## Coefficients:

## Warning in sqrt(diag(x$var.coef)): NaNs produced

##          ar1     ar2     ma1      ma2      ma3      ma4     ma5     ma6
##       0.8591  0.1292  0.1032  -0.1077  -0.1175  -0.1677  0.0527  0.0482
## s.e.     NaN     NaN     NaN   0.0828      NaN      NaN     NaN     NaN
##       intercept  datausd.train  dataemas.train  datasaham.train
##       15485.726         0.2037          0.4714           2.8436
## s.e.   2523.265         0.1647          0.4091           1.0744
## 
## sigma^2 = 5081:  log likelihood = -819.85
## AIC=1665.69   AICc=1668.47   BIC=1704.39
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE      MAPE      MASE
## Training set -4.591221 68.26456 52.26469 -0.02425419 0.2604019 0.9181358
##                      ACF1
## Training set -0.004918638

coeftest(modelx3)

## Warning in sqrt(diag(se)): NaNs produced

## 
## z test of coefficients:
## 
##                    Estimate  Std. Error z value Pr(>|z|)    
## ar1              8.5907e-01         NaN     NaN      NaN    
## ar2              1.2920e-01         NaN     NaN      NaN    
## ma1              1.0324e-01         NaN     NaN      NaN    
## ma2             -1.0766e-01  8.2821e-02 -1.2999  0.19364    
## ma3             -1.1747e-01         NaN     NaN      NaN    
## ma4             -1.6770e-01         NaN     NaN      NaN    
## ma5              5.2653e-02         NaN     NaN      NaN    
## ma6              4.8221e-02         NaN     NaN      NaN    
## intercept        1.5486e+04  2.5233e+03  6.1372  8.4e-10 ***
## datausd.train    2.0374e-01  1.6471e-01  1.2370  0.21610    
## dataemas.train   4.7136e-01  4.0911e-01  1.1522  0.24925    
## datasaham.train  2.8436e+00  1.0744e+00  2.6466  0.00813 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Arimax(4,0,5)

summary(modelx4)

## Series: datagbp.train 
## Regression with ARIMA(4,0,5) errors 
## 
## Coefficients:
##          ar1      ar2     ar3     ar4     ma1     ma2     ma3      ma4      ma5
##       0.3958  -0.3279  0.4552  0.4369  0.5711  0.8342  0.3026  -0.2426  -0.0746
## s.e.  0.9903   0.8916  0.7796  0.8515  0.9926  0.1581  0.8157   0.1811   0.1999
##       intercept  datausd.train  dataemas.train  datasaham.train
##       15167.802         0.2330          0.4407           2.6940
## s.e.   2554.521         0.1641          0.4064           0.9934
## 
## sigma^2 = 4876:  log likelihood = -818.2
## AIC=1664.4   AICc=1667.63   BIC=1706.07
## 
## Training set error measures:
##                    ME     RMSE      MAE         MPE     MAPE     MASE
## Training set -3.67632 66.62151 52.13565 -0.01960723 0.259782 0.915869
##                      ACF1
## Training set 0.0009493526

coeftest(modelx4)

## 
## z test of coefficients:
## 
##                    Estimate  Std. Error z value  Pr(>|z|)    
## ar1              3.9584e-01  9.9031e-01  0.3997  0.689369    
## ar2             -3.2788e-01  8.9159e-01 -0.3678  0.713057    
## ar3              4.5523e-01  7.7955e-01  0.5840  0.559241    
## ar4              4.3688e-01  8.5148e-01  0.5131  0.607895    
## ma1              5.7112e-01  9.9263e-01  0.5754  0.565044    
## ma2              8.3424e-01  1.5811e-01  5.2764 1.318e-07 ***
## ma3              3.0263e-01  8.1572e-01  0.3710  0.710642    
## ma4             -2.4258e-01  1.8111e-01 -1.3394  0.180428    
## ma5             -7.4572e-02  1.9992e-01 -0.3730  0.709141    
## intercept        1.5168e+04  2.5545e+03  5.9376 2.892e-09 ***
## datausd.train    2.3304e-01  1.6406e-01  1.4205  0.155473    
## dataemas.train   4.4071e-01  4.0643e-01  1.0843  0.278211    
## datasaham.train  2.6940e+00  9.9344e-01  2.7118  0.006692 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Arimax(1,0,1)

summary(modelx5)

## Series: datagbp.train 
## Regression with ARIMA(1,0,1) errors 
## 
## Coefficients:
##          ar1     ma1  intercept  datausd.train  dataemas.train  datasaham.train
##       0.9743  0.0050  16793.903         0.1634          0.1094           2.7737
## s.e.  0.0201  0.0922   2519.156         0.1706          0.3768           1.0990
## 
## sigma^2 = 5048:  log likelihood = -822.37
## AIC=1658.73   AICc=1659.55   BIC=1679.57
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE      MAPE      MASE
## Training set -2.171782 69.56307 54.73204 -0.01208502 0.2727293 0.9614798
##                     ACF1
## Training set 0.002006945

coeftest(modelx5)

## 
## z test of coefficients:
## 
##                   Estimate Std. Error z value  Pr(>|z|)    
## ar1             9.7429e-01 2.0075e-02 48.5330 < 2.2e-16 ***
## ma1             4.9793e-03 9.2163e-02  0.0540   0.95691    
## intercept       1.6794e+04 2.5192e+03  6.6665  2.62e-11 ***
## datausd.train   1.6341e-01 1.7061e-01  0.9578   0.33817    
## dataemas.train  1.0941e-01 3.7675e-01  0.2904   0.77150    
## datasaham.train 2.7737e+00 1.0990e+00  2.5239   0.01161 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

modelxaccuracy<-data.frame(
  "Model"=c("ARIMAX(2,0,5)", "ARIMAX(2,0,6)", "ARIMAX(3,0,5)", "ARIMAX(4,0,5)","ARIMAX(1,0,1)"),
  "AIC"=c(modelx1$aic, modelx2$aic, modelx3$aic,modelx4$aic,modelx5$aic),
  "BIC"=c(modelx1$bic, modelx2$bic, modelx3$bic,modelx4$bic,modelx5$bic),
  "Signifikansi"=c("Signifikan","Signifikan","Tidak Signifikan","Tidak Signifikan"," Signifikan"))

modelxaccuracy

##           Model      AIC      BIC     Signifikansi
## 1 ARIMAX(2,0,5) 1664.373 1700.094       Signifikan
## 2 ARIMAX(2,0,6) 1665.261 1703.959       Signifikan
## 3 ARIMAX(3,0,5) 1665.693 1704.390 Tidak Signifikan
## 4 ARIMAX(4,0,5) 1664.397 1706.072 Tidak Signifikan
## 5 ARIMAX(1,0,1) 1658.730 1679.567       Signifikan

Overfitting

model5x.a <- Arima(datagbp.train, order = c(1,0,2), xreg = cbind(datausd.train,dataemas.train,datasaham.train), method = "ML")
coeftest(model5x.a)

## 
## z test of coefficients:
## 
##                    Estimate  Std. Error z value  Pr(>|z|)    
## ar1              9.8135e-01  1.7965e-02 54.6251 < 2.2e-16 ***
## ma1             -1.7204e-02  9.0836e-02 -0.1894    0.8498    
## ma2             -1.0182e-01  1.1451e-01 -0.8892    0.3739    
## intercept        1.6138e+04  2.6053e+03  6.1944 5.851e-10 ***
## datausd.train    1.9170e-01  1.7075e-01  1.1227    0.2616    
## dataemas.train   2.3473e-01  4.0857e-01  0.5745    0.5656    
## datasaham.train  2.7795e+00  1.0848e+00  2.5623    0.0104 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model5x.a)

## Series: datagbp.train 
## Regression with ARIMA(1,0,2) errors 
## 
## Coefficients:
##          ar1      ma1      ma2  intercept  datausd.train  dataemas.train
##       0.9813  -0.0172  -0.1018  16138.260         0.1917          0.2347
## s.e.  0.0180   0.0908   0.1145   2605.304         0.1708          0.4086
##       datasaham.train
##                2.7795
## s.e.           1.0848
## 
## sigma^2 = 5054:  log likelihood = -821.97
## AIC=1659.94   AICc=1660.99   BIC=1683.75
## 
## Training set error measures:
##                    ME     RMSE      MAE         MPE      MAPE      MASE
## Training set -2.89548 69.35253 54.19989 -0.01573492 0.2700963 0.9521316
##                     ACF1
## Training set 0.007025035

model5x.b <- Arima(datagbp.train, order = c(2,0,1), xreg = cbind(datausd.train,dataemas.train,datasaham.train), method = "ML")
coeftest(model5x.b)

## 
## z test of coefficients:
## 
##                   Estimate Std. Error z value Pr(>|z|)    
## ar1             2.5892e-01 6.0307e-01  0.4293  0.66768    
## ar2             6.9531e-01 5.8968e-01  1.1791  0.23835    
## ma1             7.4329e-01 5.7527e-01  1.2921  0.19633    
## intercept       1.6574e+04 2.5681e+03  6.4539 1.09e-10 ***
## datausd.train   1.7478e-01 1.7253e-01  1.0131  0.31102    
## dataemas.train  1.4517e-01 3.8715e-01  0.3750  0.70768    
## datasaham.train 2.7367e+00 1.0980e+00  2.4925  0.01269 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(model5x.b)

## Series: datagbp.train 
## Regression with ARIMA(2,0,1) errors 
## 
## Coefficients:
##          ar1     ar2     ma1  intercept  datausd.train  dataemas.train
##       0.2589  0.6953  0.7433  16574.399         0.1748          0.1452
## s.e.  0.6031  0.5897  0.5753   2568.113         0.1725          0.3871
##       datasaham.train
##                2.7367
## s.e.           1.0980
## 
## sigma^2 = 5077:  log likelihood = -822.25
## AIC=1660.5   AICc=1661.56   BIC=1684.32
## 
## Training set error measures:
##                     ME     RMSE      MAE         MPE      MAPE      MASE
## Training set -2.063197 69.51066 54.49082 -0.01153924 0.2715181 0.9572424
##                     ACF1
## Training set -0.01889429

Diagnostik Model Arimax(1,0,1)

sisaan.arimax <- modelx5$residuals

Box.test(sisaan.arimax)

## 
##  Box-Pierce test
## 
## data:  sisaan.arimax
## X-squared = 0.00058403, df = 1, p-value = 0.9807

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

acf(sisaan.arimax, main = "Plot ACF Sisaan Model ARIMAX")

pacf(sisaan.arimax, lag.max = 24, main = "Plot ACF Sisaan Model ARIMAX")

Forecast

ramalanx <- forecast((modelx5), xreg = cbind(datausd.test,dataemas.test,datasaham.test))

## Warning in forecast.forecast_ARIMA((modelx5), xreg = cbind(datausd.test, : xreg
## contains different column names from the xreg used in training. Please check
## that the regressors are in the same order.

data.ramalanx <- ramalanx$mean
data.ramalan.tsx <- ts(data.ramalanx, start = 146)
plot(ramalanx,col="black",col.sub ="black",col.axis="black",
     col.lab="black",col.main="black",lwd=2)
box(col="black",lwd=2)

perbandinganx <- matrix(c(datagbp.test[1:10], data.ramalanx[1:10]), nrow = 10, ncol = 2, byrow = FALSE)
colnames(perbandinganx) <- c("Aktual", "Ramalan")
perbandinganx

##         Aktual  Ramalan
##  [1,] 19409.94 19539.05
##  [2,] 19439.03 19543.57
##  [3,] 19445.60 19560.25
##  [4,] 19490.78 19583.62
##  [5,] 19428.94 19605.98
##  [6,] 19478.72 19625.74
##  [7,] 19509.23 19646.26
##  [8,] 19431.19 19670.81
##  [9,] 19440.66 19679.05
## [10,] 19519.06 19716.02

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

error <- data.frame(datagbp.test)-data.frame(data.ramalanx[1:61]) 

## 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(datagbp.test))*100 # Relative Error
MAPE <- sapply(abs(r.error), mean, na.rm = T)

akurasiarimax <- data.frame(
  "Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"), 
  "Forecasting" = c(SSE, MSE, MAPE, RMSE, MAD))
kable(akurasiarimax)

Ukuran.Keakuratan	Forecasting
SSE	3.731259e+07
MSE	6.116817e+05
MAPE	3.670465e+00
RMSE	7.821008e+02
MAD	7.059936e+02

ARIMA VS ARIMAX

forecasterror <- data.frame(akurasifarima, akurasiarimax)
kable(forecasterror)

Ukuran.Keakuratan	Forecasting	Ukuran.Keakuratan.1	Forecasting.1
SSE	7.928205e+06	SSE	3.731259e+07
MSE	1.299706e+05	MSE	6.116817e+05
MAPE	1.635303e+00	MAPE	3.670465e+00
RMSE	3.605143e+02	RMSE	7.821008e+02
MAD	3.153171e+02	MAD	7.059936e+02

ts.plot(datagbp.ts,xlab = "Periode", ylab = "Data Inflasi(Persen)", col="black",lwd=2,main="Forecasting ARIMAX(2,1,0) vs Forecasting ARIMA(2,1,0)",gpars = list(col.main="black",col.axis="black",col.sub="black"))
lines(data.ramalan.tsx, col = "blue",lwd=2)
lines(data.ramalan.ts, col = "purple", lwd=2)
lines(datagbp.test, col = "red", lwd =2)
abline(v=2020, col=c("green"), lty=1, lwd=1)
legend(2019,9,c("Data Testing","Data Forecast ARIMAX", "Data Forecast ARIMA", "Split Data"), 
       lwd=2, col=c("red","blue","purple","green"), cex=0.8)
box(col="black",lwd=2)

korarimax <- cor(datagbp.test,data.ramalan.tsx)
korarima <- cor(datagbp.test,data.ramalan.ts)
framekorelasi <- data.frame(
  "Variabel Korelasi" = c("ARIMAX(1,0,1) vs Data Aktual", "ARIMA(3,2,2) vs Data Aktual"), 
  "Nilai Korelasi" = c(korarimax, korarima))
kable(framekorelasi)

Variabel.Korelasi	Nilai.Korelasi
ARIMAX(1,0,1) vs Data Aktual	-0.5291303
ARIMA(3,2,2) vs Data Aktual	0.7975509