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)
n <- 200

# Parameter
phi1  <-  0.5
phi2  <- -0.7
teta1 <- -0.3
teta2 <- -0.5

ar    <- numeric(4)
ar[1] <-  phi1 + 2           
ar[2] <-  phi2 - 2*phi1 - 1  
ar[3] <-  phi1 - 2*phi2      
ar[4] <-  phi2              

ma    <- numeric(2)
ma[1] <- -teta1              
ma[2] <- -teta2             

et <- rnorm(n, mean = 0, sd = sqrt(1.01))
Yt <- numeric(n)

for (t in 5:n) {
  Yt[t] <- ar[1] * Yt[t-1] +
            ar[2] * Yt[t-2] +
            ar[3] * Yt[t-3] +
            ar[4] * Yt[t-4] +
            et[t]            +
            ma[1] * et[t-1]  +
            ma[2] * et[t-2]
}

# Buang 50 data pertama
data <- ts(Yt[-(1:50)])


# Plot data
plot(data, type="l", col="blue",
     main="Time Series Data",
     ylab="Y_t", xlab="Waktu")

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(data)

# proporsi
train_size <- floor(0.8 * n)

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

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)

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

Berdasarkan hasil uji ADF diperoleh nilai p-value = 0.4373 > 0.05, maka terima \(H_0\) artinya data tidak stasioner. Sehingga perlu dilakukan differencing.

Differencing

d1 <- diff(train, differences = 1)
adf.test(d1)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  d1
## Dickey-Fuller = -1.7984, Lag order = 4, p-value = 0.6603
## alternative hypothesis: stationary

Setelah dilakukan differencing orde 1, diperoleh nilai p-value = 0.66 > 0.05, maka data masih belum stasioner. Sehingga dilanjutkan differencing orde 2.

d2 <- diff(train, differences = 2)
adf.test(d2)

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

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

Setelah dilakukan differencing orde 2, diperoleh nilai p-value = 0.01 < 0.05, maka data sudah stasioner.

Identifikasi Kandidat Model

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

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

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

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

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

library(TSA)
eacf(d2)

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

knitr::include_graphics("C:/Users/User/OneDrive/Dokumen/SEMESTER 6/Topik Statistika I/plot.eacf.png")

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

Karena data yang digunakan sudah stasioner maka \(d=0\). Oleh karena itu dapat ditentukan beberapa kandidat model yaitu \(ARIMA(2,0,0)\), \(ARIMA(2,0,2)\), \(ARIMA(2,0,3)\), \(ARIMA(3,0,2)\), dan \(ARIMA(3,0,3)\).

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

## 
## Call:
## arima(x = d2, order = c(2, 0, 0), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ar2  intercept
##       0.7566  -0.5476    -0.0325
## s.e.  0.0777   0.0772     0.1151
## 
## sigma^2 estimated as 0.9712:  log likelihood = -166.21,  aic = 338.41

arima202 <- arima(d2, order=c(2,0,2), include.mean = TRUE, method = "ML")
arima202

## 
## Call:
## arima(x = d2, order = c(2, 0, 2), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ar2      ma1     ma2  intercept
##       0.6972  -0.7064  -0.0006  0.3732    -0.0287
## s.e.  0.1220   0.0893   0.1449  0.1329     0.1204
## 
## sigma^2 estimated as 0.9226:  log likelihood = -163.29,  aic = 336.58

arima203 <- arima(d2, order=c(2,0,3), include.mean = TRUE, method = "ML")
arima203

## 
## Call:
## arima(x = d2, order = c(2, 0, 3), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ar2     ma1     ma2     ma3  intercept
##       0.6866  -0.7534  0.0276  0.4222  0.0848    -0.0286
## s.e.  0.1209   0.1157  0.1572  0.1561  0.1582     0.1271
## 
## sigma^2 estimated as 0.9202:  log likelihood = -163.15,  aic = 338.29

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

## 
## Call:
## arima(x = d2, order = c(3, 0, 2), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##          ar1      ar2     ar3      ma1     ma2  intercept
##       0.8267  -0.8322  0.0977  -0.1198  0.4073    -0.0284
## s.e.  0.3237   0.2991  0.2266   0.3095  0.1522     0.1254
## 
## sigma^2 estimated as 0.9211:  log likelihood = -163.2,  aic = 338.4

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

## 
## Call:
## arima(x = d2, order = c(3, 0, 3), include.mean = TRUE, method = "ML")
## 
## Coefficients:
##           ar1      ar2      ar3     ma1     ma2     ma3  intercept
##       -0.1776  -0.1503  -0.6166  0.9280  0.4142  0.4324    -0.0313
## s.e.   0.1520   0.1390   0.0987  0.1688  0.2099  0.1245     0.1240
## 
## sigma^2 estimated as 0.8891:  log likelihood = -161.37,  aic = 336.73

AICKandidatModel <- c(353.73, 351.9, 353.61, 353.72, 352.05)
KndidatModelARIMA <- c("ARIMA(2,0,0)", "ARIMA(2,0,2)", "ARIMA(2,0,3)", "ARIMA(3,0,2)", "ARIMA(3,0,3)")
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)    353.73
## 2   ARIMA(2,0,2)     351.9
## 3   ARIMA(2,0,3)    353.61
## 4   ARIMA(3,0,2)    353.72
## 5   ARIMA(3,0,3)    352.05

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

Penduga parameter \(\hat{\mu}=-0.0306\), \(\hat{\phi_1}=0.6972\), \(\hat{\phi_2}=-0.7064\), \(\hat{\theta_1}=-0.0006\), \(\hat{\theta_2}=0.3732\) dengan nilai dugaan \(\sigma_e^2\) sebesar 1.05. Sehingga model yang diperoleh adalah \[ \begin{aligned} Y_t &= -0.0306+0.6972Y_{t-1}-0.7064Y_{t-2}+0.0006e_{t-1}-0.3732e_{t-2}+e_t\\ \end{aligned} \]

Diagnostik Model

ARIMA202diag <- stats::arima(d2, order = c(2,0,2), method = "ML")
forecast::checkresiduals(ARIMA202diag)

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

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,2) with non-zero mean
## Q* = 1.5853, df = 6, p-value = 0.9536
## 
## Model df: 4.   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 <- arima202$residuals

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

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

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

## 
##  One Sample t-test
## 
## data:  sisaan
## t = 0.032496, df = 117, p-value = 0.9741
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1729765  0.1787476
## sample estimates:
##  mean of x 
## 0.00288558

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

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

Kesimpulan

Sehingga dapat disimpulkan bahwa semua asumsi terpenuhi