Data

Data yang digunakan dalam tugas ini adalah data time series yang digenerate berdasarkan model ARIMA(2,2,2) dengan \(\phi_1=0.5\), \(\phi_2=-0.7\), \(\theta_1=-0.3\), \(\theta_2=-0.5\), dan \(e_t \sim \mathcal{N}(0, 1.01)\).

Sehingga model yang terbentuk adalah \[ \scriptstyle \begin{aligned} \phi_p(B)(1-B)^dY_t &= \mu+\theta_q(B)e_t \\ \phi_2(B)(1-B)^2Y_t &= \mu+\theta_2(B)e_t \\ (1-\phi_1B-\phi_2B^2)(1-2B+B^2)Y_t &= \mu+(1-\theta_1B-\theta_2B)e_t \\ (1-(2+\phi_1)B+(1+2\phi_1-\phi_2)B^2-(\phi_1-2\phi_2)B^3-\phi_2B^4)Y_t &= \mu+e_t-\theta_1e_{t-1}-\theta_2e_{t-2} \\ Y_t-(2+\phi_1)Y_{t-1}+(1+2\phi_1-\phi_2)Y_{t-2}-(\phi_1-2\phi_2)Y_{t-3}-\phi_2Y_{t-4} &= \mu+e_t-\theta_1e_{t-1}-\theta_2e_{t-2} \\ Y_t-2.5Y_{t-1}+2.7Y_{t-2}-1.9Y_{t-3}+0.7Y_{t-4} &= \mu+e_t+0.3e_{t-1}+0.5e_{t-2} \\ Y_t &= \mu+e_t+2.5Y_{t-1}-2.7Y_{t-2}+1.9Y_{t-3}-0.7Y_{t-4}+0.3e_{t-1}+0.5e_{t-2} \\ \end{aligned} \] Kemudian data digenerate berdasarkan model yang terbentuk. Kemudian 50 data pertama dibuang dan diambil data yang ke 51-200 untuk analisis.

set.seed(123)

model <- list(
  ar = c(0.5, -0.7),
  ma = c(-0.3, -0.5)
)

X <- arima.sim(
  model = model,
  n = 200,
  sd = sqrt(1.01)
)

plot(X, type="l", main="Simulasi ARIMA(2,2,2)", col="red")

# buang 50 awal
X_final <- X[51:200]

plot(X_final, type="l", main="ARIMA setelah burn-in", col="red")

Splitting Data

Sebelum dilakukan analisis lebih lanjut, dilakukan splitting data terlebih dahulu yaitu dengan membagi data menjadi dua bagian bagian pertama yaitu sebagai data training sebesar 80% dan bagian kedua sebagai data testing sebesar 20%.

# jumlah data
n <- length(X_final)

# proporsi
train_size <- floor(0.8 * n)

# split
train <- X_final[1:train_size]
test  <- X_final[(train_size+1):n]

train_ts <- ts(train, start = 1, frequency = 1)
test_ts  <- ts(test, start = length(train) + 1, frequency = 1)

Uji Stasioner

Uji stasioner dilakukan menggunakan uji Augmented Dickey-Fuller (ADF) dengan hipotesis sebagai berikut: \(H_0\) : Data tidak stasioner \ \(H_1\) : Data stasioner

library(tseries)

adf.test(train_ts)

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

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

Berdasarkan hasil uji ADF diperoleh nilai p-value = 0.01 < 0.05, maka tolak \(H_0\) artinya data stasioner.

Identifikasi Kandidat Model

Identifikasi kandidat model diperoleh berdasarkan nilai \(p\) dan \(q\) dimana nilai \(d=0\).

acf(train_ts, main="ACF Data Stasioner")

pacf(train_ts, main="PACF Data Stasioner")

Berdasarkan plot PACF model yang diperoleh adalah \(ARIMA(2,0,0)\).

library(TSA)

## 
## Attaching package: 'TSA'

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

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

eacf(train_ts)

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

Berdasarkan hasil EACF model yang diperoleh adalah \(ARIMA(2,0,0)\), \(ARIMA(3,0,2)\), dan \(ARIMA(4,0,2)\).

Oleh karena itu diperoleh beberapa kandidat model yaitu \(ARIMA(2,0,0)\), \(ARIMA(3,0,2)\), dan \(ARIMA(4,0,2)\).

# ARIMA(2,0,0)
arima200 <- arima(train_ts, order=c(2,0,0), include.mean = TRUE, method = "ML")
arima200

## 
## Call:
## arima(x = train_ts, order = c(2, 0, 0), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ar2  intercept
##       0.3910  -0.8121    -0.0178
## s.e.  0.0524   0.0507     0.0724
## 
## sigma^2 estimated as 1.253:  log likelihood = -184.89,  aic = 375.79

# ARIMA(3,0,2)
arima302 <- arima(train_ts, order=c(3,0,2), include.mean = TRUE, method = "ML")
arima302

## 
## Call:
## arima(x = train_ts, order = c(3, 0, 2), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ar2      ar3     ma1      ma2  intercept
##       0.0527  -0.4538  -0.3671  0.0556  -0.8234    -0.0138
## s.e.  0.1143   0.0854   0.1114  0.0785   0.0636     0.0123
## 
## sigma^2 estimated as 0.9007:  log likelihood = -166.35,  aic = 344.69

# ARIMA(4,0,2)
arima402 <- arima(train_ts, order=c(4,0,2), include.mean = TRUE, method = "ML")
arima402

## 
## Call:
## arima(x = train_ts, order = c(4, 0, 2), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ar2      ar3     ar4     ma1      ma2  intercept
##       0.1003  -0.3702  -0.3810  0.1388  0.0462  -0.8665    -0.0132
## s.e.  0.1091   0.1002   0.0995  0.1070  0.0579   0.0574     0.0113
## 
## sigma^2 estimated as 0.8878:  log likelihood = -165.5,  aic = 345.01

AICKandidatModel <- c(375.79, 344.69, 345.01)
KndidatModelARIMA <- c("ARIMA(2,0,0)", "ARIMA(3,0,2)", "ARIMA(4,0,2)")
compmodelARIMA <- cbind(KndidatModelARIMA, AICKandidatModel)
colnames(compmodelARIMA) <- c("Kandidat Model", "Nilai AIC")
compmodelARIMA <- as.data.frame(compmodelARIMA)
compmodelARIMA

##   Kandidat Model Nilai AIC
## 1   ARIMA(2,0,0)    375.79
## 2   ARIMA(3,0,2)    344.69
## 3   ARIMA(4,0,2)    345.01

Model terbaik adalah model dengan nilai AIC terkecil yaitu model \(ARIMA(3,0,2)\), maka diperoleh \(Y_t\) dari penjabaran operator backshift dari model \(ARIMA(3,0,2)\) adalah \[ \scriptstyle \begin{aligned} \phi_p(B)(1-B)^dY_t &= \mu+\theta_q(B)e_t \\ \phi_3(B)(1-B)^0Y_t &= \mu+\theta_2(B)e_t \\ (1-\phi_1B-\phi_2B^2-\phi_3B^3)(1)Y_t &= \mu+(1-\theta_1B-\theta_2B)e_t \\ Y_t-\phi_1Y_{t-1}-\phi_2Y_{t-2}-\phi_3Y_{t-3} &= \mu+e_t-\theta_1e_{t-1}-\theta_2e_{t-2}\\ Y_t &= \mu +e_t+\phi_1Y_{t-1}+\phi_2Y_{t-2}+\phi_3Y_{t-3}-\theta_1e_{t-1}-\theta_2e_{t-2}\\ \end{aligned} \]

Penduga parameter \(\hat{\mu}=-0.0138\), \(\hat{\phi_1}=0.0527\), \(\hat{\phi_2}=-0.4538\), \(\hat{\phi_3}=-0.3671\), \(\hat{\theta_1}=0.0556\), \(\hat{\theta_2}=-0.8234\) dengan nilai dugaan \(\sigma_e^2\) sebesar 0.9007. Sehingga model yang diperoleh adalah \[ \begin{aligned} Y_t &= -0.0138+0.0527Y_{t-1}-0.4538Y_{t-2}-0.3671Y_{t-3}-0.0556e_{t-1}+0.8234e_{t-2}+e_t\\ \end{aligned} \]

Diagnostik Model

ARIMA302diag <- stats::arima(train_ts, order = c(3,0,2), method = "ML")
forecast::checkresiduals(ARIMA302diag)

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

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,0,2) with non-zero mean
## Q* = 3.7742, df = 5, p-value = 0.5824
## 
## Model df: 5.   Total lags used: 10

Berdasarkan plot ACF dan PACF terlihat bahwa tidak ada lag yang signifikan. Hal tersebut menunjukkan bahwa tidak ada gejala autokorelasi pada sisaan.

sisaan <- arima302$residuals

# Uji formal normalitas data
ks.test(sisaan,"pnorm")

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  sisaan
## D = 0.078129, p-value = 0.4565
## alternative hypothesis: two-sided

# Uji nilai tengah sisaan
t.test(sisaan, mu = 0, alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  sisaan
## t = -0.070853, df = 119, p-value = 0.9436
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1784281  0.1660999
## sample estimates:
##    mean of x 
## -0.006164081

# Uji autokorelasi
Box.test(sisaan, lag = 23 ,type = "Ljung")

## 
##  Box-Ljung test
## 
## data:  sisaan
## X-squared = 18.827, df = 23, p-value = 0.7112