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
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))
plot(data_ts, main = "Plot Deret Waktu Nilai Impor", ylab = "Nilai Impor", xlab = "Tahun")
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")
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
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
data_diff <- diff(data_boxcox)
plot(data_diff, main = "Plot Deret Waktu Setelah Differencing", ylab = "Nilai Impor", xlab = "Tahun")
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
acf(data_diff, lag.max = 100, main = "ACF setelah Differencing")
pacf(data_diff, lag.max = 100, main = "PACF setelah Differencing")
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
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
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
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
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.
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
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_m3,
main = "Forecast Nilai Impor Menggunakan Model ARIMA(1,1,1)",
xlab = "Tahun",
ylab = "Nilai Impor")