1. Memuat library yang dibutuhkan

library(TSA)
## Warning: package 'TSA' was built under R version 4.4.2
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(tseries)
## Warning: package 'tseries' was built under R version 4.4.2
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(lmtest)
## Warning: package 'lmtest' was built under R version 4.4.2
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(forecast)
## Warning: package 'forecast' was built under R version 4.4.2
## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.4.2
## Warning: package 'readr' was built under R version 4.4.2
## Warning: package 'forcats' was built under R version 4.4.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ readr::spec()   masks TSA::spec()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(rmarkdown)
## Warning: package 'rmarkdown' was built under R version 4.4.2
library(readxl)
## Warning: package 'readxl' was built under R version 4.4.2

2. Memanggil data

Nilai_Impor <- read_excel("D:/Data R SEI/Data Nilai Impor Indonesia dari Januari 2015 sd November 2024.xlsx")
Data <- as.numeric(Nilai_Impor$`Nilai Impor Indonesia (Juta US$)`)
data_ts <- ts(Data, start = c(2015, 1), frequency = 12)
paged_table(as.data.frame(Nilai_Impor))

3. Plot deret waktu

plot(data_ts, main = "Plot Deret Waktu Nilai Impor", ylab = "Nilai Impor", xlab = "Tahun")

4. Uji stasioneritas

Transformasi Box-Cox

data_boxcox <- forecast::BoxCox(data_ts, lambda = forecast::BoxCox.lambda(data_ts))
plot(data_boxcox, main = "Deret Waktu Setelah Transformasi Box-Cox", ylab = "Nilai Impor (Box-Cox)", xlab = "Tahun")

Uji ADF setelah transformasi

adf_test_boxcox <- adf.test(data_boxcox, alternative = "stationary")
print(adf_test_boxcox)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_boxcox
## Dickey-Fuller = -2.2578, Lag order = 4, p-value = 0.4695
## alternative hypothesis: stationary

Plot ACF dan PACF

acf(data_boxcox, lag.max = 100, main = "ACF setelah Transformasi Box-Cox")

pacf(data_boxcox, lag.max = 100, main = "PACF setelah Transformasi Box-Cox")

karena p-value = 0.4695 > 0.05 menyatakan bahwa data tidak stasioner, maka langkah selanjutnya adalah melakukan differencing

Differencing

data_diff <- diff(data_boxcox)
plot(data_diff, main = "Plot Deret Waktu Setelah Differencing", ylab = "Nilai Impor", xlab = "Tahun")

Uji ADF setelah differencing

adf_test_diff <- adf.test(data_diff, alternative = "stationary")
## Warning in adf.test(data_diff, alternative = "stationary"): p-value smaller
## than printed p-value
print(adf_test_diff)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data_diff
## Dickey-Fuller = -4.7989, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

Plot ACF dan PACF setelah differencing

acf(data_diff, lag.max = 100, main = "ACF setelah Differencing")

pacf(data_diff, lag.max = 100, main = "PACF setelah Differencing")

5. Identifikasi model (menggunakan ARIMA manual)

6. Estimasi parameter

Model 1: ARIMA(1,1,0)(0,1,0)[12]

m1 <- arima(data_ts, order = c(1, 1, 0), seasonal = list(order = c(0, 1, 0), period = 12), method = "ML")
cat("Model 1:\n")
## Model 1:
print(lmtest::coeftest(m1))
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.582826   0.079124 -7.3659 1.759e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model 2: ARIMA(0,1,1)(0,1,0)[12]

m2 <- arima(data_ts, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 0), period = 12), method = "ML")
cat("\nModel 2:\n")
## 
## Model 2:
print(lmtest::coeftest(m2))
## 
## z test of coefficients:
## 
##      Estimate Std. Error z value  Pr(>|z|)    
## ma1 -0.504483   0.064997 -7.7616 8.387e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model 3: ARIMA(1,1,1)(0,1,0)[12]

m3 <- arima(data_ts, order = c(1, 1, 1), seasonal = list(order = c(0, 1, 0), period = 12), method = "ML")
cat("\nModel 3:\n")
## 
## Model 3:
print(lmtest::coeftest(m3))
## 
## z test of coefficients:
## 
##     Estimate Std. Error z value  Pr(>|z|)    
## ar1 -0.45007    0.12059 -3.7324 0.0001897 ***
## ma1 -0.21992    0.11477 -1.9162 0.0553357 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

7. Diagnostik model

cat("\nDiagnostic Checks for Model 1:\n")
## 
## Diagnostic Checks for Model 1:
checkresiduals(m1)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)(0,1,0)[12]
## Q* = 75.296, df = 23, p-value = 1.802e-07
## 
## Model df: 1.   Total lags used: 24
cat("\nDiagnostic Checks for Model 2:\n")
## 
## Diagnostic Checks for Model 2:
checkresiduals(m2)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,1,1)(0,1,0)[12]
## Q* = 68.875, df = 23, p-value = 1.812e-06
## 
## Model df: 1.   Total lags used: 24
cat("\nDiagnostic Checks for Model 3:\n")
## 
## Diagnostic Checks for Model 3:
checkresiduals(m3)

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

Membandingkan AIC

cat("\nComparing AIC values:\n")
## 
## Comparing AIC values:
aic.model <- data.frame(Nama = c("m1", "m2", "m3"),
                        Model = c("ARIMA(1,1,0)", "ARIMA(0,1,1)", "ARIMA(1,1,1)"),
                        AIC = c(m1$aic, m2$aic, m3$aic))
print(aic.model)
##   Nama        Model      AIC
## 1   m1 ARIMA(1,1,0) 1893.403
## 2   m2 ARIMA(0,1,1) 1899.782
## 3   m3 ARIMA(1,1,1) 1892.084

Untuk peramalan satu periode ke depan, Model 3 (ARIMA(1,1,1)) adalah pilihan terbaik berdasarkan AIC, karena AIC terendah menunjukkan model yang lebih efisien dalam memprediksi dengan kualitas fit yang lebih baik.

Bangun model ARIMA

m3 <- auto.arima(data_ts)
summary(m3)
## Series: data_ts 
## ARIMA(0,1,2)(0,0,2)[12] 
## 
## Coefficients:
##           ma1     ma2    sma1    sma2
##       -0.7040  0.3271  0.5325  0.2169
## s.e.   0.0928  0.1029  0.1010  0.1187
## 
## sigma^2 = 2139101:  log likelihood = -1027.51
## AIC=2065.02   AICc=2065.56   BIC=2078.87
## 
## Training set error measures:
##                   ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 71.4623 1431.511 1099.949 -0.1894583 7.453007 0.4542068
##                      ACF1
## Training set -0.002239148

8. Forecast

forecast_horizon <- 1  # Desember 2024
forecast_m3 <- forecast::forecast(m3, h = forecast_horizon)
print(forecast_m3)
##          Point Forecast    Lo 80    Hi 80   Lo 95    Hi 95
## Dec 2024       20620.78 18746.43 22495.14 17754.2 23487.36

Plot forecast

plot(forecast_m3, 
     main = "Forecast Nilai Impor Menggunakan Model ARIMA(1,1,1)", 
     xlab = "Tahun", 
     ylab = "Nilai Impor")