Analisis Deret Waktu A

Simulasi dan Analisis Model SARIMA (0,0,0)(2,1,2)[12]

1. Simulasi Data Seasonal ARIMA (0,0,0)(2,1,2)[12]

set.seed(123)

e <- rnorm(700, mean=0, sd=1)

phi1 <- 0.2
phi2 <- -0.15
theta1 <- 0.15
theta2 <- 0.1

y <- numeric(length(e))
for (t in 25:length(e)) {
  y[t] <- y[t-12] + phi1 * y[t-12*1] + phi2 * y[t-12*2] + e[t] + theta1 * e[t-12*1] + theta2 * e[t-12*2]
}
y <- ts(y, frequency=12)

autoplot(y) + ggtitle("Simulated SARIMA (0,0,0)(2,1,2)[12] Time Series") + ylab("Value") + xlab("Time")

2. Plot ACF dan PACF

par(mfrow=c(1,2))
acf(y, main="ACF of Simulated Data")
pacf(y, main="PACF of Simulated Data")

par(mfrow=c(1,1))

3. Seasonal Differencing (lag=12)

y_diff <- diff(y, lag=12)
autoplot(y_diff) + ggtitle("Seasonally Differenced Time Series") + ylab("Differenced Value") + xlab("Time")

par(mfrow=c(1,2))
acf(y_diff, main="ACF of Differenced Data")
pacf(y_diff, main="PACF of Differenced Data")

par(mfrow=c(1,1))

4. Fitting Model SARIMA

fit <- Arima(y, order=c(0,0,0), seasonal=c(2,1,2), method="ML")
summary(fit)

## Series: y 
## ARIMA(0,0,0)(2,1,2)[12] 
## 
## Coefficients:
##         sar1    sar2    sma1     sma2
##       0.1040  0.8721  0.2694  -0.6166
## s.e.  0.0576  0.0560  0.0632   0.0436
## 
## sigma^2 = 1.06:  log likelihood = -1007.69
## AIC=2025.39   AICc=2025.48   BIC=2048.06
## 
## Training set error measures:
##                      ME     RMSE       MAE     MPE     MAPE      MASE
## Training set 0.01587643 1.017843 0.8004089 19.5369 42.21261 0.5908429
##                      ACF1
## Training set -0.007262761

5. Diagnostic Check

checkresiduals(fit)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0)(2,1,2)[12]
## Q* = 26.213, df = 20, p-value = 0.1589
## 
## Model df: 4.   Total lags used: 24

6. Plot Normality Residual

resid_fit <- residuals(fit)
library(ggplot2)
ggplot(data.frame(resid=resid_fit), aes(x=resid)) +
  geom_histogram(aes(y=..density..), bins=30, fill="gray", color="black") +
  stat_function(fun=dnorm, args=list(mean=mean(resid_fit), sd=sd(resid_fit)), color="darkblue", size=1) +
  ggtitle("Histogram and Normal Density of Residuals") +
  xlab("Residuals") + ylab("Density")

shapiro.test(residuals(fit))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(fit)
## W = 0.99812, p-value = 0.6475