Analisis Data Deret Waktu Nilai Impor Bulanan Periode Januari 2005 - September 2022

Pendahuluan

Impor adalah pembelian suatu barang dari luar negeri untuk dipergunakan atau dipasarkan kembali di dalam negeri. Menurut Benny, (2013), impor merupakan transportasi barang dari suatu negara ke negara lain secara legal, umumnya dalam proses perdagangan. Alih-alih menambah pendapatan negara, aktivitas impor merupakan salah satu pengeluaran negara untuk memenuhi sumber daya yang kurang di negara tersebut. Tingginya impor menyebabkan produktivitas suatu negara menurun. Hal ini mengakibatkan meningkatnya angka pengangguran dan menurunnya pendapatan sehingga daya beli masyarakat menurun (Sedyaningrum & Nuzula, 2016).

Aktivitas impor merupakan aktivitas yang sangat krusial dalam aktivitas perdagangan internasional, terutama untuk Indonesia. Meskipun Indonesia merupakan salah satu negara dengan kekayaan sumber daya alam, tidak dapat dipungkiri Indonesia masih kekurangan bahan baku tertentu yang perlu diimpor dari negara lain. Bahkan, ada beberapa produksi barang yang hanya bisa dilakukan di luar negeri padahal bahan baku produksi tersebut berasal dari Indonesia. Hal ini karena Indonesia belum memiliki teknologi yang layak untuk produksi barang tersebut sehingga Indonesia harus mengimpor barang yang telah jadi dari bahan baku tersebut (Farina & Husaini, 2015).

Kondisi nilai impor di Indonesia periode 2010-2019 mengalami fluktuasi karena terjadi peristiwa-peristiwa penting yang mengubah wajah perekonomian Indonesia diantaranya pemulihan pasca krisis ekonomi global 2008, Pilpres Periode 1 (2014), dan Pilpres periode 2 (2019), disusul oleh kasus COVID-19 yang muncul perlahan di akhir tahun 2019. Kenaikan impor tertinggi di Indonesia terjadi pada tahun 2012, tepatnya saat kegiatan impor migas dan non-migas meningkat. Diantara faktor yang mempengaruhi impor di Indonesia antara lain kurs dollar, inflasi, dan cadangan devisa (Sukirno, 2013). Oleh karena itu, perlu metode yang tepat untuk melakukan peramalan data nilai impor yang berfluktuasi agar hasil peramalan lebih akurat.

Peramalan merupakan teknik yang dapat digunakan dalam membuat suatu prediksi di masa depan berdasarkan data periode sebelumnya. Hal yang perlu diperhatikan dalam peramalan adalah pola data sebelumnya agar dapat diputuskan metode peramalan yang sesuai untuk memprediksi data yang ada. Salah satu metode yang dapat digunakan untuk melakukan peramalan adalah Autoregressive Integrated Moving Average (ARIMA) (Nasirudin et al., 2022).

Metode ARIMA merupakan salah satu teknik analisis data deret waktu yang banyak digunakan untuk peramalan data masa depan. ARIMA menggunakan nilai masa lalu dan sekarang untuk menghasilkan peramalan jangka pendek yang akurat (Bangun, 2016). Beberapa peneliti menggunakan model ARIMA untuk meramalkan kondisi masa depan seperti yang dilakukan oleh Purnama, (2020) yang melakukan peramalan nilai ekspor dan impor Indonesia menggunakan metode ARIMA Box-Jenkins.

Seasonal Autoregressive Integrated Moving Average (SARIMA) merupakan pengembangan dari model ARIMA yang memiliki unsur musiman. Beberapa peneliti menggunakan model SARIMA untuk meramalkan kondisi masa depan seperti yang dilakukan oleh Aniendyah, (2020) yang melakukan peramalan hasil ekspor dan impor menggunakan metode Seasonal Autoregressive Integrated Moving Average (SARIMA).

Pada penelitian ini, penulis ingin menerapkan model Seasonal Autoregressive Integrated Moving Average (SARIMA) pada data nilai impor bulanan karena model ini sangat sesuai digunakan pada data yang memiliki unsur musiman seperti nilai impor. Adapun tujuan dari penelitian ini adalah mendapatkan model terbaik untuk memprediksi nilai impor di masa mendatang. Hal ini diharapkan dapat membantu para pengambil keputusan dalam menentukan strategi untuk meningkatkan perekonomian.

Package

Berikut ini, packages yang digunakan dalam analisis ini menggunakan bantuan software R.

#Load Package
library(readxl)
library(ggplot2)
library(forecast)

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

library(TTR)
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(vip)

## 
## Attaching package: 'vip'

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

library(tseries)
library(FinTS)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

## 
## Attaching package: 'FinTS'

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

Data

Data yang digunakan adalah data Nilai Impor (Juta USD) bulanan pada rentang Januari 2005 - September 2022 dengan jumlah 213 amatan yang bersumber dari Tabel Dinamis Subjek Ekspor-Impor BPS https://www.bps.go.id

Nilai impor merupakan nilai suatu kelompok barang yang diimpor (melewati batas negara Indonesia) tanpa menggunakan dokumen non PEB/PIB, dalam satuan US Dolar(USD)

Input Data

setwd("C:/Users/ASUS/OneDrive/Documents/Kuliah/Semester 5/Metode Peramalan Deret Waktu/Project MPDW")
data1 <- read_excel("NILAI IMPOR 2005-2022.xlsx", sheet=2)
impor <- data1$Impor

# Membentuk Objek Time Series Impor
impor.ts<-ts(data1$Impor, frequency=12, start=2005)

Eksplorasi Data

Time Series Plot

plot(impor.ts,
     col = "navyblue",
     lwd = 1,
     type = "o",
     xlab = "Periode",
     ylab = "Nilai Impor (Juta USD)",
     main = "Time Series Plot Nilai Impor 2005 - 2022")

Plot time series di atas menunjukkan pergerakan nilai impor bulanan nasional dalam jutaan USD sejak 2005 hingga 2022 yang selalu berubah baik diamati berdasarkan bulan maupun tahun. Time series plot tersebut juga menunjukkan bahwa nilai Impor terlihat tidak stasioner baik dalam nilai tengah maupun ragam dan cenderung terlihat memiliki pola siklik. Nilai Impor ini dapat saja dipengaruhi oleh fluktuasi ekonomi jangka panjang seperti yang berhubungan dengan siklus bisnis.

datacomponents = decompose(impor.ts)
plot(datacomponents)

Grafik di atas menyajikan informasi yang lebih detail mengenai data time series impor 2010-2017 (Training), data ini memiliki trend yang tidak linear melainkan secara fluktuatif turun lalu naik secara lokal dan naik secara global. Adapun secara musiman, data impor menghasilkan pola musiman yang relatif seragam setiap tahunnya. Sedangkan berdasarkan plot sisaan, sisaan menyebar pada nilai tengah 0 dengan ragam yang tidak homogen.

Deskriptif

summary(impor)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4091    9654   12626   12213   15309   22151

Smoothing

Splitting Training dan Testing

Analisis deret waktu yang baik dibangun dengan membagi data terklebih dahulu menjadi dua kelompok, yakni data training dan data testing. Melalui data training, model analisis akan dibentuk atau dibangun yang kemudian akan diberlakukan pada data testing. Pembagian data menjadi training dan testing hendaknya memperhatikan pola data dan persentase training dan testing. Splitting yang baik menghasilkan data training dengan pola yang mirip atau serupa dengan testingnya. Untuk data impor yang terdiri dari 213 data. Data untuk training terdiri atas 144 amatan (2005-2016) dan 69 amatan untuk testing (2017-2012).

training<-data1[1:144,3]
testing<-data1[145:213,3]

# Membentuk Objek Time Series
training.ts<-ts(training, frequency=12, start=2005)
testing.ts<-ts(testing,frequency=12, start=2017)

## Plot Data Training dan testing
ts.plot(impor.ts, xlab = "Periode", ylab ="Nilai Impor (Juta USD)", 
        main = "Plot Data Training dan Data Testing")
lines(training.ts, col = "blue")
lines(testing.ts, col="Red")
legend(2007,21000,c("Data Training","Data Testing"), 
       lty=8, col=c("blue","red"), cex=0.8)

Berdasarkan hasil eksplorasi data, data impor Januari 2005 hingga September 2022 memiliki pola data trend positif. Oleh karena itu, metode pemulusan yang dapat digunakan adalah model Double Moving Average dan Double Exponensial Smooting. Dengan demikian, pemulusan data impor ini akan menggunakan kedua metode tersebut.

Double Moving Average

Penentuan M Optimum

Dalam menentukan parameter m terbaik, iterasi diperlukan untuk melihat hasil akurasi saat m=i kemudian dibandingkan nilai SSE, MSE, MAPE,atau satuan akurasi lainnya untuk menentukan m berapakah yang menghasilkan error terkecil.

akurasi.dma <- function(m){
  #Pemulusan DMA dengan n=m
  data.sma<-SMA(training.ts, n = m)
  dma <- SMA(data.sma, n = m)
  At <- 2*data.sma - dma
  Bt <- 2/(m-1)*(data.sma - dma)
  data.dma<- At+Bt
  data.ramal2<- c(NA, data.dma)
  
  t = 1:5
  f = c()
  
  for (i in t) {
    f[i] = At[length(At)] + Bt[length(Bt)]*(i)
  }
  
  data.gab <- cbind(aktual2 = c(training.ts,rep(NA,5)), pemulusan1 = c(data.sma,rep(NA,5)),pemulusan2 = c(data.dma, rep(NA,5)),At = c(At, rep(NA,5)), Bt = c(Bt,rep(NA,5)),ramalan2 = c(data.ramal2, f[-1]))
 
  #Menghitung nilai keakuratan
  error.dma = training.ts-data.ramal2[1:length(training.ts)]
  SSE.dma = sum(error.dma[(2*m):length(training.ts)]^2)
  MSE.dma = mean(error.dma[(2*m):length(training.ts)]^2)
  MAPE.dma = mean(abs((error.dma[(2*m):length(training.ts)]/training.ts[(2*m):length(training.ts)])*100))
  
  akurasi.dma <- t(matrix(c(m, SSE.dma, MSE.dma, MAPE.dma)))
  colnames(akurasi.dma) <- c("m", "SSE", "MSE", "MAPE")
  akurasi.dma
}

tabel <- rbind(akurasi.dma(2), akurasi.dma(3), akurasi.dma(4),
               akurasi.dma(5), akurasi.dma(6), akurasi.dma(7),
               akurasi.dma(8), akurasi.dma(9), akurasi.dma(10),
               akurasi.dma(11))
tabel

##        m       SSE     MSE      MAPE
##  [1,]  2 309401369 2194336 10.439217
##  [2,]  3 225355885 1621265  9.431641
##  [3,]  4 216266266 1578586  9.324884
##  [4,]  5 214972409 1592388  9.574156
##  [5,]  6 227304531 1709057  9.687653
##  [6,]  7 265080125 2023512 10.638508
##  [7,]  8 298191688 2311563 11.482700
##  [8,]  9 327450009 2578347 12.130971
##  [9,] 10 357273415 2858187 12.402245
## [10,] 11 370520139 3012359 12.551286

Berdasarkan hasil tabel di atas, parameter m terbaik yang dapat digunakan dalam pemulusan DMA data Impor adalah m=4 sebab memiliki nilai MSE dan MAPE yang lebih kecil dibandingkan parameter m lainnya. Oleh karena itu, pemulusan DMA data Impor akan menggunakan parameter m=4.

Selanjutnya, pemodelan DMA dilakukan dengan menggunakan parameter m=4.

Pemulusan DMA

pemulusan_sma2 <- TTR::SMA(training.ts,  n=4)
df_dma <- TTR::SMA(pemulusan_sma2, n=4)
At <- 2*pemulusan_sma2-df_dma
Bt <- 2/3*(pemulusan_sma2-df_dma)
pemulusan_dma <- At+Bt
ramal_dma <- c(NA, pemulusan_dma)

fcast = c()
for(i in 1:69){
  fcast[i]=At[length(At)] + Bt[length(Bt)]*(i)
}

df_dma <- cbind(df_aktual=c(training.ts,rep(NA,69)), pemulusan_dma=c(pemulusan_dma,rep(NA,69)), ramal_dma=c(ramal_dma,fcast[-1]))

ts.plot(df_dma[,1], xlab="periode waktu", ylab="Yt", col="blue", lty=3)
points(df_dma[,1])
lines(df_dma[,2],col="red",lwd=2)
lines(c(df_dma[,3]), col= "green")
title("Plot Pemulusan DMA",cex.main=1,font.main=4 ,col.main="black")
legend("topleft", c("Data aktual","Pemulusan DMA","Peramalan DMA"),lty=0.7,col=c ("blue","red","Green"))

Plot hasil pemulusan di atas menunjukkan pemulusan dengan metode Double Moving Average (DMA) cukup baik memuluskan data sebab garis hasil pemulusan tidak jauh berbeda (menyimpang) dari keadaan data aktual.

Akurasi Pemulusan DMA Pada Data Training dan Testing

# Ukuran Keakuratan
error.dma_train <- df_dma[, 1] - df_dma[, 3]
error.dma_test <- testing.ts-fcast

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

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

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

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

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

akurasi <- data.frame(
          "Ukuran Keakuratan" = c("SSE", "MSE", "MAPE", "RMSE", "MAD"),
          "Akurasi Training" = c(SSE.dma_train, MSE.dma_train, MAPE.dma_train, RMSE.dma_train, MAD.dma_train),
          "Akurasi Testing" = c(SSE.dma_test, MSE.dma_test, MAPE.dma_test, RMSE.dma_test, MAD.dma_test))

akurasi

##   Ukuran.Keakuratan Akurasi.Training Akurasi.Testing
## 1               SSE     2.162663e+08    7.943187e+09
## 2               MSE     1.578586e+06    1.151186e+08
## 3              MAPE     9.324884e+00    6.362226e+01
## 4              RMSE     1.256418e+03    1.072934e+04
## 5               MAD     9.787806e+02    9.102710e+03

Double Exponensial Smooting

Penentuan Parameter α dan β Optimum

pemulusan2 <- function(x,alfa=NULL,beta=NULL,method=c("SES","DES"),
                       metrik=c("SSE","MSE","RMSE")){
  df.master <- data.frame()
  if(method=="SES"){
    for(i in alfa){
      df_ses <- HoltWinters(x, alpha = i, beta=F, gamma=F)
      sse <- df_ses$SSE 
      mse <- sse/length(x)
      rmse <-sqrt(mse)
      ak <- data.frame(Alpha=i,SSE=sse,MSE=mse,RMSE=rmse)
      df.master <- rbind(df.master,ak)
      datases <- data.frame(x, c(NA, df_ses$fitted[,1]))
      colnames(datases) = c("y","yhat")
    }
  }
  else if(method=="DES"){
    for (i in alfa) {
      for (j in beta) {
        df_des <- HoltWinters(x, alpha = i, beta=j, gamma=F)
        sse <- df_des$SSE 
        mse <- sse/length(x)
        rmse <-sqrt(mse)
        ak <- data.frame(Alpha=i,Beta=j,SSE=sse,MSE=mse,RMSE=rmse)
        df.master <- rbind(df.master,ak)
        datades <- data.frame(x, c(NA, NA, df_des$fitted[,1]))
        colnames(datades) = c("y","yhat")
      }
    }
  }
  if(method=="SES"){
    opt <- df.master$Alpha[which.min(df.master[,metrik])]
    return(list(Akurasi=df.master,n.optimum=paste("Alpha optimum adalah",opt)))
  }
  else if(method=="DES"){
    opt <- df.master$Alpha[which.min(df.master[,metrik])]
    opt2 <- df.master$Beta[which.min(df.master[,metrik])]
    return(list(Akurasi=df.master,n.optimum=paste("Alpha optimum adalah",opt,
                                                  ",Beta optimum adalah",opt2)))
  }
}
pemulusan2(training.ts,alfa=c(seq(0.1, 0.9, by = 0.1)),
           beta=c(seq(0.1, 0.9, by = 0.1)),method = "DES",metrik = "RMSE")

## $Akurasi
##    Alpha Beta       SSE     MSE     RMSE
## 1    0.1  0.1 411043344 2854468 1689.517
## 2    0.1  0.2 465321444 3231399 1797.609
## 3    0.1  0.3 502713807 3491068 1868.440
## 4    0.1  0.4 498204852 3459756 1860.042
## 5    0.1  0.5 504045089 3500313 1870.912
## 6    0.1  0.6 526639428 3657218 1912.385
## 7    0.1  0.7 556343976 3863500 1965.579
## 8    0.1  0.8 590834539 4103018 2025.591
## 9    0.1  0.9 624764322 4338641 2082.941
## 10   0.2  0.1 279482599 1940851 1393.144
## 11   0.2  0.2 294547445 2045468 1430.199
## 12   0.2  0.3 304300104 2113195 1453.683
## 13   0.2  0.4 314419482 2183469 1477.656
## 14   0.2  0.5 319639996 2219722 1489.873
## 15   0.2  0.6 317394271 2204127 1484.630
## 16   0.2  0.7 311369046 2162285 1470.471
## 17   0.2  0.8 306386363 2127683 1458.658
## 18   0.2  0.9 304149218 2112147 1453.323
## 19   0.3  0.1 220817691 1533456 1238.328
## 20   0.3  0.2 228316607 1585532 1259.179
## 21   0.3  0.3 232551990 1614944 1270.805
## 22   0.3  0.4 233532622 1621754 1273.481
## 23   0.3  0.5 232268355 1612975 1270.029
## 24   0.3  0.6 230531181 1600911 1265.271
## 25   0.3  0.7 228617276 1587620 1260.008
## 26   0.3  0.8 225989304 1569370 1252.745
## 27   0.3  0.9 222604189 1545862 1243.327
## 28   0.4  0.1 192219835 1334860 1155.361
## 29   0.4  0.2 197040259 1368335 1169.759
## 30   0.4  0.3 199422866 1384881 1176.810
## 31   0.4  0.4 200344548 1391282 1179.526
## 32   0.4  0.5 200796316 1394419 1180.855
## 33   0.4  0.6 201314116 1398015 1182.377
## 34   0.4  0.7 202326982 1405048 1185.347
## 35   0.4  0.8 204447316 1419773 1191.542
## 36   0.4  0.9 208235980 1446083 1202.532
## 37   0.5  0.1 178885178 1242258 1114.566
## 38   0.5  0.2 183559325 1274718 1129.034
## 39   0.5  0.3 186961716 1298345 1139.450
## 40   0.5  0.4 190120154 1320279 1149.034
## 41   0.5  0.5 193799504 1345830 1160.099
## 42   0.5  0.6 198559379 1378885 1174.259
## 43   0.5  0.7 204828404 1422419 1192.652
## 44   0.5  0.8 212829893 1477985 1215.724
## 45   0.5  0.9 222549467 1545482 1243.174
## 46   0.6  0.1 174603964 1212528 1101.148
## 47   0.6  0.2 180514792 1253575 1119.632
## 48   0.6  0.3 186119432 1292496 1136.880
## 49   0.6  0.4 192317842 1335541 1155.656
## 50   0.6  0.5 199691494 1386746 1177.602
## 51   0.6  0.6 208562733 1448352 1203.475
## 52   0.6  0.7 219010069 1520903 1233.249
## 53   0.6  0.8 230897972 1603458 1266.277
## 54   0.6  0.9 243948561 1694087 1301.571
## 55   0.7  0.1 176399258 1224995 1106.795
## 56   0.7  0.2 184312847 1279950 1131.349
## 57   0.7  0.3 192646500 1337823 1156.643
## 58   0.7  0.4 202108082 1403528 1184.706
## 59   0.7  0.5 213052546 1479532 1216.360
## 60   0.7  0.6 225568015 1566445 1251.577
## 61   0.7  0.7 239553837 1663568 1289.794
## 62   0.7  0.8 254821811 1769596 1330.262
## 63   0.7  0.9 271202614 1883351 1372.353
## 64   0.8  0.1 182904959 1270173 1127.020
## 65   0.8  0.2 193399321 1343051 1158.901
## 66   0.8  0.3 204907426 1422968 1192.882
## 67   0.8  0.4 217973982 1513708 1230.328
## 68   0.8  0.5 232824109 1616834 1271.548
## 69   0.8  0.6 249499551 1732636 1316.296
## 70   0.8  0.7 267982220 1860988 1364.180
## 71   0.8  0.8 288305546 2002122 1414.964
## 72   0.8  0.9 310627948 2157139 1468.720
## 73   0.9  0.1 193711506 1345219 1159.836
## 74   0.9  0.2 207474774 1440797 1200.332
## 75   0.9  0.3 222894169 1547876 1244.137
## 76   0.9  0.4 240463516 1669886 1292.241
## 77   0.9  0.5 260450501 1808684 1344.873
## 78   0.9  0.6 283063408 1965718 1402.041
## 79   0.9  0.7 308578023 2142903 1463.866
## 80   0.9  0.8 337423052 2343216 1530.757
## 81   0.9  0.9 370222608 2570990 1603.431
## 
## $n.optimum
## [1] "Alpha optimum adalah 0.6 ,Beta optimum adalah 0.1"

Berdasrkan fungsi optimasi tersebut, parameter Alpha optimum adalah 0.6 dan parameter Beta optimum bagi DES adalah 0.1. Parameter inilah yang kemudian digunakan dalam pemulusan Double Exponensial Smooting

DES <- HoltWinters(training.ts, gamma = FALSE, beta = 0.1, alpha = 0.6)
plot(DES, xlab= "Periode", ylab="Nilai Impor (Juta USD)", 
     main = "Pemulusan menggunakan DES (alpha = 0.6 dan beta = 0.1)",
     lwd=2)
legend(2006,17000,c("data aktual","data pemulusan"), lty=8, col=c("black","red"), cex=0.8)

Forecast

ramalandes<- forecast(DES, h=69)
ramalandes$mean

##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2017 12645.05 12734.56 12824.07 12913.57 13003.08 13092.59 13182.10 13271.61
## 2018 13719.15 13808.66 13898.16 13987.67 14077.18 14166.69 14256.20 14345.71
## 2019 14793.25 14882.75 14972.26 15061.77 15151.28 15240.79 15330.29 15419.80
## 2020 15867.34 15956.85 16046.36 16135.87 16225.38 16314.88 16404.39 16493.90
## 2021 16941.44 17030.95 17120.46 17209.97 17299.47 17388.98 17478.49 17568.00
## 2022 18015.54 18105.05 18194.56 18284.06 18373.57 18463.08 18552.59 18642.10
##           Sep      Oct      Nov      Dec
## 2017 13361.12 13450.62 13540.13 13629.64
## 2018 14435.21 14524.72 14614.23 14703.74
## 2019 15509.31 15598.82 15688.33 15777.84
## 2020 16583.41 16672.92 16762.43 16851.93
## 2021 17657.51 17747.01 17836.52 17926.03
## 2022 18731.60

Akurasi Pemulusan DES

selisihdestest<-ramalandes$mean-testing.ts

#SSE
SSEtraindes <- DES$SSE
SSEtestingdes<-sum(selisihdestest^2)

#RMSE
rmse.des.train <- sqrt(DES$SSE/length(training.ts))
rmse.des.test <- sqrt(mean((selisihdestest)^2))

#MAPE
MAPEtraindes <- mean(abs((DES$fitted[,1]-training.ts)/training.ts)*100)
MAPEtestingdes <- sum(abs(selisihdestest/testing.ts))*100/69

akurasi <- matrix(c(SSEtraindes,SSEtestingdes,rmse.des.train,rmse.des.test, MAPEtraindes, MAPEtestingdes), 
                  nrow=2, ncol=3)
ab <- matrix(c("Training", "Testing"))
akurasi <- cbind(ab,akurasi)
colnames(akurasi) <- c("Data","SSE", "RMSE", "MAPE")
akurasi <- data.frame(akurasi)
akurasi

##       Data              SSE             RMSE             MAPE
## 1 Training 174603963.562378 1101.14827554637 8.48539767889533
## 2  Testing 506662084.251334 2709.78389493807 15.7376262496705

Berdasarkan tabel akurasi di atas, DES yang dibangun saat alpha sama dengan 0.6 dan beta sama dengan 0.1 menghasilkan MAPE yang terbilang rendah baik pada data training maupun pada data testing.

Berdasarkan tabel di atas, pemulusan DES merupakan pemulusan terbaik bagi data impor karena menghasilkan SSE, RMSE, dan MAPE yang lebih kecil dibandingkan DMA.

Pemodelan SARIMA

Transformasi Logaritma

ln.impor <- log(data1[,3])
data1 <- data.frame(data1,ln.impor)

Splitting Training dan Testing Data Transformasi

Analisis deret waktu yang baik dibangun dengan membagi data terklebih dahulu menjadi dua kelompok, yakni data training dan data testing. Melalui data training, model analisis akan dibentuk atau dibangun yang kemudian akan diberlakukan pada data testing. Untuk data impor yang terdiri dari 213 data. Data untuk training terdiri atas 144 amatan (Januari 2005- Desember 2016) dan 69 amatan untuk testing (Januari 2017 - September 2022).

training<-data1[1:144,5]
testing<-data1[145:213,5]

# Membentuk Objek Time Series
training.ts<-ts(training, frequency=12, start=2005)
testing.ts<-ts(testing,frequency=12, start=2017)

## Plot Data Training dan testing
impor1.ts<-ts(data1$Impor.1, frequency=12, start=2005)
ts.plot(impor1.ts, xlab = "Periode", ylab ="Nilai Impor (Juta USD)", 
        main = "Plot Data Training dan Data Testing")
lines(training.ts, col = "blue")
lines(testing.ts, col="Red")
legend(2006,9.8,c("Data Training","Data Testing"), 
       lty=8, col=c("blue","red"), cex=0.8)

Berdasarkan plot time series di atas, data training dan data testing terlihat memiliki pola yang tidak berbeda jauh dan masih dapat dikategorikan berpola sama. Dengan demikian, splitting data kita cukup baik karena tetap menjaga pola data pada data training dan testing.

Uji Kestasioneran

Plot ACF

acf.impor <- acf(training.ts, lag=50,xaxt="n")
axis(1, at=0:50/12, labels=0:50)

Berdasarkan plot ACF di atas terlihat bahwa antar lag menurun secara perlahan. Hal ini mengindikasikan bahwa data tidak stasioner.

1. Uji Kestasioneran Dalam Rataan (Augmented Dickey - Fuller)

Uji kestasioneran dalam rataan dapat menggunakan uji Augmented Dickey - Fuller.

Hipotesis yang digunakan dalam uji ini yaitu:

H0: Data tidak stasioner

H1: Data stasioner

Penolakan terhadap H0 dapat diputuskan apabila p-value yang dihasilkan uji Augmented Dickey-Fuller lebih kecil dibandingkan alpha (5%).

Hasil perhitungan ADF dapat dilihat sebagai berikut.

#Uji formal - Augmented Dickey-Fuller Test
adf.test(training.ts)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  training.ts
## Dickey-Fuller = -1.8707, Lag order = 5, p-value = 0.6302
## alternative hypothesis: stationary

Karena nilai p-value = 0.6302 > 0.05, maka belum cukup bukti untuk menolak H0. Oleh karena itu dapat dikatakan bahwa data memang tidak stasioner dalam rataan.

2. Uji Kestasioneran Dalam Ragam (Barlett)

Uji kestasioneran dalam rataan dapat menggunakan uji Barlet.

Hipotesis yang digunakan dalam uji ini yaitu:

H0: Data deret waktu stasioner

H1: Data deret waktu tidak stasioner

Penolakan terhadap H0 dapat diputuskan apabila p-value yang dihasilkan uji Barlett lebih kecil dibandingkan alpha (5%).

Hasil perhitungan ADF dapat dilihat sebagai berikut.

#Uji formal - Barlett
bartlett.test(Impor.1~Tahun, data=data1)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  Impor.1 by Tahun
## Bartlett's K-squared = 34.737, df = 17, p-value = 0.006734

Karena nilai p-value = 0.006734 < 0.05, maka cukup bukti untuk menolak H0. Oleh karena itu dapat dikatakan bahwa data tidak stasioner dalam ragam.

Penanganan Ketakstasioneran

Untuk membuat kesimpulan statistik tentang struktur proses stokastik berdasarkan catatan yang diamati dari proses itu, beberapa asumsi penyederhanaan (dan mungkin masuk akal) tentang struktur terkadang diperlukan. Salah satu asumsi paling penting adalah asumsi stasioneritas data (Cyrer JD dan Chan KS. 2008). Ketidakstasionaran data impor 2010-2019 di atas harus ditangani. Salah satu metode untuk menangani masalah ini adalah metode differencing. Hasil differencing pertama data training dapat dilihat di bawah ini.

#differencing pertama
train_ts_diff <- diff(training.ts, differences=1)

#identifikasi stasioneritas komponen non musiman
plot.ts(train_ts_diff,lty=1,xlab="Periode", ylab="Nilai Impor Differencing 1")

acf1 <- acf(train_ts_diff, lag.max = 60, xaxt="n")
axis(1, at=0:60/12, labels=0:60)

adf.test(train_ts_diff)

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

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

#uji homogenitas dengan kelompok = tahun (data differencing 1)
training.diff<-as.data.frame(train_ts_diff)
tahun_training <- data.frame("tahun"=data1[2:144,4])
data2 <- data.frame("training"= training.diff$x,"tahun_training"=tahun_training$tahun)
bartlett.test(training ~ tahun_training, data=data2)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  training by tahun_training
## Bartlett's K-squared = 10.315, df = 11, p-value = 0.5023

Pada plot time-series terlihat bahwa data yang sudah dilakukan differencing sudah menyebar di sekitar suatu konstanta sumbu y, sehingga dapat dikatakan bahwa data sudah terlihat stasioner. Akan tetapi, dari plot tersebut masih diragukan apabila ragam pada data tersebut sudah homogen.

Berdasarkan plot ACF terlihat bahwa autokorelasi antar lag sudah tidak menurun secara perlahan namun berupa tails off. Begitu pula dengan hasil pengujian ADF memperlihatkan bahwa data sudah stasioner (p-value < 5%). Selain itu, uji Bartlett juga menunjukkan bahwa ragam data sudah homogen, yang artinya data sudah stasioner (p-value > 5%). Oleh karena itu, perlu dilakukan transformasi logaritma pada data untuk membuat ragam homogen.

#identifikasi stasioneritas pada komponen musiman
acf1$lag <- acf1$lag * 12
     acf1.1 <- as.data.frame(cbind(acf1$acf,acf1$lag))
     acf1.2 <- acf1.1[which(acf1.1$V2%%12==0),]
     barplot(height=acf1.2$V1,names.arg=acf1.2$V2,ylab="ACF",xlab="lag")

Berdasarkan plot ACF di atas terlihat bahwa komponen musiman tidak menurun secara perlahan yang artinya sudah stasioner. Oleh karena itu, kedua komponen (musiman dan non musiman) sudah stasioner.

Identifikasi/Spesifikasi Model

Komponen Non Musiman

Penentuan parameter model ARIMA dapat dilihat dari beberapa sudut pandang. Untuk parameter p dan q dapat dilihat dari plot ACF, plot PACF, dan plot EACF. Apabila dilakukan differencing ordo 1, maka nilai $d=1$, dan seterusnya.

Pemodelan akan didasarkan pada data training, kemudian validasi model akan didasarkan pada testing. Karena data dalam keadaan tidak stasioner, maka data training dan testing akan dikenakan proses differencing. Kemudian, plot ACF, PACF, dan EACF dari data training dapat dilihat sebagai berikut.

acf.diff <- acf(train_ts_diff, lag.max = 60, xaxt="n")
axis(1, at=0:60/12, labels=0:60)

pacf.diff <- 
  pacf(train_ts_diff, lag.max = 60, xaxt="n") 
axis(1, at=0:60/12, labels=0:60)

eacf(train_ts_diff)

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

Berdasarkan plot ACF (cut off setelah lag ke-1), PACF (tails off), dan EACF, diperoleh kandidat model untuk komponen non musiman yaitu

ARIMA(0,1,1)
ARIMA(1,1,0)
ARIMA(1,1,1)
ARIMA(2,1,1)
ARIMA(3,1,1)

Komponen Musiman

#ACF
acf.diff$lag <- acf.diff$lag * 12
acf.diff <- as.data.frame(cbind(acf.diff$acf,acf.diff$lag))
acf.difff <- acf.diff[which(acf.diff$V2%%12==0),]
barplot(height = acf.difff$V1, names.arg=acf.difff$V2, ylab="ACF", xlab="Lag")

# PACF
pacf.diff$lag <- pacf.diff$lag * 12
pacf.diff <- as.data.frame(cbind(pacf.diff$acf,pacf.diff$lag))
pacf.difff <- pacf.diff[which(pacf.diff$V2%%12==0),]
barplot(height = pacf.difff$V1, names.arg=pacf.difff$V2, ylab="PACF", xlab="Lag")

Berdasarkan plot ACF dan PACF diperoleh kandidat model untuk komponen musiman adalah

ARIMA(0,0,1)12
ARIMA(1,0,0)12
ARIMA(1,0,1)12

Oleh karena itu, model tentatif yang diperoleh berdasarkan hasil identifikasi menggunakan plot ACF, PACF, dan EACF adalah:

ARIMA(0,1,1)x(0,0,1)12
ARIMA(1,1,0)x(0,0,1)12
ARIMA(1,1,1)x(0,0,1)12
ARIMA(2,1,1)x(0,0,1)12
ARIMA(3,1,1)x(0,0,1)12
ARIMA(0,1,1)x(1,0,0)12
ARIMA(1,1,0)x(1,0,0)12
ARIMA(1,1,1)x(1,0,0)12
ARIMA(2,1,1)x(1,0,0)12

10.ARIMA(3,1,1)x(1,0,0)12

11.ARIMA(0,1,1)x(1,0,1)12

12.ARIMA(1,1,0)x(1,0,1)12

13.ARIMA(1,1,1)x(1,0,1)12

14.ARIMA(2,1,1)x(1,0,1)12

15.ARIMA(3,1,1)x(1,0,1)12

Pendugaan Parameter

printstatarima <- function (x, digits = 4,se=TRUE){
       if (length(x$coef) > 0) {
         cat("\nCoefficients:\n")
         coef <- round(x$coef, digits = digits)
         if (se && nrow(x$var.coef)) {
           ses <- rep(0, length(coef))
           ses[x$mask] <- round(sqrt(diag(x$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(x$residuals)-1, lower.tail = FALSE)
           coef <- rbind(coef, t=round(statt,digits=digits),sign.=round(pval,digits=digits))
           coef <- t(coef)
         }
         print.default(coef, print.gap = 2)
       }
     }

1. ARIMA(0,1,1)x(0,0,1)12

mod1.impor <- Arima(training.ts, order = c(0,1,1), seasonal=c(0,0,1))
printstatarima(mod1.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -0.2687  0.0695  -3.8662  0.0002
## sma1   0.2098  0.0788   2.6624  0.0086

summary(mod1.impor)

## Series: training.ts 
## ARIMA(0,1,1)(0,0,1)[12] 
## 
## Coefficients:
##           ma1    sma1
##       -0.2687  0.2098
## s.e.   0.0695  0.0788
## 
## sigma^2 = 0.01049:  log likelihood = 123.67
## AIC=-241.34   AICc=-241.17   BIC=-232.45
## 
## Training set error measures:
##                       ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.009038194 0.1013284 0.0787431 0.09446829 0.8619154 0.3706507
##                     ACF1
## Training set -0.06032249

2. ARIMA(1,1,0)x(0,0,1)12

mod2.impor <- Arima(training.ts, order = c(1,1,0), seasonal=c(0,0,1))
printstatarima(mod2.impor)

## 
## Coefficients:
##                  s.e.        t  sign.
## ar1   -0.3376  0.0790  -4.2734  0.000
## sma1   0.2106  0.0769   2.7386  0.007

summary(mod2.impor)

## Series: training.ts 
## ARIMA(1,1,0)(0,0,1)[12] 
## 
## Coefficients:
##           ar1    sma1
##       -0.3376  0.2106
## s.e.   0.0790  0.0769
## 
## sigma^2 = 0.01021:  log likelihood = 125.55
## AIC=-245.09   AICc=-244.92   BIC=-236.2
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.008883085 0.09999139 0.07813673 0.09293898 0.8546978 0.3677965
##                   ACF1
## Training set 0.0151872

3. ARIMA(1,1,1)x(0,0,1)12

mod3.impor <- Arima(training.ts, order = c(1,1,1), seasonal=c(0,0,1))
printstatarima(mod3.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.4696  0.1909  -2.4599  0.0151
## ma1    0.1497  0.2121   0.7058  0.4815
## sma1   0.2206  0.0780   2.8282  0.0054

summary(mod3.impor)

## Series: training.ts 
## ARIMA(1,1,1)(0,0,1)[12] 
## 
## Coefficients:
##           ar1     ma1    sma1
##       -0.4696  0.1497  0.2206
## s.e.   0.1909  0.2121  0.0780
## 
## sigma^2 = 0.01024:  log likelihood = 125.79
## AIC=-243.59   AICc=-243.3   BIC=-231.74
## 
## Training set error measures:
##                       ME       RMSE       MAE        MPE      MAPE      MASE
## Training set 0.008444578 0.09979705 0.0784178 0.08821354 0.8573058 0.3691194
##                      ACF1
## Training set -0.003184908

4. ARIMA(2,1,1)x(0,0,1)12

mod4.impor <- Arima(training.ts, order = c(2,1,1), seasonal=c(0,0,1))
printstatarima(mod4.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1    0.2138  0.2883   0.7416  0.4596
## ar2    0.2675  0.1104   2.4230  0.0166
## ma1   -0.5225  0.2878  -1.8155  0.0715
## sma1   0.2303  0.0773   2.9793  0.0034

summary(mod4.impor)

## Series: training.ts 
## ARIMA(2,1,1)(0,0,1)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1    sma1
##       0.2138  0.2675  -0.5225  0.2303
## s.e.  0.2883  0.1104   0.2878  0.0773
## 
## sigma^2 = 0.01023:  log likelihood = 126.38
## AIC=-242.76   AICc=-242.32   BIC=-227.94
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.007206821 0.09936332 0.07874904 0.07490873 0.8609939 0.3706786
##                     ACF1
## Training set -0.01671797

5. ARIMA(3,1,1)x(0,0,1)12

mod5.impor <- Arima(training.ts, order = c(3,1,1), seasonal=c(0,0,1))
printstatarima(mod5.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1    0.0739  0.4061   0.1820  0.8559
## ar2    0.2230  0.1436   1.5529  0.1227
## ar3    0.0581  0.0960   0.6052  0.5460
## ma1   -0.3937  0.3997  -0.9850  0.3263
## sma1   0.2266  0.0772   2.9352  0.0039

summary(mod5.impor)

## Series: training.ts 
## ARIMA(3,1,1)(0,0,1)[12] 
## 
## Coefficients:
##          ar1     ar2     ar3      ma1    sma1
##       0.0739  0.2230  0.0581  -0.3937  0.2266
## s.e.  0.4061  0.1436  0.0960   0.3997  0.0772
## 
## sigma^2 = 0.01028:  log likelihood = 126.54
## AIC=-241.09   AICc=-240.47   BIC=-223.31
## 
## Training set error measures:
##                       ME       RMSE       MAE        MPE      MAPE      MASE
## Training set 0.007086245 0.09925194 0.0782665 0.07362963 0.8561099 0.3684073
##                      ACF1
## Training set -0.007573443

6. ARIMA(0,1,1)x(1,0,0)12

mod6.impor <- Arima(training.ts, order = c(0,1,1), seasonal=c(1,0,0))
printstatarima(mod6.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -0.2706  0.0683  -3.9619  0.0001
## sar1   0.2423  0.0865   2.8012  0.0058

summary(mod6.impor)

## Series: training.ts 
## ARIMA(0,1,1)(1,0,0)[12] 
## 
## Coefficients:
##           ma1    sar1
##       -0.2706  0.2423
## s.e.   0.0683  0.0865
## 
## sigma^2 = 0.0104:  log likelihood = 124.17
## AIC=-242.34   AICc=-242.16   BIC=-233.45
## 
## Training set error measures:
##                       ME      RMSE        MAE        MPE     MAPE      MASE
## Training set 0.008313194 0.1009111 0.07821074 0.08669396 0.856158 0.3681448
##                     ACF1
## Training set -0.06251728

7. ARIMA(1,1,0)x(1,0,0)12

mod7.impor <- Arima(training.ts, order = c(1,1,0), seasonal=c(1,0,0))
printstatarima(mod7.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.3444  0.0783  -4.3985  0.0000
## sar1   0.2519  0.0852   2.9566  0.0036

summary(mod7.impor)

## Series: training.ts 
## ARIMA(1,1,0)(1,0,0)[12] 
## 
## Coefficients:
##           ar1    sar1
##       -0.3444  0.2519
## s.e.   0.0783  0.0852
## 
## sigma^2 = 0.01009:  log likelihood = 126.24
## AIC=-246.49   AICc=-246.32   BIC=-237.6
## 
## Training set error measures:
##                       ME       RMSE       MAE        MPE      MAPE      MASE
## Training set 0.008095158 0.09941876 0.0774354 0.08448085 0.8471548 0.3644952
##                   ACF1
## Training set 0.0200676

8. ARIMA(1,1,1)x(1,0,0)12

mod8.impor <- Arima(training.ts, order = c(1,1,1), seasonal=c(1,0,0))
printstatarima(mod8.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.4849  0.1834  -2.6439  0.0091
## ma1    0.1594  0.2040   0.7814  0.4359
## sar1   0.2641  0.0864   3.0567  0.0027

summary(mod8.impor)

## Series: training.ts 
## ARIMA(1,1,1)(1,0,0)[12] 
## 
## Coefficients:
##           ar1     ma1    sar1
##       -0.4849  0.1594  0.2641
## s.e.   0.1834  0.2040  0.0864
## 
## sigma^2 = 0.01012:  log likelihood = 126.55
## AIC=-245.09   AICc=-244.8   BIC=-233.24
## 
## Training set error measures:
##                       ME       RMSE        MAE       MPE      MAPE      MASE
## Training set 0.007605038 0.09917789 0.07771942 0.0791964 0.8497753 0.3658321
##                       ACF1
## Training set -0.0004801932

9. ARIMA(2,1,1)x(1,0,0)12

mod9.impor <- Arima(training.ts, order = c(2,1,1), seasonal=c(1,0,0))
printstatarima(mod9.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1    0.2054  0.2638   0.7786  0.4375
## ar2    0.2814  0.1062   2.6497  0.0090
## ma1   -0.5188  0.2630  -1.9726  0.0505
## sar1   0.2805  0.0867   3.2353  0.0015

summary(mod9.impor)

## Series: training.ts 
## ARIMA(2,1,1)(1,0,0)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1    sar1
##       0.2054  0.2814  -0.5188  0.2805
## s.e.  0.2638  0.1062   0.2630  0.0867
## 
## sigma^2 = 0.01008:  log likelihood = 127.28
## AIC=-244.55   AICc=-244.11   BIC=-229.74
## 
## Training set error measures:
##                       ME       RMSE        MAE       MPE    MAPE      MASE
## Training set 0.006259876 0.09862273 0.07826242 0.0647195 0.85582 0.3683881
##                     ACF1
## Training set -0.01654813

10.ARIMA(3,1,1)x(1,0,0)12

mod10.impor <- Arima(training.ts, order = c(3,1,1), seasonal=c(1,0,0))
printstatarima(mod10.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1    0.0612  0.3791   0.1614  0.8720
## ar2    0.2356  0.1375   1.7135  0.0888
## ar3    0.0645  0.0967   0.6670  0.5058
## ma1   -0.3868  0.3730  -1.0370  0.3015
## sar1   0.2786  0.0866   3.2171  0.0016

summary(mod10.impor)

## Series: training.ts 
## ARIMA(3,1,1)(1,0,0)[12] 
## 
## Coefficients:
##          ar1     ar2     ar3      ma1    sar1
##       0.0612  0.2356  0.0645  -0.3868  0.2786
## s.e.  0.3791  0.1375  0.0967   0.3730  0.0866
## 
## sigma^2 = 0.01012:  log likelihood = 127.48
## AIC=-242.95   AICc=-242.34   BIC=-225.18
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.00613422 0.0984846 0.0778802 0.06338274 0.8521004 0.3665889
##                      ACF1
## Training set -0.006035975

11.ARIMA(0,1,1)x(1,0,1)12

mod11.impor <- Arima(training.ts, order = c(0,1,1), seasonal=c(1,0,1))
printstatarima(mod11.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ma1   -0.2756  0.0704  -3.9148  0.0001
## sar1   0.3288  0.2978   1.1041  0.2714
## sma1  -0.0919  0.3117  -0.2948  0.7685

summary(mod11.impor)

## Series: training.ts 
## ARIMA(0,1,1)(1,0,1)[12] 
## 
## Coefficients:
##           ma1    sar1     sma1
##       -0.2756  0.3288  -0.0919
## s.e.   0.0704  0.2978   0.3117
## 
## sigma^2 = 0.01047:  log likelihood = 124.21
## AIC=-240.43   AICc=-240.14   BIC=-228.58
## 
## Training set error measures:
##                       ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.008134257 0.1008739 0.0781378 0.08478943 0.8553199 0.3678015
##                     ACF1
## Training set -0.06431237

12.ARIMA(1,1,0)x(1,0,1)12

mod12.impor <- Arima(training.ts, order = c(1,1,0), seasonal=c(1,0,1))
printstatarima(mod12.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.3609  0.0826  -4.3692  0.0000
## sar1   0.4459  0.2964   1.5044  0.1347
## sma1  -0.2059  0.3209  -0.6416  0.5221

summary(mod12.impor)

## Series: training.ts 
## ARIMA(1,1,0)(1,0,1)[12] 
## 
## Coefficients:
##           ar1    sar1     sma1
##       -0.3609  0.4459  -0.2059
## s.e.   0.0826  0.2964   0.3209
## 
## sigma^2 = 0.01013:  log likelihood = 126.45
## AIC=-244.91   AICc=-244.62   BIC=-233.05
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.007559821 0.09924023 0.07701405 0.07878029 0.8424047 0.3625119
##                    ACF1
## Training set 0.01946021

13.ARIMA(1,1,1)x(1,0,1)12

mod13.impor <- Arima(training.ts, order = c(1,1,1), seasonal=c(1,0,1))
printstatarima(mod13.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1   -0.4824  0.1769  -2.7270  0.0072
## ma1    0.1424  0.1993   0.7145  0.4761
## sar1   0.4248  0.2815   1.5091  0.1335
## sma1  -0.1722  0.3044  -0.5657  0.5725

summary(mod13.impor)

## Series: training.ts 
## ARIMA(1,1,1)(1,0,1)[12] 
## 
## Coefficients:
##           ar1     ma1    sar1     sma1
##       -0.4824  0.1424  0.4248  -0.1722
## s.e.   0.1769  0.1993  0.2815   0.3044
## 
## sigma^2 = 0.01016:  log likelihood = 126.71
## AIC=-243.42   AICc=-242.98   BIC=-228.61
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE     MASE
## Training set 0.007208284 0.09904052 0.07725627 0.07496752 0.8446906 0.363652
##                       ACF1
## Training set -0.0006875571

14.ARIMA(2,1,1)x(1,0,1)12

mod14.impor <- Arima(training.ts, order = c(2,1,1), seasonal=c(1,0,1))
printstatarima(mod14.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1    0.2010  0.2589   0.7764  0.4388
## ar2    0.2893  0.1084   2.6688  0.0085
## ma1   -0.5280  0.2584  -2.0433  0.0429
## sar1   0.4320  0.2599   1.6622  0.0987
## sma1  -0.1633  0.2813  -0.5805  0.5625

summary(mod14.impor)

## Series: training.ts 
## ARIMA(2,1,1)(1,0,1)[12] 
## 
## Coefficients:
##          ar1     ar2      ma1    sar1     sma1
##       0.2010  0.2893  -0.5280  0.4320  -0.1633
## s.e.  0.2589  0.1084   0.2584  0.2599   0.2813
## 
## sigma^2 = 0.01012:  log likelihood = 127.45
## AIC=-242.89   AICc=-242.28   BIC=-225.12
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.005918812 0.09847979 0.07790866 0.06106489 0.8518978 0.3667229
##                     ACF1
## Training set -0.01822714

15.ARIMA(3,1,1)x(1,0,1)12

mod15.impor <- Arima(training.ts, order = c(3,1,1), seasonal=c(1,0,1))
printstatarima(mod15.impor)

## 
## Coefficients:
##                  s.e.        t   sign.
## ar1    0.0434  0.3644   0.1191  0.9054
## ar2    0.2391  0.1383   1.7289  0.0860
## ar3    0.0728  0.0956   0.7615  0.4476
## ma1   -0.3873  0.3576  -1.0831  0.2806
## sar1   0.4566  0.2638   1.7309  0.0856
## sma1  -0.1919  0.2887  -0.6647  0.5073

summary(mod15.impor)

## Series: training.ts 
## ARIMA(3,1,1)(1,0,1)[12] 
## 
## Coefficients:
##          ar1     ar2     ar3      ma1    sar1     sma1
##       0.0434  0.2391  0.0728  -0.3873  0.4566  -0.1919
## s.e.  0.3644  0.1383  0.0956   0.3576  0.2638   0.2887
## 
## sigma^2 = 0.01016:  log likelihood = 127.71
## AIC=-241.41   AICc=-240.58   BIC=-220.67
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE      MASE
## Training set 0.005703391 0.09829342 0.07757298 0.05876581 0.8486729 0.3651428
##                      ACF1
## Training set -0.006030904

Hanya model ARIMA(0,1,1)x(0,0,1)12, ARIMA(1,1,0)x(0,0,1)12, ARIMA(0,1,1)x(1,0,0)12,dan ARIMA(1,1,0)x(1,0,0)12 yang memiliki seluruh parameter signifikan dalam model yang ditunjukkan oleh nilai pr(>|z|) < alpha (5%). Sehingga model-model tersebut yang akan dipertimbangkan untuk digunakan dalam peramalan.

Perbandingan MAPE dan AIC

Selanjutnya untuk mengetahui model yang terbaik diantara keenam model tersebut, perbandingan nilai AIC dan MAPE dapat digunakan untuk menentukannya model terbaik. Model dengan nilai mape terkecil dan aic terkecil merupakan model yang akan dipilih.

akurasi<-data.frame("Model ARIMA"=c("ARIMA(0,1,1)x(0,0,1)12",
                                    "ARIMA(1,1,0)x(0,0,1)12",
                                    "ARIMA(0,1,1)x(1,0,0)12",
                                    "ARIMA(1,1,0)x(1,0,0)12"),
                     "MAPE"=c(0.8619154,0.8546978,0.856158,0.8471548), 
                    "AIC"=c(-241.34,-245.09,-242.34,-246.49))

akurasi %>% knitr::kable(digits = 3)

Model.ARIMA	MAPE	AIC
ARIMA(0,1,1)x(0,0,1)12	0.862	-241.34
ARIMA(1,1,0)x(0,0,1)12	0.855	-245.09
ARIMA(0,1,1)x(1,0,0)12	0.856	-242.34
ARIMA(1,1,0)x(1,0,0)12	0.847	-246.49

Model ARIMA(1,1,0)x(1,0,0)12 memiliki MAPE dan AIC terkecil. Oleh karena itu, model ARIMA(1,1,0)x(1,0,0)12 yang akan digunakan.

Diagnostik Model

#diagnostik model
checkresiduals(mod7.impor)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)(1,0,0)[12]
## Q* = 30.816, df = 22, p-value = 0.09993
## 
## Model df: 2.   Total lags used: 24

tsdisplay(residuals(mod7.impor), lag.max=60, main='ARIMA(1,1,0)x(1,0,0)12 Model Residuals')

#Uji Formal
resid.mod7 <- mod7.impor$residuals
norm.resid <- tseries::jarque.bera.test(resid.mod7)
stat.ljungbox <- 30.816
p.value.ljungbox <- 0.09993

results <- data.frame("Statistics" = c(stat.ljungbox,
                                       norm.resid$statistic), 
                      "p-value" = c(p.value.ljungbox,
                                    norm.resid$p.value))
rownames(results) <- c( "Ljung-Box",
                        "Jarque bera")

results %>% knitr::kable(digits = 4)

	Statistics	p.value
Ljung-Box	30.8160	0.0999
Jarque bera	3.9458	0.1391

Secara eksploratif, berdasarkan plot ACF dan PACF terlihat bahwa tidak terdapat autokorelasi yang signifikan sehingga tidak ada gejala autokorelasi pada sisaan. Selanjutnya dilakukan pengujian menggunakan Ljung-Box Test dan diperoleh hasil p-value = 0.0999 > ?? = 0.05 yang berarti bahwa H0 tidak ditolak (tidak terdapat autokorelasi pada sisaan model).

Selanjutnya, apabila dilihat melalui histogram sisaan, terlihat bahwa sisaan menyebar normal, selain itu, berdasarkan pengujian menggunakan Jarque bera Test diperoleh p-value = 0.1391 > 0.05 yang berarti bahwa tak tolak H0 sehingga dapat disimpulkan sisaan model menyebar normal.

Overfitting

overfit1 <-  Arima(training.ts, order = c(1, 1, 0),seasonal=c(1,0,1) , include.drift = F, method = "ML")
overfit2 <-  Arima(training.ts, order = c(1, 1, 0), seasonal=c(2,0,0), include.drift = F, method = "ML")

lmtest::coeftest(overfit1)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.360864   0.082557 -4.3711 1.236e-05 ***
## sar1  0.445998   0.296383  1.5048    0.1324    
## sma1 -0.205985   0.320896 -0.6419    0.5209    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

lmtest::coeftest(overfit2)

## 
## z test of coefficients:
## 
##       Estimate Std. Error z value  Pr(>|z|)    
## ar1  -0.364761   0.081482 -4.4766 7.585e-06 ***
## sar1  0.233544   0.087318  2.6746  0.007481 ** 
## sar2  0.075018   0.095208  0.7879  0.430733    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Validasi Model

#validasi model
ramalan_arima = forecast(mod7.impor, 69)
ramalan_arima %>%knitr::kable(digits = 3, 
                                caption = "Hasil peramalan model ARIMA(1,1,0)x(1,0,0)12")

Hasil peramalan model ARIMA(1,1,0)x(1,0,0)12
	Point Forecast	Lo 80	Hi 80	Lo 95	Hi 95
Jan 2017	9.421	9.292	9.550	9.224	9.618
Feb 2017	9.413	9.259	9.567	9.178	9.649
Mar 2017	9.440	9.257	9.623	9.159	9.720
Apr 2017	9.429	9.222	9.635	9.113	9.744
May 2017	9.436	9.209	9.664	9.088	9.784
Jun 2017	9.457	9.210	9.704	9.079	9.835
Jul 2017	9.383	9.118	9.648	8.978	9.788
Aug 2017	9.463	9.181	9.745	9.032	9.894
Sep 2017	9.440	9.142	9.737	8.985	9.895
Oct 2017	9.444	9.132	9.757	8.966	9.922
Nov 2017	9.469	9.142	9.796	8.969	9.969
Dec 2017	9.471	9.130	9.811	8.950	9.992
Jan 2018	9.462	9.098	9.826	8.906	10.019
Feb 2018	9.460	9.078	9.842	8.876	10.045
Mar 2018	9.467	9.066	9.868	8.854	10.080
Apr 2018	9.464	9.046	9.882	8.824	10.104
May 2018	9.466	9.031	9.901	8.800	10.132
Jun 2018	9.471	9.020	9.923	8.781	10.162
Jul 2018	9.453	8.985	9.920	8.738	10.167
Aug 2018	9.473	8.990	9.955	8.735	10.210
Sep 2018	9.467	8.970	9.964	8.707	10.227
Oct 2018	9.468	8.957	9.979	8.686	10.250
Nov 2018	9.474	8.949	9.999	8.671	10.277
Dec 2018	9.475	8.936	10.013	8.651	10.298
Jan 2019	9.472	8.919	10.026	8.626	10.319
Feb 2019	9.472	8.904	10.040	8.604	10.340
Mar 2019	9.474	8.892	10.055	8.584	10.363
Apr 2019	9.473	8.878	10.068	8.563	10.383
May 2019	9.473	8.865	10.082	8.543	10.403
Jun 2019	9.475	8.854	10.096	8.525	10.425
Jul 2019	9.470	8.836	10.104	8.501	10.439
Aug 2019	9.475	8.829	10.121	8.487	10.463
Sep 2019	9.474	8.815	10.132	8.467	10.480
Oct 2019	9.474	8.804	10.144	8.449	10.499
Nov 2019	9.476	8.794	10.157	8.433	10.518
Dec 2019	9.476	8.782	10.169	8.415	10.536
Jan 2020	9.475	8.770	10.180	8.397	10.554
Feb 2020	9.475	8.758	10.192	8.379	10.571
Mar 2020	9.475	8.748	10.203	8.362	10.589
Apr 2020	9.475	8.736	10.214	8.345	10.605
May 2020	9.475	8.725	10.225	8.329	10.622
Jun 2020	9.476	8.715	10.236	8.312	10.639
Jul 2020	9.474	8.703	10.246	8.295	10.654
Aug 2020	9.476	8.694	10.257	8.280	10.671
Sep 2020	9.475	8.683	10.267	8.264	10.687
Oct 2020	9.475	8.673	10.278	8.249	10.702
Nov 2020	9.476	8.664	10.288	8.234	10.718
Dec 2020	9.476	8.654	10.298	8.218	10.733
Jan 2021	9.476	8.644	10.308	8.203	10.748
Feb 2021	9.476	8.634	10.318	8.188	10.763
Mar 2021	9.476	8.624	10.327	8.173	10.778
Apr 2021	9.476	8.615	10.337	8.159	10.793
May 2021	9.476	8.605	10.346	8.144	10.807
Jun 2021	9.476	8.596	10.356	8.130	10.822
Jul 2021	9.476	8.586	10.365	8.116	10.835
Aug 2021	9.476	8.578	10.374	8.102	10.850
Sep 2021	9.476	8.568	10.383	8.088	10.864
Oct 2021	9.476	8.560	10.392	8.074	10.877
Nov 2021	9.476	8.551	10.401	8.061	10.891
Dec 2021	9.476	8.542	10.410	8.048	10.904
Jan 2022	9.476	8.533	10.419	8.034	10.918
Feb 2022	9.476	8.525	10.427	8.021	10.931
Mar 2022	9.476	8.516	10.436	8.008	10.944
Apr 2022	9.476	8.507	10.444	7.995	10.957
May 2022	9.476	8.499	10.453	7.982	10.970
Jun 2022	9.476	8.491	10.461	7.969	10.983
Jul 2022	9.476	8.482	10.469	7.956	10.995
Aug 2022	9.476	8.474	10.478	7.944	11.008
Sep 2022	9.476	8.466	10.486	7.932	11.020

accuracy(ramalan_arima,testing.ts)

##                       ME       RMSE       MAE        MPE      MAPE      MASE
## Training set 0.008095158 0.09941876 0.0774354 0.08448085 0.8471548 0.3644952
## Test set     0.127612982 0.23535819 0.1924961 1.28739532 1.9887593 0.9060962
##                   ACF1 Theil's U
## Training set 0.0200676        NA
## Test set     0.6441583  1.457076

Terlihat bahwa MAPE dari proses peramalan tersebut relatif kecil, sehingga model tersebut adalah model yang baik untuk digunakan dalam peramalan nilai impor.

plot(ramalan_arima)
lines(testing.ts,col="red")
legend("topleft", legend = c("Data Transformasi","Peramalan Arima (1,1,0)(1,0,0)[12]"),lty=0.6,col= c("red","blue"))

Transformasi Balik

forecast_arima <- cbind(ramalan_arima$mean,ramalan_arima$lower,ramalan_arima$upper)
forecast_arima <- cbind(data1[145:213,1:5],exp(forecast_arima))
forecast_arima <- as.data.frame(forecast_arima)
forecast_arima%>%knitr::kable(digits = 3)

	No	Periode	Impor	Tahun	Impor.1	ramalan_arima$mean\| ramalan_arima$lower.80%	ramalan_arima$lower.95%\| ramalan_arima$upper.80%	ramalan_arima$upper.95%
145	145	2017-01-01	11973.8	2017	9.390	12342.82	10851.640	10136.630	14038.92	15029.19
146	146	2017-02-01	11359.4	2017	9.338	12250.96	10502.812	9680.763	14290.09	15503.54
147	147	2017-03-01	13283.2	2017	9.494	12580.71	10472.461	9503.446	15113.39	16654.42
148	148	2017-04-01	11950.6	2017	9.389	12441.06	10121.965	9074.805	15291.50	17056.02
149	149	2017-05-01	13772.5	2017	9.530	12534.98	9983.015	8849.675	15739.31	17754.97
150	150	2017-06-01	9991.6	2017	9.210	12797.16	9997.354	8772.467	16381.06	18668.33
151	151	2017-07-01	13889.8	2017	9.539	11884.75	9119.442	7926.464	15488.58	17819.69
152	152	2017-08-01	13509.2	2017	9.511	12873.77	9713.988	8368.574	17061.38	19804.33
153	153	2017-09-01	12788.2	2017	9.456	12579.14	9342.492	7981.319	16937.10	19825.64
154	154	2017-10-01	14249.2	2017	9.564	12637.54	9245.792	7836.083	17273.54	20381.04
155	155	2017-11-01	15113.5	2017	9.623	12947.55	9337.667	7854.074	17952.99	21344.22
156	156	2017-12-01	15104.5	2017	9.623	12976.57	9230.861	7707.959	18242.21	21846.42
157	157	2018-01-01	15309.4	2018	9.636	12862.67	8938.844	7372.503	18508.90	22441.25
158	158	2018-02-01	14185.5	2018	9.560	12838.49	8759.765	7154.941	18816.34	23036.77
159	159	2018-03-01	14463.6	2018	9.579	12924.66	8655.563	7000.349	19299.37	23862.66
160	160	2018-04-01	16162.3	2018	9.690	12888.38	8481.987	6796.904	19583.88	24439.10
161	161	2018-05-01	17662.9	2018	9.779	12912.81	8355.888	6636.317	19954.88	25125.49
162	162	2018-06-01	11267.9	2018	9.330	12980.31	8264.585	6507.752	20386.81	25890.44
163	163	2018-07-01	18297.1	2018	9.814	12740.72	7986.047	6236.532	20326.20	26028.25
164	164	2018-08-01	16818.1	2018	9.730	12999.84	8026.024	6217.716	21055.99	27179.74
165	165	2018-09-01	14610.1	2018	9.589	12924.25	7863.067	6044.316	21243.16	27635.28
166	166	2018-10-01	17667.6	2018	9.779	12939.34	7760.792	5920.836	21573.39	28277.53
167	167	2018-11-01	16901.8	2018	9.735	13018.57	7700.739	5832.036	22008.68	29060.71
168	168	2018-12-01	15364.0	2018	9.640	13025.91	7601.649	5715.951	22320.73	29684.36
169	169	2019-01-01	15005.2	2019	9.616	12997.02	7471.742	5573.785	22608.18	30306.60
170	170	2019-02-01	12465.1	2019	9.431	12990.86	7364.433	5453.193	22915.88	30947.45
171	171	2019-03-01	13746.6	2019	9.529	13012.77	7275.277	5347.779	23275.01	31664.01
172	172	2019-04-01	15399.2	2019	9.642	13003.56	7172.845	5235.020	23573.97	32300.25
173	173	2019-05-01	14606.7	2019	9.589	13009.76	7082.187	5132.870	23898.54	32974.51
174	174	2019-06-01	11495.4	2019	9.350	13026.86	7000.542	5039.145	24240.84	33676.15
175	175	2019-07-01	15518.5	2019	9.650	12965.87	6880.193	4919.466	24434.46	34173.18
176	176	2019-08-01	14169.4	2019	9.559	13031.79	6829.956	4851.595	24865.11	35004.48
177	177	2019-09-01	14263.4	2019	9.565	13012.66	6737.460	4755.166	25132.53	35609.57
178	178	2019-10-01	14759.1	2019	9.600	13016.49	6659.412	4670.454	25442.03	36276.77
179	179	2019-11-01	15340.5	2019	9.638	13036.52	6591.839	4594.429	25782.00	36990.61
180	180	2019-12-01	14506.8	2019	9.582	13038.37	6517.146	4514.711	26084.89	37654.47
181	181	2020-01-01	14268.7	2020	9.566	13031.08	6437.612	4432.032	26377.63	38314.02
182	182	2020-02-01	11548.1	2020	9.354	13029.52	6363.870	4354.899	26676.92	38983.32
183	183	2020-03-01	13352.2	2020	9.499	13035.05	6295.249	4282.326	26990.61	39677.65
184	184	2020-04-01	12535.2	2020	9.436	13032.73	6224.775	4209.623	27286.45	40348.51
185	185	2020-05-01	8438.6	2020	9.041	13034.30	6157.905	4140.396	27589.39	41033.00
186	186	2020-06-01	10760.3	2020	9.284	13038.61	6094.002	4074.152	27897.15	41727.77
187	187	2020-07-01	10464.3	2020	9.256	13023.20	6022.547	4003.823	28161.49	42360.48
188	188	2020-08-01	10742.4	2020	9.282	13039.85	5967.446	3945.268	28494.22	43099.15
189	189	2020-09-01	11570.1	2020	9.356	13035.03	5903.934	3881.991	28779.44	43769.27
190	190	2020-10-01	10786.0	2020	9.286	13035.99	5844.469	3822.203	29076.57	44460.51
191	191	2020-11-01	12664.4	2020	9.447	13041.04	5788.157	3765.252	29382.20	45167.96
192	192	2020-12-01	14438.4	2020	9.578	13041.51	5731.079	3708.545	29676.95	45861.90
193	193	2021-01-01	13329.9	2021	9.498	13039.67	5673.784	3652.267	29968.19	46555.48
194	194	2021-02-01	13265.0	2021	9.493	13039.28	5618.492	3598.031	30261.28	47254.39
195	195	2021-03-01	16787.5	2021	9.728	13040.67	5565.073	3545.644	30558.30	47962.84
196	196	2021-04-01	16204.3	2021	9.693	13040.09	5511.939	3494.084	30850.10	48666.21
197	197	2021-05-01	14234.8	2021	9.563	13040.48	5460.289	3444.080	31143.80	49375.79
198	198	2021-06-01	17218.5	2021	9.754	13041.57	5409.960	3395.499	31438.77	50090.58
199	199	2021-07-01	15263.1	2021	9.633	13037.69	5358.579	3346.830	31721.33	50788.74
200	200	2021-08-01	16678.9	2021	9.722	13041.88	5311.479	3301.383	32023.22	51521.03
201	201	2021-09-01	16234.1	2021	9.695	13040.67	5263.096	3255.662	32311.58	52234.83
202	202	2021-10-01	16293.6	2021	9.699	13040.91	5216.203	3211.373	32603.28	52957.20
203	203	2021-11-01	19328.2	2021	9.869	13042.18	5170.586	3168.357	32897.33	53686.65
204	204	2021-12-01	21352.0	2021	9.969	13042.30	5125.347	3126.045	33188.30	54414.31
205	205	2022-01-01	18211.1	2022	9.810	13041.84	5080.604	3084.464	33478.20	55143.94
206	206	2022-02-01	16638.5	2022	9.719	13041.74	5036.827	3043.922	33768.67	55877.55
207	207	2022-03-01	21962.4	2022	9.997	13042.09	4993.977	3004.365	34060.24	56616.31
208	208	2022-04-01	19757.4	2022	9.891	13041.94	4951.682	2965.556	34350.39	57355.92
209	209	2022-05-01	18609.3	2022	9.831	13042.04	4910.200	2927.634	34641.12	58099.76
210	210	2022-06-01	21003.9	2022	9.952	13042.31	4869.481	2890.553	34932.26	58847.55
211	211	2022-07-01	21345.0	2022	9.969	13041.34	4828.975	2853.974	35219.99	59592.85
212	212	2022-08-01	22150.6	2022	10.006	13042.39	4789.880	2818.593	35513.21	60350.69
213	213	2022-09-01	19808.1	2022	9.894	13042.09	4750.923	2783.643	35802.73	61105.54

Peramalan (1 tahun ke depan)

imporr.ts <- ts(data1$Impor, frequency=12,start=2005)
model.impor<-Arima(imporr.ts,order=c(1,1,0),seasonal=c(1,0,0))
ramalan_arima.impor = forecast(model.impor, 12)
ramalan_arima.impor%>%knitr::kable(digits = 3, 
                                caption = "Hasil peramalan 1 tahun ke depan ARIMA(1,1,0)x(1,0,0)12")

Hasil peramalan 1 tahun ke depan ARIMA(1,1,0)x(1,0,0)12
	Point Forecast	Lo 80	Hi 80	Lo 95	Hi 95
Oct 2022	20753.96	19131.39	22376.53	18272.45	23235.47
Nov 2022	21880.95	20015.08	23746.82	19027.35	24734.55
Dec 2022	23067.72	20836.00	25299.43	19654.60	26480.83
Jan 2023	21417.50	18932.34	23902.66	17616.78	25218.22
Feb 2023	20660.26	17921.92	23398.60	16472.33	24848.19
Mar 2023	23318.23	20357.61	26278.85	18790.35	27846.11
Apr 2023	22217.63	19046.51	25388.75	17367.82	27067.44
May 2023	21638.88	18271.93	25005.83	16489.57	26788.18
Jun 2023	22841.72	19289.09	26394.35	17408.44	28274.99
Jul 2023	23012.41	19283.60	26741.23	17309.68	28715.14
Aug 2023	23416.91	19519.76	27314.06	17456.73	29377.09
Sep 2023	22241.24	18182.78	26299.71	16034.35	28448.13

plot(ramalan_arima.impor)

Perbandingan DMA, DES, dan SARIMA

Untuk mengetahui metode manakah yang paling baik digunakan pada data ini, maka perbandingan nilai MAPE testing dari hasil pemulusan DMA, DES, dan Sarima yang terkecil menunjukkan bahwa metode tersebut adalah metode yang terbaik.

baris <- c("DMA (n = 4)", "DES (α = 0.6, β = 0.1)", "SARIMA")
MAPE <- c(MAPE.dma_test, MAPEtestingdes,0.8471548)
dmavsdes <- data.frame(baris, MAPE)
kableExtra::kable(dmavsdes, align = c(rep('c', times = 3)), 
                  col.names = c(" ", "MAPE"), digits = 3)

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

	MAPE
DMA (n = 4)	63.622
DES (α = 0.6, β = 0.1)	15.738
SARIMA	0.847

Berdasarkan hasil perbandingan MAPE DMA, DES, dan SARIMA (ARIMA(1,1,0)x(1,0,0)12) di dapatkan bahwa metode terbaik adalah SARIMA dengan nilai MAPEnya yang paling kecil.

Daftar Pustaka

Aniendyah, M. W. (2020). Peramalan hasil ekspor dan impor menggunakan metode Seasonal Autoregressive Integrated Moving Average (SARIMA) (Doctoral dissertation, Universitas Negeri Malang).

Bangun, R. H. B. (2016). Penerapan autoregressive integrated moving average (ARIMA) pada peramalan produksi kedelai di Sumatera Utara. Jurnal Agrica, 9(2), 90-100. 10.31289/agrica.v9i2.484

Benny, J. (2013). Ekspor dan impor pengaruhnya terhadap posisi cadangan devisa di Indonesia. Jurnal EMBA: Jurnal Riset Ekonomi, Manajemen, Bisnis Dan Akuntansi, 1(4). https://doi.org/10.35794/emba.1.4.2013.2920

Farina, F., & Husaini, A., (2015). Pengaruh dampak perkembangan tingkat ekspor dan impor terhadap nilai tukar negara ASEAN per dollar Amerika Serikat. Jurnal Administasi Bisnis, 50(6).

Nasirudin, F., Pindianti, M., Said, D. I. S., & Widodo, E. (2022). Peramalan jumlah produksi kopi di Jawa Timur pada tahun 2020-2021 menggunakan metode seasonal autoregressive integrated moving average (SARIMA). AGRIUM: Jurnal Ilmu Pertanian, 25(1), 34-43.

Purnama, S. (2020). Peramalan nilai ekspor dan impor Indonesia menggunakan metode ARIMA Box-Jenkins (Doctoral dissertation, UNIMED). Sedyaningrum, M., & Nuzula, N. F. (2016). Pengaruh Jumlah Nilai Ekspor, Impor dan Pertumbuhan Ekonomi terhadap Nilai Tukar dan Daya Beli Masyarakat di Indonesia. Jurnal Administrasi Bisnis, 34(1).

Sukirno, S. (2013). Teori Pengantar Makroekonomi Edisi Ketiga. Jakarta: Rajawali Pers.