Import Library

rm(list=ls())
set.seed(42)
library(tseries)

## Warning: package 'tseries' was built under R version 4.5.3

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

library(forecast)

## Warning: package 'forecast' was built under R version 4.5.3

library(TSA)

## Warning: package 'TSA' was built under R version 4.5.3

## 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

Generate Data ARIMA(2,2,2)

n <- 200
phi1 <- 0.5
phi2 <- -0.7
theta1 <- -0.3
theta2 <- -0.5
sd_e <- sqrt(1.16)
e <- rnorm(n, mean = 0, sd = sd_e)
x <- numeric(n)

for (t in 3:n) {
  x[t] <- phi1 * x[t-1] +
          phi2 * x[t-2] +
          e[t] +
          theta1 * e[t-1] +
          theta2 * e[t-2]
}

x_d1 <- cumsum(x)        
data_sim <- cumsum(x_d1)

Hapus 50 Data AWal

data_trim <- data_sim[51:200] 

Splitting Data (80% Train, 20% Test)

data_trim <- data_sim[51:200]

n <- length(data_trim)
train_size <- floor(0.8 * n)

train <- ts(data_trim[1:train_size])
test  <- ts(data_trim[(train_size+1):n],
            start = train_size + 1)

ts.plot(train, main="Data Training")

Uji ADF

cat("ADF Data Asli:\n")

## ADF Data Asli:

print(adf.test(train))

## Warning in adf.test(train): p-value greater than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train
## Dickey-Fuller = 0.82052, Lag order = 4, p-value = 0.99
## alternative hypothesis: stationary

# Differencing ke-1
train_diff1 <- diff(train)
cat("\nADF Differencing 1:\n")

## 
## ADF Differencing 1:

print(adf.test(train_diff1))

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_diff1
## Dickey-Fuller = -4.1669, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

# Differencing ke-2
train_diff2 <- diff(train_diff1)
cat("\nADF Differencing 2:\n")

## 
## ADF Differencing 2:

print(adf.test(train_diff2))

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  train_diff2
## Dickey-Fuller = -8.2592, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary

ts.plot(train_diff1, main="Data Setelah Differencing 1")

ts.plot(train_diff2, main="Data Setelah Differencing 2")

ACF, PACF, EACF

par(mfrow=c(1,2))
acf(train_diff2, main="ACF")
pacf(train_diff2, main="PACF")

par(mfrow=c(1,1))
eacf(as.matrix(train_diff2))

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

Kandidat Model

model <- forecast::Arima(train, order=c(2,2,2))
model$coef

##        ar1        ar2        ma1        ma2 
##  0.4836874 -0.6429071 -0.3262104 -0.4903318

Model yang diperoleh yaitu

\[ ∇^2Y_t​=ϕ_1​∇^2Y_{t−1}​+ϕ_2​∇^2Y_{t−2}​+ε_{t}​+θ_1​ε_{t−1}​+θ_2​ε_{t−2} \\∇^2Y_t​=0.48​∇^2Y_{t−1}​-0.64​∇^2Y_{t−2}​+ε_{t}​-0.33​ε_{t−1}​-0.49​ε_{t−2} \]

AIC, BIC, AICc

calc_AICc <- function(model, n) {
  k <- length(model$coef)
  aic <- AIC(model)
  return(aic + (2*k*(k+1))/(n-k-1))
}

n_train <- length(train)

result <- data.frame(
  Model = c("ARIMA(2,2,2)"),
  AIC   = c(AIC(model)),
  BIC   = c(BIC(model)),
  AICc  = c(calc_AICc(model, n_train))
)
     
print(result)

##          Model      AIC     BIC     AICc
## 1 ARIMA(2,2,2) 342.4406 356.294 342.7884

Berdasarkan hasil estimasi, model yang digunakan adalah ARIMA(2,2,2) dengan nilai AIC sebesar 342.44, BIC sebesar 356.29, dan AICc sebesar 342.79

Diagnostik Model

ts.plot(residuals(model), main="Residual")

acf(residuals(model), main="ACF Residual")

cat("\nLjung-Box Test:\n")

## 
## Ljung-Box Test:

print(Box.test(residuals(model), lag=20, type="Ljung-Box"))

## 
##  Box-Ljung test
## 
## data:  residuals(model)
## X-squared = 4.1712, df = 20, p-value = 0.9999

print(shapiro.test(residuals(model)))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model)
## W = 0.98913, p-value = 0.4598

hist(residuals(model), main="Histogram Residual", col="lightblue")

qqnorm(residuals(model))
qqline(residuals(model), col="red")