library(quantmod)
library(rugarch)
library(tseries)
library(FinTS)
library(ggplot2)
library(dplyr)
library(lubridate)
library(moments)
library(forecast)
library(gridExtra)
library(scales)

getSymbols("BRIS.JK",
           src  = "yahoo",
           from = "2022-01-01",
           to   = "2026-03-31")

harga   <- as.numeric(Cl(BRIS.JK))
tanggal <- index(Cl(BRIS.JK))

cat("Jumlah data harga  :", length(harga), "\n")

cat("Rentang tanggal    :", format(min(tanggal)), "–", format(max(tanggal)), "\n")

cat("Contoh harga (5 pertama):\n")

print(head(harga, 5))

# r_t = ln(P_t / P_{t-1}) x 100  (dalam persen)
log_return     <- diff(log(harga)) * 100
tanggal_return <- tanggal[-1]

cat("Jumlah observasi return:", length(log_return), "\n")

cat("Contoh return (5 pertama):\n")

print(head(log_return, 5))

df_harga  <- data.frame(Tanggal = tanggal,        Harga  = harga)
df_return <- data.frame(Tanggal = tanggal_return,  Return = log_return)

# Plot 1: Harga Penutupan
p1 <- ggplot(df_harga, aes(x = Tanggal, y = Harga)) +
  geom_line(color = "#1565C0", linewidth = 0.6) +
  scale_y_continuous(labels = comma) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  labs(title    = "Harga Penutupan Saham BRIS",
       subtitle = paste("Periode:", format(min(tanggal), "%d %b %Y"),
                        "–", format(max(tanggal), "%d %b %Y")),
       x = "Tanggal", y = "Harga (IDR)") +
  theme_bw(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))
print(p1)

# Plot 2: Log Return Harian
p2 <- ggplot(df_return, aes(x = Tanggal, y = Return)) +
  geom_line(color = "#C62828", linewidth = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  labs(title = "Log Return Harian Saham BRIS (%)",
       x = "Tanggal", y = "Log Return (%)") +
  theme_bw(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))
print(p2)

# Plot 3: Histogram Return
p3 <- ggplot(df_return, aes(x = Return)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 50, fill = "#1E88E5",
                 color = "white", alpha = 0.8) +
  stat_function(fun  = dnorm,
                args = list(mean = mean(log_return), sd = sd(log_return)),
                color = "red", linewidth = 1, linetype = "dashed") +
  labs(title    = "Histogram Log Return BRIS",
       subtitle = "Garis merah = distribusi normal",
       x = "Log Return (%)", y = "Densitas") +
  theme_bw(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))
print(p3)

skew_val <- skewness(log_return, na.rm = TRUE)
kurt_val <- kurtosis(log_return, na.rm = TRUE)

desc_stats <- data.frame(
  Statistik = c("Jumlah Obs.", "Mean (%)", "Minimum (%)",
                "Maksimum (%)", "Standar Deviasi (%)",
                "Skewness", "Kurtosis"),
  Nilai = c(
    length(log_return),
    round(mean(log_return, na.rm = TRUE), 6),
    round(min(log_return,  na.rm = TRUE), 6),
    round(max(log_return,  na.rm = TRUE), 6),
    round(sd(log_return,   na.rm = TRUE), 6),
    round(skew_val, 4),
    round(kurt_val, 4)
  )
)
print(desc_stats, row.names = FALSE)

cat("\nInterpretasi:\n")

if (abs(skew_val) > 0.5) {
  arah <- ifelse(skew_val > 0, "kanan (positif)", "kiri (negatif)")
  cat("  - Return menceng ke", arah, "→ distribusi tidak simetris\n")
} else {
  cat("  - Distribusi relatif simetris\n")
}

if (kurt_val > 3) {
  cat("  - Kurtosis > 3 → distribusi leptokurtik (ekor tebal / fat-tailed)\n")
  cat("    Cocok menggunakan distribusi Student-t dalam model GARCH\n")
} else {
  cat("  - Kurtosis mendekati 3 → distribusi mendekati normal\n")
}

jb_test <- jarque.bera.test(log_return)

cat("Uji Normalitas Jarque-Bera:\n")

cat("  Statistik:", round(jb_test$statistic, 4), "\n")

cat("  p-value  :", format(jb_test$p.value, scientific = TRUE, digits = 4), "\n")

if (jb_test$p.value < 0.05) {
  cat("  Kesimpulan: Data TIDAK berdistribusi normal\n")
} else {
  cat("  Kesimpulan: Data berdistribusi normal\n")
}

adf_harga  <- adf.test(harga)
adf_return <- adf.test(log_return)

cat("ADF Test — Harga Penutupan:\n")

cat("  Statistik:", round(adf_harga$statistic, 4),
    "| p-value:", round(adf_harga$p.value, 4), "\n")

cat("  Kesimpulan:",
    ifelse(adf_harga$p.value < 0.05, "STASIONER", "TIDAK STASIONER"), "\n\n")

cat("ADF Test — Log Return:\n")

cat("  Statistik:", round(adf_return$statistic, 4),
    "| p-value:", round(adf_return$p.value, 4), "\n")

cat("  Kesimpulan:",
    ifelse(adf_return$p.value < 0.05, "STASIONER", "TIDAK STASIONER"), "\n")

# Differencing otomatis jika return tidak stasioner
if (adf_return$p.value >= 0.05) {
  cat("\nReturn tidak stasioner → dilakukan differencing...\n")
  log_return     <- diff(log_return)
  tanggal_return <- tanggal_return[-1]
  adf_diff <- adf.test(log_return)
  cat("ADF setelah differencing:\n")
  cat("  Statistik:", round(adf_diff$statistic, 4),
      "| p-value:", round(adf_diff$p.value, 4), "\n")
}

# 7. Uji Efek ARCH (ARCH-LM)
arch_test <- ArchTest(log_return, lags = 12)

cat("ARCH-LM Test (lag = 12):\n")

cat("  Chi-squared:", round(arch_test$statistic, 4), "\n")

cat("  p-value    :", format(arch_test$p.value, scientific = TRUE, digits = 4), "\n")

if (arch_test$p.value < 0.05) {
  cat("  Kesimpulan: Terdapat efek ARCH → Pemodelan GARCH tepat digunakan\n")
} else {
  cat("  Kesimpulan: Tidak ada efek ARCH signifikan\n")
}

n       <- length(log_return)
n_train <- round(n * 0.80)
n_test  <- n - n_train

train_return  <- log_return[1:n_train]
test_return   <- log_return[(n_train + 1):n]
train_tanggal <- tanggal_return[1:n_train]
test_tanggal  <- tanggal_return[(n_train + 1):n]

cat("Total observasi return :", n, "\n")

cat("Data training (80%)    :", n_train, "obs. |",
    format(min(train_tanggal)), "–", format(max(train_tanggal)), "\n")

cat("Data testing  (20%)    :", n_test,  "obs. |",
    format(min(test_tanggal)),  "–", format(max(test_tanggal)), "\n")

# Fungsi bantu estimasi
estimasi_garch <- function(p, q, dist = "norm", data) {
  spec <- ugarchspec(
    variance.model     = list(model = "sGARCH", garchOrder = c(p, q)),
    mean.model         = list(armaOrder = c(0, 0), include.mean = TRUE),
    distribution.model = dist
  )
  fit <- tryCatch(
    ugarchfit(spec = spec, data = data, solver = "hybrid"),
    error = function(e) NULL
  )
  return(fit)
}

# Daftar model kandidat
model_config <- list(
  list(p = 1, q = 1, dist = "norm", nama = "GARCH(1,1)-Normal"),
  list(p = 1, q = 1, dist = "std",  nama = "GARCH(1,1)-Student-t"),
  list(p = 1, q = 2, dist = "norm", nama = "GARCH(1,2)-Normal"),
  list(p = 1, q = 2, dist = "std",  nama = "GARCH(1,2)-Student-t"),
  list(p = 2, q = 1, dist = "norm", nama = "GARCH(2,1)-Normal"),
  list(p = 2, q = 1, dist = "std",  nama = "GARCH(2,1)-Student-t")
)

hasil_ic   <- data.frame()
model_list <- list()

for (cfg in model_config) {
  cat("  Estimasi:", cfg$nama, "... ")
  fit <- estimasi_garch(cfg$p, cfg$q, cfg$dist, train_return)
  if (!is.null(fit)) {
    ic  <- infocriteria(fit)
    row <- data.frame(Model = cfg$nama,
                      AIC   = round(ic[1], 6),
                      BIC   = round(ic[2], 6),
                      p = cfg$p, q = cfg$q, dist = cfg$dist)
    hasil_ic           <- rbind(hasil_ic, row)
    model_list[[cfg$nama]] <- fit
    cat("OK\n")
  } else {
    cat("GAGAL\n")
  }
}

cat("\nPerbandingan AIC & BIC semua model:\n")

print(hasil_ic[, c("Model", "AIC", "BIC")], row.names = FALSE)

# Pilih model terbaik (AIC terkecil)
best_idx  <- which.min(hasil_ic$AIC)
best_nama <- hasil_ic$Model[best_idx]
best_fit  <- model_list[[best_nama]]
best_cfg  <- model_config[[best_idx]]   # simpan konfigurasi untuk rolling forecast

cat("\nModel Terbaik (AIC terkecil):", best_nama, "\n")

cat("  AIC =", hasil_ic$AIC[best_idx],
    "| BIC =", hasil_ic$BIC[best_idx], "\n")

cat("Ringkasan model terbaik:", best_nama, "\n\n")

coef_best <- coef(best_fit)
print(round(coef_best, 6))

cat("\nLog-Likelihood:", round(likelihood(best_fit), 4), "\n")

mu    <- ifelse("mu"    %in% names(coef_best), coef_best["mu"],    NA)
omega <- ifelse("omega" %in% names(coef_best), coef_best["omega"], NA)
alpha <- coef_best[grep("^alpha", names(coef_best))]
beta  <- coef_best[grep("^beta",  names(coef_best))]

cat("\nInterpretasi parameter:\n")

cat("  μ (Mean)  =", round(mu, 6),
    "→ rata-rata return harian\n")

cat("  ω (Omega) =", round(omega, 8),
    "→ komponen volatilitas jangka panjang\n")

cat("  α (Alpha) =", round(sum(alpha), 6),
    "→ pengaruh shock/guncangan terhadap volatilitas\n")

cat("  β (Beta)  =", round(sum(beta), 6),
    "→ persistensi volatilitas dari periode sebelumnya\n")

persistensi <- sum(alpha) + sum(beta)
cat("α + β =", round(persistensi, 6), "\n\n")

if (persistensi >= 0.95) {
  cat("  Sangat mendekati 1 → Volatilitas sangat PERSISTEN\n")
  cat("  Guncangan pada volatilitas akan bertahan sangat lama\n")
} else if (persistensi >= 0.90) {
  cat("  Persistensi TINGGI → Volatilitas lambat kembali ke rata-rata\n")
} else {
  cat("  Persistensi MODERAT → Volatilitas cukup cepat kembali ke rata-rata\n")
}

half_life <- log(0.5) / log(persistensi)
cat("\n  Half-life volatilitas ≈", round(half_life, 1), "hari perdagangan\n")

cat("  (waktu yang dibutuhkan shock untuk memudar setengahnya)\n")

std_resid <- residuals(best_fit, standardize = TRUE)

lb_resid   <- Box.test(std_resid,   lag = 10, type = "Ljung-Box")
lb_sq      <- Box.test(std_resid^2, lag = 10, type = "Ljung-Box")
arch_resid <- ArchTest(std_resid, lags = 12)

cat("Ljung-Box (residual standar) — lag 10:\n")

cat("  p-value:", round(lb_resid$p.value, 4),
    ifelse(lb_resid$p.value > 0.05,
           "→ Tidak ada autokorelasi ",
           "→ Masih ada autokorelasi "), "\n\n")

cat("Ljung-Box (kuadrat residual) — lag 10:\n")

cat("  p-value:", round(lb_sq$p.value, 4),
    ifelse(lb_sq$p.value > 0.05,
           "→ Tidak ada autokorelasi pada kuadrat ",
           "→ Masih ada efek ARCH pada residual ️"), "\n\n")

cat("ARCH-LM (residual standar) — lag 12:\n")

cat("  p-value:", round(arch_resid$p.value, 4),
    ifelse(arch_resid$p.value > 0.05,
           "→ Tidak ada efek ARCH tersisa ",
           "→ Masih ada efek ARCH ️"), "\n\n")

if (lb_resid$p.value > 0.05 &&
    lb_sq$p.value    > 0.05 &&
    arch_resid$p.value > 0.05) {
  cat("  Model VALID: residual bersih dari autokorelasi & efek ARCH\n")
} else {
  cat("  Pertimbangkan model alternatif atau spesifikasi lebih lanjut\n")
}

# --- Rolling 1-step forecast pada data testing ---
# Gunakan ARMA(0,0) konsisten dengan saat estimasi
best_spec_roll <- ugarchspec(
  variance.model     = list(model = "sGARCH",
                           garchOrder = c(best_cfg$p, best_cfg$q)),
  mean.model         = list(armaOrder = c(0, 0), include.mean = TRUE),
  distribution.model = best_cfg$dist
)

cat("Menjalankan rolling 1-step forecast pada data testing...\n")

roll_fc <- ugarchroll(
  spec             = best_spec_roll,
  data             = log_return,
  n.ahead          = 1,
  forecast.length  = n_test,
  refit.every      = 20,
  refit.window     = "recursive",
  solver           = "hybrid",
  calculate.VaR    = FALSE
)

roll_df    <- as.data.frame(roll_fc)
pred_vol   <- roll_df$Sigma
actual_ret <- roll_df$Realized

cat("Selesai! Jumlah forecast:", length(pred_vol), "\n")

# --- Forecast 10 hari ke depan (out-of-sample) ---
cat("\nForecast 10 hari ke depan (out-of-sample):\n")

fc_10 <- ugarchforecast(best_fit, n.ahead = 10)

vol_10hari <- as.numeric(sigma(fc_10))
ret_10hari <- as.numeric(fitted(fc_10))

# Buat tanggal hari kerja (hindari weekend)
tanggal_fc <- seq(max(tanggal_return) + 1, by = "day", length.out = 20)
tanggal_fc <- tanggal_fc[!weekdays(tanggal_fc) %in% 
                           c("Saturday", "Sunday")][1:10]

fc10_df <- data.frame(
  Tanggal     = tanggal_fc,
  Return      = round(ret_10hari, 6),
  Volatilitas = round(vol_10hari, 6),
  Lower_95    = round(ret_10hari - 1.96 * vol_10hari, 6),
  Upper_95    = round(ret_10hari + 1.96 * vol_10hari, 6)
)

cat("\nForecast 10 Hari + Interval Prediksi 95%:\n")

print(fc10_df, row.names = FALSE)

# --- Evaluasi Model ---
actual_vol <- abs(actual_ret)

rmse <- sqrt(mean((pred_vol - actual_vol)^2))
mae  <- mean(abs(pred_vol - actual_vol))

# MAPE berbasis harga (Nusrang et al., 2025)
pred_mean   <- roll_df$Mu
harga_awal  <- harga[n_train + 1]

harga_ramal <- numeric(n_test)
harga_ramal[1] <- harga_awal * exp(pred_mean[1] / 100)
for (i in 2:n_test) {
  harga_ramal[i] <- harga_ramal[i - 1] * exp(pred_mean[i] / 100)
}

harga_aktual <- harga[(n_train + 2):(n_train + n_test + 1)]
mape <- mean(abs((harga_aktual - harga_ramal) / harga_aktual)) * 100

cat("\nEvaluasi Model:\n")

cat("  RMSE :", round(rmse, 6), "(akurasi prediksi volatilitas)\n")

cat("  MAE  :", round(mae,  6), "(akurasi prediksi volatilitas)\n")

cat("  MAPE :", round(mape, 4), "% (akurasi rekonstruksi harga)\n")

# MAPE Volatilitas
# Proxy: sqrt(return^2) = |return|, tapi filter hari return = 0
idx_valid <- which(actual_vol > 0)

mape_vol <- mean(abs((actual_vol[idx_valid] - pred_vol[idx_valid]) / 
                       actual_vol[idx_valid])) * 100

cat("\nEvaluasi Model:\n")

cat("  RMSE :", round(rmse, 6), "\n")

cat("  MAE  :", round(mae,  6), "\n")

cat("  MAPE :", round(mape_vol, 4), "%\n")

# 15. Visualisasi Hasil
actual_vol <- abs(actual_ret)
# --- Plot 1: Volatilitas GARCH periode training ---
vol_train    <- as.numeric(sigma(best_fit))
df_vol_train <- data.frame(
  Tanggal     = train_tanggal,
  Volatilitas = vol_train
)

p_vol <- ggplot(df_vol_train, aes(x = Tanggal, y = Volatilitas)) +
  geom_line(color = "#E65100", linewidth = 0.6) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  labs(title    = paste("Volatilitas GARCH —", best_nama),
       subtitle = "Periode Training",
       x = "Tanggal", y = "Volatilitas Bersyarat (%)") +
  theme_bw(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))
print(p_vol)

# --- Plot 2: Return vs Volatilitas ---
df_rv <- data.frame(
  Tanggal     = train_tanggal,
  Return      = train_return,
  Volatilitas = vol_train
)

p_rv <- ggplot(df_rv, aes(x = Tanggal)) +
  geom_line(aes(y = Return,      color = "Return"),      linewidth = 0.4) +
  geom_line(aes(y = Volatilitas, color = "Volatilitas"), linewidth = 0.7) +
  scale_color_manual(values = c("Return"      = "#1565C0",
                                "Volatilitas" = "#B71C1C")) +
  scale_x_date(date_breaks = "1 year", date_labels = "%Y") +
  labs(title  = "Return vs Volatilitas GARCH — BRIS",
       x = "Tanggal", y = "(%)", color = "") +
  theme_bw(base_size = 12) +
  theme(plot.title      = element_text(face = "bold"),
        legend.position = "top")
print(p_rv)

# --- Plot 3: Prediksi vs Aktual (data testing) ---
df_test_fc <- data.frame(
  Tanggal  = test_tanggal,
  Aktual   = actual_vol,
  Prediksi = pred_vol
)

p_pred <- ggplot(df_test_fc, aes(x = Tanggal)) +
  geom_line(aes(y = Aktual,   color = "Aktual |r|"), linewidth = 0.5) +
  geom_line(aes(y = Prediksi, color = "Prediksi σ"), linewidth = 0.7) +
  scale_color_manual(values = c("Aktual |r|" = "#37474F",
                                "Prediksi σ" = "#E53935")) +
  scale_x_date(date_breaks = "3 months", date_labels = "%b %Y") +
  labs(title    = "Prediksi vs Aktual Volatilitas — Data Testing",
       subtitle = paste("Model:", best_nama),
       x = "Tanggal", y = "Volatilitas / |Return| (%)", color = "") +
  theme_bw(base_size = 12) +
  theme(plot.title      = element_text(face = "bold"),
        legend.position = "top",
        axis.text.x     = element_text(angle = 45, hjust = 1))
print(p_pred)

# --- Plot 4: Interval Prediksi ±2σ pada data testing ---
df_test_fc$Upper <- pred_vol * 2
df_test_fc$Lower <- -pred_vol * 2

p_interval <- ggplot(df_test_fc, aes(x = Tanggal)) +
  geom_ribbon(aes(ymin = Lower, ymax = Upper),
              fill = "#FFCDD2", alpha = 0.6) +
  geom_line(data = data.frame(Tanggal = test_tanggal, Return = actual_ret),
            aes(y = Return), color = "#1A237E", linewidth = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_date(date_breaks = "3 months", date_labels = "%b %Y") +
  labs(title    = "Interval Prediksi 95% (±2σ GARCH) — Data Testing",
       subtitle = "Area merah muda = interval prediksi | Biru = return aktual",
       x = "Tanggal", y = "Return (%)") +
  theme_bw(base_size = 12) +
  theme(plot.title  = element_text(face = "bold"),
        axis.text.x = element_text(angle = 45, hjust = 1))
print(p_interval)

# --- Plot 5: Forecast Return 10 Hari ke Depan ---
p_fc10 <- ggplot(fc10_df, aes(x = Tanggal)) +
  geom_ribbon(aes(ymin = Lower_95, ymax = Upper_95),
              fill = "#C8E6C9", alpha = 0.7) +
  geom_line(aes(y = Return), color = "#2E7D32",
            linewidth = 1) +
  geom_point(aes(y = Return), color = "#1B5E20", size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_date(date_breaks = "1 day", date_labels = "%d %b") +
  labs(title    = "Forecast Return 10 Hari ke Depan",
       subtitle = "Hijau = interval prediksi 95% | Garis = return prediksi",
       x = "Tanggal", y = "Return (%)") +
  theme_bw(base_size = 12) +
  theme(plot.title = element_text(face = "bold"))
print(p_fc10)