1 Library & Data

pkgs <- c("forecast", "tseries", "FinTS", "rugarch",
          "ggplot2", "dplyr", "lubridate", "nortest",
          "lmtest", "zoo", "xts", "gridExtra")

for (p in pkgs) {
  if (!requireNamespace(p, quietly = TRUE)) install.packages(p)
}

library(forecast)
library(tseries)
library(FinTS)
library(rugarch)
library(ggplot2)
library(dplyr)
library(lubridate)
library(nortest)
library(lmtest)
library(zoo)
library(xts)
library(gridExtra)
df <- read.csv("C:\\Users\\aspir\\Documents\\Ekolan\\Data_Historis_GOTO (1).csv",
               stringsAsFactors = FALSE)
df$Date  <- as.Date(df$Date)
df       <- df[order(df$Date), ]

close_ts <- df$Close
n        <- length(close_ts)

cat("Jumlah observasi :", n, "\n")
## Jumlah observasi : 481
cat("Periode          :", format(min(df$Date)), "hingga", format(max(df$Date)), "\n")
## Periode          : 2022-04-11 hingga 2024-04-05
cat("Statistik deskriptif harga Close:\n")
## Statistik deskriptif harga Close:
summary(close_ts)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    56.0    88.0   112.0   157.4   218.0   404.0
ggplot(df, aes(x = Date, y = Close)) +
  geom_line(color = "#2c7bb6", linewidth = 0.7) +
  labs(title = "Harga Penutupan Saham GOTO (GOTO.JK)",
       subtitle = paste("Periode:", format(min(df$Date)), "–", format(max(df$Date))),
       x = "Tanggal", y = "Harga (IDR)") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold"))


2 Langkah 1 — Uji Stasioneritas Mean (ADF)

Uji Augmented Dickey-Fuller (ADF) digunakan untuk memeriksa apakah deret waktu bersifat stasioner. Hipotesis nol (H₀): data tidak stasioner (terdapat unit root).

2.1 Level (Data Asli)

adf_level <- adf.test(close_ts)
print(adf_level)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  close_ts
## Dickey-Fuller = -2.2508, Lag order = 7, p-value = 0.472
## alternative hypothesis: stationary
cat("\n-- Keputusan --\n")
## 
## -- Keputusan --
if (adf_level$p.value < 0.05) {
  cat("H0 ditolak: Data STASIONER pada level (p =",
      round(adf_level$p.value, 4), ")\n")
} else {
  cat("H0 tidak ditolak: Data TIDAK STASIONER pada level (p =",
      round(adf_level$p.value, 4), ") → perlu differencing\n")
}
## H0 tidak ditolak: Data TIDAK STASIONER pada level (p = 0.472 ) → perlu differencing

2.2 First Difference (d = 1)

close_diff <- diff(close_ts, differences = 1)
adf_diff   <- adf.test(close_diff)
print(adf_diff)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  close_diff
## Dickey-Fuller = -7.1166, Lag order = 7, p-value = 0.01
## alternative hypothesis: stationary
cat("\n-- Keputusan --\n")
## 
## -- Keputusan --
if (adf_diff$p.value < 0.05) {
  cat("H0 ditolak: Data STASIONER setelah differencing 1 kali (p =",
      round(adf_diff$p.value, 4), ") → d = 1\n")
} else {
  cat("H0 tidak ditolak: Data masih TIDAK STASIONER (p =",
      round(adf_diff$p.value, 4), ") → coba d = 2\n")
}
## H0 ditolak: Data STASIONER setelah differencing 1 kali (p = 0.01 ) → d = 1
df_diff <- data.frame(
  Date  = df$Date[-1],
  Diff  = close_diff
)

ggplot(df_diff, aes(x = Date, y = Diff)) +
  geom_line(color = "#d7191c", linewidth = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  labs(title = "Return/Differencing Pertama Harga Penutupan GOTO",
       x = "Tanggal", y = "ΔClose (IDR)") +
  theme_minimal(base_size = 13)


3 Langkah 2 — Correlogram (ACF & PACF)

ACF (Autocorrelation Function) dan PACF (Partial Autocorrelation Function) digunakan untuk mengidentifikasi orde AR (p) dan MA (q) yang sesuai pada model ARIMA.

3.1 Data Asli

par(mfrow = c(1, 2))
acf(close_ts,  lag.max = 12, main = "ACF — Harga Close (Level)")
pacf(close_ts, lag.max = 12, main = "PACF — Harga Close (Level)")

par(mfrow = c(1, 1))

3.2 Setelah Differencing (d = 1)

par(mfrow = c(1, 2))
acf(close_diff,  lag.max = 12, main = "ACF — ΔClose (d=1)")
pacf(close_diff, lag.max = 12, main = "PACF — ΔClose (d=1)")

par(mfrow = c(1, 1))

Pembacaan: Berdasarkan pola ACF dan PACF pada data yang sudah di-differencing, kandidat model ARIMA diidentifikasi. Lag yang signifikan menunjukkan kandidat orde AR dan MA.


4 Langkah 3 — Estimasi Model ARIMA (Kandidat)

Beberapa kandidat model ARIMA diestimasi berdasarkan hasil correlogram. Parameter yang disertakan: dengan konstanta (c) dan tanpa konstanta. Model dengan MA lag ke-5 saja [5] membutuhkan parameter fixed.

y <- close_ts   # data asli; differencing ditangani oleh order d=1

# ---- Model dengan konstanta (include.constant = TRUE) ----
fit_c110 <- Arima(y, order = c(1, 1, 0), include.constant = TRUE)
fit_c011 <- Arima(y, order = c(0, 1, 1), include.constant = TRUE)
fit_c111 <- Arima(y, order = c(1, 1, 1), include.constant = TRUE)

# ARIMA(0,1,[5]) dengan konstanta: MA hanya pada lag 5
fit_c015 <- Arima(y, order = c(0, 1, 5),
                  include.constant = TRUE,
                  fixed = c(0, 0, 0, 0, NA, NA))

# ARIMA(1,1,[5]) dengan konstanta: AR1 + MA hanya pada lag 5
fit_c115 <- Arima(y, order = c(1, 1, 5),
                  include.constant = TRUE,
                  fixed = c(NA, 0, 0, 0, 0, NA, NA))

# ---- Model tanpa konstanta ----
fit_110  <- Arima(y, order = c(1, 1, 0), include.constant = FALSE)
fit_011  <- Arima(y, order = c(0, 1, 1), include.constant = FALSE)
fit_111  <- Arima(y, order = c(1, 1, 1), include.constant = FALSE)

fit_015  <- Arima(y, order = c(0, 1, 5),
                  include.constant = FALSE,
                  fixed = c(0, 0, 0, 0, NA))

fit_115  <- Arima(y, order = c(1, 1, 5),
                  include.constant = FALSE,
                  fixed = c(NA, 0, 0, 0, 0, NA))

5 Langkah 4 — Pemilihan Model ARIMA Terbaik (AIC)

model_list <- list(
  "c ARIMA(1,1,0)"   = fit_c110,
  "c ARIMA(0,1,1)"   = fit_c011,
  "c ARIMA(1,1,1)"   = fit_c111,
  "c ARIMA(0,1,[5])" = fit_c015,
  "c ARIMA(1,1,[5])" = fit_c115,
  "ARIMA(1,1,0)"     = fit_110,
  "ARIMA(0,1,1)"     = fit_011,
  "ARIMA(1,1,1)"     = fit_111,
  "ARIMA(0,1,[5])"   = fit_015,
  "ARIMA(1,1,[5])"   = fit_115
)

aic_df <- data.frame(
  Model  = names(model_list),
  AIC    = sapply(model_list, AIC),
  BIC    = sapply(model_list, BIC),
  LogLik = sapply(model_list, logLik)
)

aic_df <- aic_df[order(aic_df$AIC), ]
rownames(aic_df) <- NULL
aic_df$AIC_per_obs <- aic_df$AIC / n

print(aic_df)
##               Model      AIC      BIC    LogLik AIC_per_obs
## 1  c ARIMA(1,1,[5]) 3319.835 3336.530 -1655.918    6.901944
## 2    ARIMA(1,1,[5]) 3319.844 3332.366 -1656.922    6.901963
## 3    c ARIMA(1,1,0) 3322.908 3335.430 -1658.454    6.908333
## 4    c ARIMA(0,1,1) 3322.957 3335.478 -1658.478    6.908434
## 5      ARIMA(1,1,0) 3323.286 3331.634 -1659.643    6.909118
## 6      ARIMA(0,1,1) 3323.494 3331.842 -1659.747    6.909552
## 7    c ARIMA(1,1,1) 3324.862 3341.557 -1658.431    6.912395
## 8      ARIMA(1,1,1) 3325.283 3337.805 -1659.642    6.913271
## 9  c ARIMA(0,1,[5]) 3333.567 3346.088 -1663.783    6.930492
## 10   ARIMA(0,1,[5]) 3334.428 3342.776 -1665.214    6.932283
best_model_name <- aic_df$Model[1]
cat("\n✔ Model terbaik berdasarkan AIC:", best_model_name, "\n")
## 
## ✔ Model terbaik berdasarkan AIC: c ARIMA(1,1,[5])
cat("  AIC =", round(aic_df$AIC[1], 6), "\n")
##   AIC = 3319.835
ggplot(aic_df, aes(x = reorder(Model, AIC), y = AIC,
                   fill = (Model == best_model_name))) +
  geom_col(show.legend = FALSE) +
  scale_fill_manual(values = c("steelblue", "#e63946")) +
  coord_flip() +
  labs(title = "Perbandingan AIC Kandidat Model ARIMA",
       subtitle = paste("Model terbaik:", best_model_name),
       x = NULL, y = "AIC") +
  theme_minimal(base_size = 12)

# Gunakan model dengan AIC terkecil secara dinamis — tidak hardcode
best_fit <- model_list[[best_model_name]]

cat("=== Ringkasan Model Terbaik:", best_model_name, "===\n")
## === Ringkasan Model Terbaik: c ARIMA(1,1,[5]) ===
summary(best_fit)
## Series: y 
## ARIMA(1,1,5) with drift 
## 
## Coefficients:
##          ar1  ma1  ma2  ma3  ma4     ma5    drift
##       0.1797    0    0    0    0  0.0981  -0.6626
## s.e.  0.0449    0    0    0    0  0.0429   0.4650
## 
## sigma^2 = 58.43:  log likelihood = -1655.92
## AIC=3319.84   AICc=3319.92   BIC=3336.53
## 
## Training set error measures:
##                       ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.007967425 7.612315 4.740044 0.07831046 3.106406 0.9773287
##                      ACF1
## Training set 0.0004522709
coeftest(best_fit)
## 
## z test of coefficients:
## 
##        Estimate Std. Error z value Pr(>|z|)    
## ar1    0.179675   0.044917  4.0002 6.33e-05 ***
## ma5    0.098117   0.042902  2.2870  0.02219 *  
## drift -0.662595   0.464956 -1.4251  0.15414    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Simpan orde ARMA untuk digunakan pada spesifikasi mean model ARCH
p_arma <- best_fit$arma[1]   # AR order
q_arma <- best_fit$arma[2]   # MA order
cat("\nOrde ARMA yang akan digunakan di mean model ARCH:",
    "AR =", p_arma, "| MA =", q_arma, "\n")
## 
## Orde ARMA yang akan digunakan di mean model ARCH: AR = 1 | MA = 5

6 Langkah 5 — Uji Diagnostik Residual Model Terbaik

6.1 5a. Uji Normalitas

resid_arima <- residuals(best_fit)

# Shapiro-Wilk (untuk n ≤ 5000)
if (length(resid_arima) <= 5000) {
  sw_test <- shapiro.test(resid_arima)
  cat("Shapiro-Wilk Test:\n")
  print(sw_test)
  cat("Keputusan:", ifelse(sw_test$p.value < 0.05,
      "H0 ditolak — Residual TIDAK berdistribusi normal\n",
      "H0 tidak ditolak — Residual berdistribusi normal\n"))
}
## Shapiro-Wilk Test:
## 
##  Shapiro-Wilk normality test
## 
## data:  resid_arima
## W = 0.86077, p-value < 2.2e-16
## 
## Keputusan: H0 ditolak — Residual TIDAK berdistribusi normal
# Kolmogorov-Smirnov
ks_test <- ks.test(scale(resid_arima), "pnorm")
cat("\nKolmogorov-Smirnov Test:\n")
## 
## Kolmogorov-Smirnov Test:
print(ks_test)
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  scale(resid_arima)
## D = 0.11948, p-value = 2.17e-06
## alternative hypothesis: two-sided
# Jarque-Bera
jb_test <- jarque.bera.test(resid_arima)
cat("\nJarque-Bera Test:\n")
## 
## Jarque-Bera Test:
print(jb_test)
## 
##  Jarque Bera Test
## 
## data:  resid_arima
## X-squared = 1692, df = 2, p-value < 2.2e-16
cat("Keputusan:", ifelse(jb_test$p.value < 0.05,
    "H0 ditolak — Residual TIDAK normal\n",
    "H0 tidak ditolak — Residual normal\n"))
## Keputusan: H0 ditolak — Residual TIDAK normal
par(mfrow = c(1, 2))
hist(resid_arima, breaks = 30, probability = TRUE,
     main = "Histogram Residual", xlab = "Residual",
     col = "#cce5ff", border = "white")
curve(dnorm(x, mean(resid_arima), sd(resid_arima)),
      add = TRUE, col = "red", lwd = 2)
qqnorm(resid_arima, main = "Q-Q Plot Residual")
qqline(resid_arima, col = "red", lwd = 2)

par(mfrow = c(1, 1))

6.2 5b. Uji White Noise (Ljung-Box)

cat("=== Uji Ljung-Box (White Noise) ===\n\n")
## === Uji Ljung-Box (White Noise) ===
for (lag in c(1, 6, 12, 24)) {
  lb <- Box.test(resid_arima, lag = lag, type = "Ljung-Box", fitdf = 2)
  cat(sprintf("Lag %-3d: Q = %.4f, p-value = %.4f → %s\n",
              lag, lb$statistic, lb$p.value,
              ifelse(lb$p.value > 0.05,
                     "Tidak ada autokorelasi (White Noise ✓)",
                     "Terdapat autokorelasi (bukan White Noise ✗)")))
}
## Lag 1  : Q = 0.0001, p-value = NaN → NA
## Lag 6  : Q = 3.9172, p-value = 0.4173 → Tidak ada autokorelasi (White Noise ✓)
## Lag 12 : Q = 10.0066, p-value = 0.4399 → Tidak ada autokorelasi (White Noise ✓)
## Lag 24 : Q = 32.2665, p-value = 0.0730 → Tidak ada autokorelasi (White Noise ✓)
par(mfrow = c(1, 2))
acf(resid_arima,  lag.max = 30, main = "ACF Residual")
pacf(resid_arima, lag.max = 30, main = "PACF Residual")

par(mfrow = c(1, 1))

6.3 5c. Uji Heteroskedastisitas — Efek ARCH (ARCH-LM)

Uji ARCH-LM (Lagrange Multiplier) digunakan untuk mendeteksi adanya heteroskedastisitas kondisional (efek ARCH) pada residual model ARIMA.

cat("=== Uji ARCH-LM pada Residual ARIMA ===\n\n")
## === Uji ARCH-LM pada Residual ARIMA ===
for (lag in 1:12) {
  arch_lm <- ArchTest(resid_arima, lags = lag)
  cat(sprintf("Lag %-3d: Chi² = %.4f, p-value = %.4f → %s\n",
              lag, arch_lm$statistic, arch_lm$p.value,
              ifelse(arch_lm$p.value < 0.05,
                     "Terdapat efek ARCH ✗",
                     "Tidak ada efek ARCH ✓")))
}
## Lag 1  : Chi² = 29.5696, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 2  : Chi² = 30.7980, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 3  : Chi² = 42.5202, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 4  : Chi² = 44.1316, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 5  : Chi² = 44.0574, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 6  : Chi² = 44.7274, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 7  : Chi² = 46.3282, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 8  : Chi² = 52.0703, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 9  : Chi² = 78.4336, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 10 : Chi² = 78.6379, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 11 : Chi² = 85.7399, p-value = 0.0000 → Terdapat efek ARCH ✗
## Lag 12 : Chi² = 87.5614, p-value = 0.0000 → Terdapat efek ARCH ✗

7 Langkah 6 — Pemodelan ARCH

Karena terdapat efek ARCH pada residual, dilakukan pemodelan volatilitas menggunakan model ARCH. Pemilihan orde terbaik mempertimbangkan dua kriteria: (1) nilai AIC terkecil dan (2) signifikansi statistik semua parameter alpha. Model yang lolos keduanya diprioritaskan.

7.1 Estimasi & Perbandingan ARCH(1)–ARCH(5)

# Inisialisasi vektor penyimpan hasil
aic_arch <- rep(NA, 5)
sig_arch <- rep(FALSE, 5)   # TRUE jika semua alpha signifikan (p < 0.05)

# Mean model mengikuti orde ARMA dari model ARIMA terbaik
cat(sprintf("Mean model yang digunakan: ARMA(%d, %d)\n\n", p_arma, q_arma))
## Mean model yang digunakan: ARMA(1, 5)
for (q in 1:5) {
  spec <- ugarchspec(
    variance.model     = list(model = "sGARCH", garchOrder = c(q, 0)),
    mean.model         = list(armaOrder = c(p_arma, q_arma),
                               include.mean = TRUE),
    distribution.model = "norm"
  )

  tryCatch({
    fit_arch_q <- ugarchfit(spec = spec, data = close_ts,
                            solver = "hybrid", out.sample = 0)

    # AIC (per observasi, skala rugarch)
    ic        <- infocriteria(fit_arch_q)
    aic_arch[q] <- ic["Akaike", ]

    # Signifikansi parameter alpha1, ..., alphaq
    coef_q   <- coef(fit_arch_q)
    se_q     <- sqrt(diag(vcov(fit_arch_q)))
    pval_q   <- 2 * pnorm(abs(coef_q / se_q), lower.tail = FALSE)

    alpha_names <- paste0("alpha", 1:q)
    # Ambil hanya alpha yang benar-benar ada di koefisien
    alpha_names <- alpha_names[alpha_names %in% names(pval_q)]
    alpha_pval  <- pval_q[alpha_names]
    all_sig     <- all(alpha_pval < 0.05, na.rm = TRUE)
    sig_arch[q] <- all_sig

    cat(sprintf(
      "ARCH(%d): AIC = %.6f | Alpha signifikan: %s\n",
      q, aic_arch[q],
      ifelse(all_sig,
             paste0("Semua signifikan ✓ (p: ",
                    paste(round(alpha_pval, 4), collapse = ", "), ")"),
             paste0("Ada yang tidak signifikan ✗ (p: ",
                    paste(round(alpha_pval, 4), collapse = ", "), ")"))
    ))
  }, error = function(e) {
    cat(sprintf("ARCH(%d): Gagal diestimasi — %s\n", q, e$message))
    aic_arch[q] <<- NA
    sig_arch[q] <<- FALSE
  })
}
## ARCH(1): AIC = 6.656117 | Alpha signifikan: Semua signifikan ✓ (p: 0)
## ARCH(2): AIC = 6.650132 | Alpha signifikan: Semua signifikan ✓ (p: 0, NaN)
## ARCH(3): AIC = 6.491010 | Alpha signifikan: Ada yang tidak signifikan ✗ (p: 0, 0.2365, 0)
## ARCH(4): AIC = 6.492193 | Alpha signifikan: Ada yang tidak signifikan ✗ (p: 0, 0.2834, 0, 0.3218)
## ARCH(5): AIC = 6.458745 | Alpha signifikan: Ada yang tidak signifikan ✗ (p: 0, 0.2419, 0, 1, 0.0044)

7.2 Pemilihan ARCH Terbaik: AIC + Signifikansi

# Tabel ringkasan perbandingan seluruh orde
arch_crit_df <- data.frame(
  Orde          = paste0("ARCH(", 1:5, ")"),
  AIC           = round(aic_arch, 6),
  Semua_Alpha_Sig = ifelse(sig_arch, "Ya ✓", "Tidak ✗"),
  Keterangan    = ifelse(
    is.na(aic_arch), "Gagal estimasi",
    ifelse(sig_arch, "Memenuhi kedua kriteria", "Hanya lolos AIC")
  )
)
print(arch_crit_df)
##      Orde      AIC Semua_Alpha_Sig              Keterangan
## 1 ARCH(1) 6.656117            Ya ✓ Memenuhi kedua kriteria
## 2 ARCH(2) 6.650132            Ya ✓ Memenuhi kedua kriteria
## 3 ARCH(3) 6.491010         Tidak ✗         Hanya lolos AIC
## 4 ARCH(4) 6.492193         Tidak ✗         Hanya lolos AIC
## 5 ARCH(5) 6.458745         Tidak ✗         Hanya lolos AIC
# Pilih: AIC terkecil di antara model yang semua alpha-nya signifikan
valid_idx <- which(!is.na(aic_arch) & sig_arch)

if (length(valid_idx) > 0) {
  best_q <- valid_idx[which.min(aic_arch[valid_idx])]
  cat(sprintf(
    "\n✔ Model ARCH terbaik (AIC terkecil & semua alpha signifikan): ARCH(%d)\n",
    best_q))
  cat("  AIC =", round(aic_arch[best_q], 6), "\n")
  cat("  Catatan: model ini memenuhi kriteria AIC sekaligus signifikansi parameter.\n")
} else {
  # Fallback: pilih berdasarkan AIC saja dan beri peringatan
  best_q <- which.min(aic_arch)
  cat(sprintf(
    "\n⚠ Tidak ada model dengan semua alpha signifikan.\n  Dipilih berdasarkan AIC terkecil saja: ARCH(%d)\n",
    best_q))
  cat("  AIC =", round(aic_arch[best_q], 6), "\n")
}
## 
## ✔ Model ARCH terbaik (AIC terkecil & semua alpha signifikan): ARCH(2)
##   AIC = 6.650132 
##   Catatan: model ini memenuhi kriteria AIC sekaligus signifikansi parameter.
# Warnai: merah = terpilih, biru = valid tapi bukan terbaik,
#         abu = tidak memenuhi signifikansi
arch_df <- data.frame(
  Orde   = paste0("ARCH(", 1:5, ")"),
  AIC    = aic_arch,
  Status = ifelse(
    is.na(aic_arch), "Gagal",
    ifelse(1:5 == best_q, "Terpilih",
           ifelse(sig_arch, "Memenuhi", "Tidak Signifikan"))
  )
)

ggplot(arch_df[!is.na(arch_df$AIC), ],
       aes(x = Orde, y = AIC, fill = Status)) +
  geom_col(show.legend = TRUE) +
  scale_fill_manual(
    values = c("Terpilih"       = "#e63946",
               "Memenuhi"       = "#2a9d8f",
               "Tidak Signifikan" = "steelblue",
               "Gagal"          = "gray70")) +
  labs(title    = "Perbandingan AIC Model ARCH(1)–ARCH(5)",
       subtitle = paste0("Model terpilih: ARCH(", best_q,
                         ") — AIC terkecil & semua alpha signifikan"),
       x = NULL, y = "AIC (per obs)", fill = "Status") +
  theme_minimal(base_size = 12)

7.3 Estimasi Model ARCH Terbaik

# Spesifikasi menggunakan orde terpilih (best_q) secara dinamis
spec_best <- ugarchspec(
  variance.model     = list(model = "sGARCH", garchOrder = c(best_q, 0)),
  mean.model         = list(armaOrder = c(p_arma, q_arma),
                             include.mean = TRUE),
  distribution.model = "norm"
)

fit_best_arch <- ugarchfit(spec = spec_best, data = close_ts,
                           solver = "hybrid")

cat(sprintf("=== Model Terpilih: %s-ARCH(%d) ===\n\n",
            best_model_name, best_q))
## === Model Terpilih: c ARIMA(1,1,[5])-ARCH(2) ===
show(fit_best_arch)
## 
## *---------------------------------*
## *          GARCH Model Fit        *
## *---------------------------------*
## 
## Conditional Variance Dynamics    
## -----------------------------------
## GARCH Model  : sGARCH(2,0)
## Mean Model   : ARFIMA(1,0,5)
## Distribution : norm 
## 
## Optimal Parameters
## ------------------------------------
##          Estimate  Std. Error    t value Pr(>|t|)
## mu     373.590626   18.874356   19.79356 0.000000
## ar1      1.000000    0.000975 1025.23198 0.000000
## ma1      0.281670    0.041367    6.80901 0.000000
## ma2      0.181271    0.141182    1.28395 0.199159
## ma3     -0.085357    0.018641   -4.57888 0.000005
## ma4      0.011037    0.028132    0.39232 0.694820
## ma5      0.018698    0.080655    0.23183 0.816673
## omega   16.078988   10.761944    1.49406 0.135160
## alpha1   0.487211    0.055800    8.73143 0.000000
## alpha2   0.511788    0.402208    1.27245 0.203215
## 
## Robust Standard Errors:
##          Estimate  Std. Error    t value Pr(>|t|)
## mu     373.590626  248.908404   1.500916  0.13338
## ar1      1.000000    0.005933 168.558712  0.00000
## ma1      0.281670    0.432208   0.651699  0.51460
## ma2      0.181271    1.837221   0.098666  0.92140
## ma3     -0.085357    0.637731  -0.133845  0.89353
## ma4      0.011037    0.583793   0.018906  0.98492
## ma5      0.018698    1.116192   0.016752  0.98664
## omega   16.078988  134.895126   0.119196  0.90512
## alpha1   0.487211    0.905995   0.537764  0.59074
## alpha2   0.511788    5.054367   0.101257  0.91935
## 
## LogLikelihood : -1589.357 
## 
## Information Criteria
## ------------------------------------
##                    
## Akaike       6.6501
## Bayes        6.7369
## Shibata      6.6493
## Hannan-Quinn 6.6843
## 
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
##                          statistic   p-value
## Lag[1]                       3.377 6.611e-02
## Lag[2*(p+q)+(p+q)-1][17]    12.928 2.306e-09
## Lag[4*(p+q)+(p+q)-1][29]    18.022 1.477e-01
## d.o.f=6
## H0 : No serial correlation
## 
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
##                         statistic p-value
## Lag[1]                    0.01556  0.9007
## Lag[2*(p+q)+(p+q)-1][5]   2.83118  0.4387
## Lag[4*(p+q)+(p+q)-1][9]   4.47556  0.5102
## d.o.f=2
## 
## Weighted ARCH LM Tests
## ------------------------------------
##             Statistic Shape Scale P-Value
## ARCH Lag[3]     2.703 0.500 2.000  0.1002
## ARCH Lag[5]     2.775 1.440 1.667  0.3241
## ARCH Lag[7]     3.274 2.315 1.543  0.4634
## 
## Nyblom stability test
## ------------------------------------
## Joint Statistic:  4.8799
## Individual Statistics:               
## mu     0.005941
## ar1    0.121712
## ma1    0.506718
## ma2    0.059753
## ma3    0.243392
## ma4    0.309561
## ma5    0.197394
## omega  2.905915
## alpha1 1.206375
## alpha2 0.560335
## 
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic:          2.29 2.54 3.05
## Individual Statistic:     0.35 0.47 0.75
## 
## Sign Bias Test
## ------------------------------------
##                    t-value   prob sig
## Sign Bias           1.2711 0.2043    
## Negative Sign Bias  0.6080 0.5435    
## Positive Sign Bias  0.9334 0.3511    
## Joint Effect        1.8301 0.6084    
## 
## 
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
##   group statistic p-value(g-1)
## 1    20     72.22    3.914e-08
## 2    30     85.88    1.546e-07
## 3    40    105.36    5.132e-08
## 4    50    103.93    7.968e-06
## 
## 
## Elapsed time : 0.5622752

7.4 Uji Signifikansi Parameter

coef_arch <- coef(fit_best_arch)
se_arch   <- sqrt(diag(vcov(fit_best_arch)))
z_arch    <- coef_arch / se_arch
p_arch    <- 2 * pnorm(abs(z_arch), lower.tail = FALSE)

sig_df <- data.frame(
  Parameter  = names(coef_arch),
  Estimasi   = round(coef_arch, 6),
  Std_Error  = round(se_arch,   6),
  z_value    = round(z_arch,    4),
  p_value    = round(p_arch,    4),
  Signifikan = ifelse(p_arch < 0.05, "Ya ✓", "Tidak ✗")
)

cat(sprintf("=== Signifikansi Parameter ARCH(%d) ===\n\n", best_q))
## === Signifikansi Parameter ARCH(2) ===
print(sig_df)
##        Parameter   Estimasi Std_Error   z_value p_value Signifikan
## mu            mu 373.590626       NaN       NaN     NaN       <NA>
## ar1          ar1   1.000000  0.000975 1025.2320       0       Ya ✓
## ma1          ma1   0.281670  0.041367    6.8090       0       Ya ✓
## ma2          ma2   0.181271       NaN       NaN     NaN       <NA>
## ma3          ma3  -0.085357       NaN       NaN     NaN       <NA>
## ma4          ma4   0.011037       NaN       NaN     NaN       <NA>
## ma5          ma5   0.018698       NaN       NaN     NaN       <NA>
## omega      omega  16.078988       NaN       NaN     NaN       <NA>
## alpha1    alpha1   0.487211  0.055800    8.7314       0       Ya ✓
## alpha2    alpha2   0.511788       NaN       NaN     NaN       <NA>
n_sig    <- sum(p_arch < 0.05)
n_total  <- length(p_arch)
alpha_idx <- grep("^alpha", names(p_arch))
n_alpha_sig <- sum(p_arch[alpha_idx] < 0.05)

cat(sprintf(
  "\nRingkasan: %d dari %d parameter signifikan pada α = 5%%.\n",
  n_sig, n_total))
## 
## Ringkasan: NA dari 10 parameter signifikan pada α = 5%.
cat(sprintf(
  "Parameter alpha (ARCH): %d dari %d signifikan.\n",
  n_alpha_sig, length(alpha_idx)))
## Parameter alpha (ARCH): NA dari 2 signifikan.

7.5 Uji ARCH-LM pada Residual Model Terbaik

resid_arch <- residuals(fit_best_arch, standardize = TRUE)

cat(sprintf("=== Uji ARCH-LM pada Residual %s-ARCH(%d) ===\n\n",
            best_model_name, best_q))
## === Uji ARCH-LM pada Residual c ARIMA(1,1,[5])-ARCH(2) ===
for (lag in c(1, 5, 10, 15)) {
  arch_lm2 <- ArchTest(resid_arch, lags = lag)
  cat(sprintf("Lag %-3d: Chi² = %.4f, p-value = %.4f → %s\n",
              lag, arch_lm2$statistic, arch_lm2$p.value,
              ifelse(arch_lm2$p.value > 0.05,
                     "Tidak ada efek ARCH (model cukup ✓)",
                     "Masih terdapat efek ARCH ✗")))
}
## Lag 1  : Chi² = 0.0154, p-value = 0.9011 → Tidak ada efek ARCH (model cukup ✓)
## Lag 5  : Chi² = 4.1271, p-value = 0.5313 → Tidak ada efek ARCH (model cukup ✓)
## Lag 10 : Chi² = 17.8188, p-value = 0.0581 → Tidak ada efek ARCH (model cukup ✓)
## Lag 15 : Chi² = 21.4241, p-value = 0.1238 → Tidak ada efek ARCH (model cukup ✓)
sigma_arch <- sigma(fit_best_arch)

vol_df <- data.frame(
  Date  = df$Date,
  Close = close_ts,
  Sigma = sigma_arch
)

p1 <- ggplot(vol_df, aes(x = Date, y = Close)) +
  geom_line(color = "#2c7bb6") +
  labs(title = "Harga Penutupan GOTO", x = NULL, y = "Harga (IDR)") +
  theme_minimal(base_size = 11)

p2 <- ggplot(vol_df, aes(x = Date, y = Sigma)) +
  geom_line(color = "#e63946") +
  labs(title = paste0("Estimasi Volatilitas Kondisional — ARCH(", best_q, ")"),
       x = "Tanggal", y = "σ_t (Std Dev Kondisional)") +
  theme_minimal(base_size = 11)

grid.arrange(p1, p2, ncol = 1)


8 Langkah 7 — Forecasting

Dilakukan peramalan harga saham GOTO menggunakan model gabungan terbaik untuk 30 hari ke depan.

8.1 Forecast Mean (ARIMA)

h <- 30  # horizon peramalan (hari kerja)

fc_arima <- forecast(best_fit, h = h, level = c(80, 95))
autoplot(fc_arima) +
  labs(title   = paste0("Forecast Harga GOTO — ", best_model_name),
       subtitle = paste("Horizon:", h, "hari ke depan"),
       x = "Waktu", y = "Harga (IDR)") +
  theme_minimal(base_size = 12)

tanggal_fc <- seq(max(df$Date) + 1, by = "day", length.out = h * 2)
tanggal_fc <- tanggal_fc[!weekdays(tanggal_fc) %in% c("Saturday", "Sunday")][1:h]

fc_df <- data.frame(
  Tanggal  = tanggal_fc,
  Forecast = round(as.numeric(fc_arima$mean),        2),
  Lower_80 = round(as.numeric(fc_arima$lower[, 1]),  2),
  Upper_80 = round(as.numeric(fc_arima$upper[, 1]),  2),
  Lower_95 = round(as.numeric(fc_arima$lower[, 2]),  2),
  Upper_95 = round(as.numeric(fc_arima$upper[, 2]),  2)
)

print(fc_df)
##       Tanggal Forecast Lower_80 Upper_80 Lower_95 Upper_95
## 1  2024-04-08    67.17    57.37    76.96    52.19    82.15
## 2  2024-04-09    66.74    51.59    81.89    43.57    89.91
## 3  2024-04-10    65.94    46.69    85.18    36.50    95.37
## 4  2024-04-11    65.44    42.80    88.09    30.81   100.08
## 5  2024-04-12    64.71    39.11    90.31    25.56   103.87
## 6  2024-04-15    64.04    35.37    92.71    20.20   107.88
## 7  2024-04-16    63.38    31.87    94.88    15.19   111.56
## 8  2024-04-17    62.71    28.59    96.84    10.52   114.90
## 9  2024-04-18    62.05    25.49    98.61     6.14   117.96
## 10 2024-04-19    61.39    22.55   100.23     1.99   120.79
## 11 2024-04-22    60.72    19.73   101.72    -1.97   123.42
## 12 2024-04-23    60.06    17.02   103.10    -5.76   125.89
## 13 2024-04-24    59.40    14.41   104.39    -9.41   128.21
## 14 2024-04-25    58.74    11.87   105.60   -12.94   130.41
## 15 2024-04-26    58.07     9.41   106.74   -16.35   132.50
## 16 2024-04-29    57.41     7.01   107.81   -19.67   134.49
## 17 2024-04-30    56.75     4.67   108.83   -22.90   136.40
## 18 2024-05-01    56.09     2.38   109.79   -26.05   138.22
## 19 2024-05-02    55.42     0.14   110.71   -29.13   139.97
## 20 2024-05-03    54.76    -2.06   111.58   -32.13   141.66
## 21 2024-05-06    54.10    -4.21   112.41   -35.08   143.28
## 22 2024-05-07    53.44    -6.33   113.20   -37.97   144.84
## 23 2024-05-08    52.77    -8.42   113.96   -40.81   146.35
## 24 2024-05-09    52.11   -10.47   114.69   -43.60   147.82
## 25 2024-05-10    51.45   -12.49   115.39   -46.34   149.23
## 26 2024-05-13    50.79   -14.48   116.05   -49.04   150.61
## 27 2024-05-14    50.12   -16.45   116.70   -51.69   151.94
## 28 2024-05-15    49.46   -18.39   117.31   -54.31   153.23
## 29 2024-05-16    48.80   -20.31   117.91   -56.89   154.49
## 30 2024-05-17    48.13   -22.21   118.48   -59.44   155.71

8.2 Forecast Volatilitas (ARCH)

fc_arch  <- ugarchforecast(fit_best_arch, n.ahead = h)

sigma_fc <- as.numeric(sigma(fc_arch))
mean_fc  <- as.numeric(fitted(fc_arch))

fc_vol_df <- data.frame(
  Hari     = 1:h,
  Tanggal  = tanggal_fc,
  Mean_FC  = round(mean_fc,                    4),
  Sigma_FC = round(sigma_fc,                   4),
  Lower_95 = round(mean_fc - 1.96 * sigma_fc,  4),
  Upper_95 = round(mean_fc + 1.96 * sigma_fc,  4)
)

print(fc_vol_df)
##    Hari    Tanggal Mean_FC Sigma_FC Lower_95 Upper_95
## 1     1 2024-04-08 68.1168   4.1302  60.0216  76.2120
## 2     2 2024-04-09 67.8786   4.9762  58.1252  77.6320
## 3     3 2024-04-10 67.9273   6.0724  56.0254  79.8293
## 4     4 2024-04-11 67.9386   6.8351  54.5418  81.3353
## 5     5 2024-04-12 67.9226   7.5969  53.0327  82.8124
## 6     6 2024-04-15 67.9226   8.2527  51.7473  84.0978
## 7     7 2024-04-16 67.9226   8.8768  50.5240  85.3211
## 8     8 2024-04-17 67.9226   9.4513  49.3981  86.4471
## 9     9 2024-04-18 67.9226   9.9964  48.3297  87.5155
## 10   10 2024-04-19 67.9226  10.5110  47.3210  88.5241
## 11   11 2024-04-22 67.9226  11.0022  46.3583  89.4869
## 12   12 2024-04-23 67.9226  11.4716  45.4382  90.4069
## 13   13 2024-04-24 67.9226  11.9225  44.5545  91.2907
## 14   14 2024-04-25 67.9226  12.3566  43.7037  92.1414
## 15   15 2024-04-26 67.9226  12.7756  42.8823  92.9628
## 16   16 2024-04-29 67.9226  13.1811  42.0875  93.7576
## 17   17 2024-04-30 67.9226  13.5743  41.3170  94.5281
## 18   18 2024-05-01 67.9226  13.9561  40.5686  95.2765
## 19   19 2024-05-02 67.9225  14.3275  39.8407  96.0044
## 20   20 2024-05-03 67.9225  14.6893  39.1316  96.7135
## 21   21 2024-05-06 67.9225  15.0421  38.4400  97.4051
## 22   22 2024-05-07 67.9225  15.3866  37.7647  98.0803
## 23   23 2024-05-08 67.9225  15.7234  37.1047  98.7404
## 24   24 2024-05-09 67.9225  16.0529  36.4589  99.3862
## 25   25 2024-05-10 67.9225  16.3755  35.8265 100.0185
## 26   26 2024-05-13 67.9225  16.6917  35.2068 100.6383
## 27   27 2024-05-14 67.9225  17.0018  34.5989 101.2461
## 28   28 2024-05-15 67.9225  17.3062  34.0024 101.8427
## 29   29 2024-05-16 67.9225  17.6051  33.4165 102.4285
## 30   30 2024-05-17 67.9225  17.8988  32.8408 103.0042
hist_df <- data.frame(
  Tanggal = df$Date,
  Close   = close_ts,
  Type    = "Aktual"
)

fc_plot_df <- data.frame(
  Tanggal  = tanggal_fc,
  Mean     = mean_fc,
  Lower_95 = mean_fc - 1.96 * sigma_fc,
  Upper_95 = mean_fc + 1.96 * sigma_fc
)

hist_recent <- tail(hist_df, 60)

ggplot() +
  geom_line(data = hist_recent,
            aes(x = Tanggal, y = Close),
            color = "#2c7bb6", linewidth = 0.8) +
  geom_ribbon(data = fc_plot_df,
              aes(x = Tanggal, ymin = Lower_95, ymax = Upper_95),
              fill = "#e63946", alpha = 0.2) +
  geom_line(data = fc_plot_df,
            aes(x = Tanggal, y = Mean),
            color = "#e63946", linewidth = 0.9, linetype = "dashed") +
  geom_vline(xintercept = as.numeric(max(df$Date)),
             linetype = "dotted", color = "gray50") +
  annotate("text",
           x = max(df$Date) + 2,
           y = max(hist_recent$Close),
           label = "← Forecast", hjust = 0,
           color = "#e63946", size = 3.5) +
  labs(title    = paste0("Forecast Harga Penutupan GOTO — ",
                         best_model_name, "-ARCH(", best_q, ")"),
       subtitle = paste("Horizon:", h, "hari kerja | Pita: 95% CI"),
       x = "Tanggal", y = "Harga (IDR)") +
  theme_minimal(base_size = 12)


9 Ringkasan Hasil Analisis

cat("============================================================\n")
## ============================================================
cat("  RINGKASAN ANALISIS DERET WAKTU SAHAM GOTO\n")
##   RINGKASAN ANALISIS DERET WAKTU SAHAM GOTO
cat("============================================================\n\n")
## ============================================================
cat("[1] Stasioneritas:\n")
## [1] Stasioneritas:
cat("    • Data level tidak stasioner → differencing 1 kali (d = 1)\n")
##     • Data level tidak stasioner → differencing 1 kali (d = 1)
cat("    • ADF setelah differencing: p =",
    round(adf_diff$p.value, 4), "→ stasioner\n\n")
##     • ADF setelah differencing: p = 0.01 → stasioner
cat("[2] Model ARIMA Terbaik:", best_model_name, "\n")
## [2] Model ARIMA Terbaik: c ARIMA(1,1,[5])
cat("    • Dipilih berdasarkan AIC terkecil di antara 10 kandidat model\n")
##     • Dipilih berdasarkan AIC terkecil di antara 10 kandidat model
cat("    • AIC =", round(aic_df$AIC[1], 4), "\n\n")
##     • AIC = 3319.835
cat("[3] Diagnostik Residual ARIMA:\n")
## [3] Diagnostik Residual ARIMA:
cat("    • Uji White Noise (Ljung-Box) : lihat output Bagian 5b\n")
##     • Uji White Noise (Ljung-Box) : lihat output Bagian 5b
cat("    • Uji Normalitas (JB/SW)      : lihat output Bagian 5a\n")
##     • Uji Normalitas (JB/SW)      : lihat output Bagian 5a
cat("    • Efek ARCH terdeteksi        : perlu pemodelan volatilitas\n\n")
##     • Efek ARCH terdeteksi        : perlu pemodelan volatilitas
cat("[4] Pemilihan Model ARCH:\n")
## [4] Pemilihan Model ARCH:
cat("    • Kandidat: ARCH(1) hingga ARCH(5)\n")
##     • Kandidat: ARCH(1) hingga ARCH(5)
cat("    • Kriteria: AIC terkecil + semua parameter alpha signifikan\n")
##     • Kriteria: AIC terkecil + semua parameter alpha signifikan
valid_str <- if (length(valid_idx) > 0) {
  paste0("ARCH(", valid_idx, ")", collapse = ", ")
} else "tidak ada"
cat("    • Orde yang memenuhi kedua kriteria:", valid_str, "\n")
##     • Orde yang memenuhi kedua kriteria: ARCH(1), ARCH(2)
cat(sprintf("    • Model terpilih: ARCH(%d) — AIC = %.6f\n\n",
            best_q, aic_arch[best_q]))
##     • Model terpilih: ARCH(2) — AIC = 6.650132
cat("[5] Diagnostik Residual ARCH:\n")
## [5] Diagnostik Residual ARCH:
cat(sprintf("    • Uji ARCH-LM pada residual standar ARCH(%d): lihat Bagian 6\n",
            best_q))
##     • Uji ARCH-LM pada residual standar ARCH(2): lihat Bagian 6
cat("    • Jika semua lag p > 0.05 → model volatilitas sudah memadai\n\n")
##     • Jika semua lag p > 0.05 → model volatilitas sudah memadai
cat("[6] Forecast:", h, "hari ke depan\n")
## [6] Forecast: 30 hari ke depan
cat(sprintf("    • Model gabungan: %s-ARCH(%d)\n", best_model_name, best_q))
##     • Model gabungan: c ARIMA(1,1,[5])-ARCH(2)
cat("    • Interval kepercayaan 95% disertakan\n")
##     • Interval kepercayaan 95% disertakan
cat("============================================================\n")
## ============================================================

Analisis dilakukan dengan R. Paket utama: forecast, tseries, rugarch, FinTS.