library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(tseries)
## Warning: package 'tseries' was built under R version 4.3.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.3.3
library(stats)

set.seed(123)

Parameter SARIMA(0,0,0)(0,1,2)^12 p=0, d=0, q=0 (non-seasonal) P=0, D=1, Q=2, s=12 (seasonal)

Simulate SARIMA model Karena D=1, kita perlu mensimulasi dengan seasonal differencing

n <- 120  # 10 tahun data bulanan
seasonal_period <- 12

# Parameter untuk seasonal MA(2)
seasonal_ma_coef <- c(0.6, -0.3)  # theta1 dan theta2

# Generate seasonal MA(2) process
seasonal_noise <- arima.sim(model = list(ma = seasonal_ma_coef), n = n)

# Apply seasonal integration (karena D=1)
seasonal_random_walk <- cumsum(seasonal_noise)

# Time series object
ts_data <- ts(seasonal_random_walk, frequency = 12, start = c(2014, 1))

ANALISIS DATA ORIGINAL (TANPA DIFFERENCING)

# 1. Time Series Plot - Original Data
par(mfrow = c(2, 2))

plot(ts_data, main = "Time Series Plot - Data Original", 
     ylab = "Nilai", xlab = "Waktu", col = "blue", lwd = 2)
grid()

# 2. ACF - Original Data
acf(ts_data, main = "ACF - Data Original", lag.max = 36)

# 3. PACF - Original Data
pacf(ts_data, main = "PACF - Data Original", lag.max = 36)

# 4. Decomposition
decompose_data <- decompose(ts_data)
plot(decompose_data$trend, main = "Trend Component", col = "red", lwd = 2)

ANALISIS DENGAN SEASONAL DIFFERENCING

# Seasonal differencing (lag-12 differencing)
seasonal_diff <- diff(ts_data, lag = 12)

# Remove NA values
seasonal_diff <- na.omit(seasonal_diff)

par(mfrow = c(2, 2))

# Time Series Plot - Setelah Seasonal Differencing
plot(seasonal_diff, main = "Time Series Plot - Setelah Seasonal Differencing", 
     ylab = "Differenced Values", xlab = "Waktu", col = "green", lwd = 2)
grid()
abline(h = 0, col = "red", lty = 2)

# ACF - Setelah Seasonal Differencing
acf(seasonal_diff, main = "ACF - Setelah Seasonal Differencing", lag.max = 36)

# PACF - Setelah Seasonal Differencing
pacf(seasonal_diff, main = "PACF - Setelah Seasonal Differencing", lag.max = 36)

# Histogram untuk melihat distribusi
hist(seasonal_diff, main = "Histogram - Data Setelah Seasonal Differencing", 
     col = "lightblue", breaks = 20)

ESTIMASI MODEL SARIMA

# Fit SARIMA(0,0,0)(0,1,2)^12 model
sarima_model <- Arima(ts_data, order = c(0,0,0), seasonal = c(0,1,2))
print(summary(sarima_model))
## Series: ts_data 
## ARIMA(0,0,0)(0,1,2)[12] 
## 
## Coefficients:
##          sma1    sma2
##       -1.1128  0.9999
## s.e.   0.1573  0.2587
## 
## sigma^2 = 8.167:  log likelihood = -288.24
## AIC=582.48   AICc=582.71   BIC=590.53
## 
## Training set error measures:
##                     ME     RMSE      MAE       MPE     MAPE    MASE      ACF1
## Training set 0.6732497 2.685856 2.100646 -187.1066 273.1632 0.62397 0.8664228

DIAGNOSTIK MODEL

# 1. Residual plot
plot(sarima_model$residuals, main = "Residuals", 
     ylab = "Residuals", col = "red")
abline(h = 0, col = "blue", lty = 2)

# 2. ACF of residuals
acf(sarima_model$residuals, main = "ACF of Residuals", lag.max = 24)

# 3. PACF of residuals
pacf(sarima_model$residuals, main = "PACF of Residuals", lag.max = 24)

# 4. Q-Q plot of residuals
qqnorm(sarima_model$residuals, main = "Q-Q Plot of Residuals")
qqline(sarima_model$residuals, col = "red")

Ljung-Box test untuk residuals

ljung_box <- Box.test(sarima_model$residuals, lag = 20, type = "Ljung-Box")
print(ljung_box)
## 
##  Box-Ljung test
## 
## data:  sarima_model$residuals
## X-squared = 382.21, df = 20, p-value < 2.2e-16

Shapiro-Wilk test untuk normalitas residuals

shapiro_test <- shapiro.test(sarima_model$residuals)
print(shapiro_test)
## 
##  Shapiro-Wilk normality test
## 
## data:  sarima_model$residuals
## W = 0.99298, p-value = 0.81

FORECASTING

# Forecast untuk 12 periode ke depan
forecast_result <- forecast(sarima_model, h = 12)
print(forecast_result)
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2024      12.030537 8.022837 16.03824 5.901289 18.15979
## Feb 2024      12.609187 8.601487 16.61689 6.479939 18.73844
## Mar 2024      12.379128 8.371428 16.38683 6.249880 18.50838
## Apr 2024      12.044756 8.037056 16.05246 5.915507 18.17400
## May 2024      13.100002 9.092302 17.10770 6.970754 19.22925
## Jun 2024      12.486345 8.478645 16.49404 6.357096 18.61559
## Jul 2024      11.070602 7.062902 15.07830 4.941353 17.19985
## Aug 2024       9.754269 5.746569 13.76197 3.625020 15.88352
## Sep 2024       9.655127 5.647427 13.66283 3.525879 15.78438
## Oct 2024       9.726312 5.718612 13.73401 3.597063 15.85556
## Nov 2024       9.743469 5.735769 13.75117 3.614221 15.87272
## Dec 2024      12.618369 8.610669 16.62607 6.489121 18.74762
# Plot forecast
par(mfrow = c(1, 1))
plot(forecast_result, main = "SARIMA(0,0,0)(0,1,2)^12 Forecast", 
     ylab = "Nilai", xlab = "Waktu")
grid()

SUMMARY STATISTIK

cat("Data Original:\n")
## Data Original:
cat("Mean:", mean(ts_data), "\n")
## Mean: 4.991389
cat("Variance:", var(ts_data), "\n")
## Variance: 10.70972
cat("Standard Deviation:", sd(ts_data), "\n")
## Standard Deviation: 3.272571
cat("\nData Setelah Seasonal Differencing:\n")
## 
## Data Setelah Seasonal Differencing:
cat("Mean:", mean(seasonal_diff), "\n")
## Mean: 0.2289271
cat("Variance:", var(seasonal_diff), "\n")
## Variance: 17.05882
cat("Standard Deviation:", sd(seasonal_diff), "\n")
## Standard Deviation: 4.130233
cat("\nModel Information:\n")
## 
## Model Information:
cat("AIC:", sarima_model$aic, "\n")
## AIC: 582.4838
cat("BIC:", BIC(sarima_model), "\n")
## BIC: 590.5302
cat("Log-likelihood:", sarima_model$loglik, "\n")
## Log-likelihood: -288.2419

INTERPRETASI ACF DAN PACF

cat("SARIMA(0,0,0)(0,1,2)^12 berarti:\n")
## SARIMA(0,0,0)(0,1,2)^12 berarti:
cat("- Tidak ada komponen ARIMA non-seasonal\n")
## - Tidak ada komponen ARIMA non-seasonal
cat("- Seasonal differencing order = 1 (D=1)\n")
## - Seasonal differencing order = 1 (D=1)
cat("- Seasonal MA order = 2 (Q=2)\n")
## - Seasonal MA order = 2 (Q=2)
cat("- Seasonal period = 12 (bulanan)\n\n")
## - Seasonal period = 12 (bulanan)
cat("Ekspektasi untuk ACF dan PACF setelah seasonal differencing:\n")
## Ekspektasi untuk ACF dan PACF setelah seasonal differencing:
cat("- ACF: Pola cut-off setelah lag 2 pada seasonal lags (12, 24)\n")
## - ACF: Pola cut-off setelah lag 2 pada seasonal lags (12, 24)
cat("- PACF: Pola tailing off pada seasonal lags\n")
## - PACF: Pola tailing off pada seasonal lags
cat("- Model ini cocok untuk data dengan pola seasonal yang kuat\n")
## - Model ini cocok untuk data dengan pola seasonal yang kuat