########################################
### EDVINA KRESNANINGRUM 
### NRP 5003221053
### Model Deret Waktu untuk Data Musiman (SARIMA) 
### ARIMA(0,0,0) (2,0,2)^12)
########################################

library(tseries)
## Warning: package 'tseries' was built under R version 4.3.3
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(TSA)
## Warning: package 'TSA' was built under R version 4.3.3
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
library(forecast)
## Warning: package 'forecast' was built under R version 4.3.3
## Registered S3 methods overwritten by 'forecast':
##   method       from
##   fitted.Arima TSA 
##   plot.Arima   TSA
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
##Menentukan Parameter dari generate data
#1. Simulasi data musiman ARIMA (0,0,0)(2,0,2)^12
set.seed(123) #agar hasil konsisten
datarandom <- arima.sim(n = 100, model = list(order = c(0,0,0), 
                                        seasonal = list(order = c(2,0,2), period = 12)))

#2. simpan data
write.csv(datarandom, file = "datarandom_arima.csv", row.names = FALSE)
getwd()  #cek di folder mana file disimpan
## [1] "C:/Users/kresn/Downloads"
#3.panggil data 
data <- read.csv("C:/Users/kresn/OneDrive/Documents/datarandom_arima.csv")

#4. melihat estimasi parameter aktual dari datarandom
fit <- Arima(data, order = c(0,0,0), 
             seasonal = list(order = c(2,0,2), period = 12))
summary(fit)
## Series: data 
## ARIMA(0,0,0)(2,0,2)[12] with non-zero mean 
## 
## Coefficients:
##         sar1    sar2     sma1     sma2    mean
##       0.4421  0.4247  -0.7032  -0.2967  0.0821
## s.e.  0.4983  0.3889   0.5461   0.4911  0.0678
## 
## sigma^2 = 0.7919:  log likelihood = -129.96
## AIC=271.92   AICc=272.83   BIC=287.55
## 
## Training set error measures:
##                        ME      RMSE       MAE      MPE     MAPE      MASE
## Training set -0.004230376 0.8673351 0.6917656 131.7239 135.2947 0.6861125
##                     ACF1
## Training set -0.05351072
##Model SARIMA
#Parameter hasil estimasi
set.seed(123)
n <- 124  # agar cukup untuk lag 24 + 100 data utama
at <- rnorm(n, mean = 0, sd = sqrt(0.83))  #pakai sigma^2 dari output
phi1 <- -0.6551  #SAR(1)
phi2 <- -0.7266  #SAR(2)
theta1 <- 0.6552 #SMA(1)
theta2 <- 0.9999 #SMA(2)

#Inisialisasi output y dan loop simulasi
y <- numeric(n)
for (t in 25:n) {
  y[t] <- phi1 * y[t-12] +
    phi2 * y[t-24] +
    at[t] +
    theta1 * at[t-12] +
    theta2 * at[t-24]
}
y_final <- y[25:n]

## ANALISIS
#1. Time Series Plot
data_timeseries <- ts(y_final, frequency = 12)
plot.ts(data_timeseries, main = "Time Series Plot", ylab = "y", col = "blue")

#2. Plot ACF & PACF Data
acf(data_timeseries, main = "ACF")

pacf(data_timeseries, main = "PACF")