PERSIAPAN: INSTALASI DAN LOAD PAKET

required_packages <- c(
  "quantmod",   
  "ggplot2",    
  "dplyr",      
  "tidyr",      
  "quadprog",   
  "gridExtra",  
  "ggcorrplot", 
  "scales",     
  "moments",    
  "tseries"     
)

new_pkg <- required_packages[
  !(required_packages %in% installed.packages()[, "Package"])
]
if (length(new_pkg) > 0) {
  install.packages(new_pkg, repos = "https://cloud.r-project.org")
}

library(quantmod)
## Warning: package 'quantmod' was built under R version 4.5.1
## Loading required package: xts
## Warning: package 'xts' was built under R version 4.5.1
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.5.1
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: TTR
## Warning: package 'TTR' was built under R version 4.5.1
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(ggplot2)
library(dplyr)
## 
## ######################### Warning from 'xts' package ##########################
## #                                                                             #
## # The dplyr lag() function breaks how base R's lag() function is supposed to  #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or       #
## # source() into this session won't work correctly.                            #
## #                                                                             #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop           #
## # dplyr from breaking base R's lag() function.                                #
## #                                                                             #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning.  #
## #                                                                             #
## ###############################################################################
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:xts':
## 
##     first, last
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(quadprog)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library(ggcorrplot)
## Warning: package 'ggcorrplot' was built under R version 4.5.3
library(scales)
library(moments)
library(tseries)
## Warning: package 'tseries' was built under R version 4.5.1

KONFIGURASI PARAMETER PENELITIAN

Daftar 10 saham LQ45 top berdasarkan kapitalisasi pasar Sumber: IDX Company Fact Sheet LQ45 Periode Feb-Jul 2026 Dipilih saham yang konsisten terdaftar di LQ45 sepanjang 2021-2025

tickers <- c(
  "BBCA.JK",   # Bank Central Asia Tbk
  "BBRI.JK",   # Bank Rakyat Indonesia (Persero) Tbk
  "BMRI.JK",   # Bank Mandiri (Persero) Tbk
  "TLKM.JK",   # Telkom Indonesia (Persero) Tbk
  "ASII.JK",   # Astra International Tbk
  "BBNI.JK",   # Bank Negara Indonesia (Persero) Tbk
  "ANTM.JK",   # Aneka Tambang Tbk
  "UNTR.JK",   # United Tractors Tbk
  "ICBP.JK",   # Indofood CBP Sukses Makmur Tbk
  "ADRO.JK"    # Alamtri Resources Indonesia Tbk (d/h Adaro Energy)
)

start_date <- "2021-01-01"   # awal periode data
end_date   <- "2025-12-31"   # akhir periode data
cl         <- 0.95           # confidence level VaR
alpha      <- 1 - cl         # tingkat signifikansi (0.05)
window     <- 250            # rolling window = 1 tahun perdagangan
n_sim      <- 10000          # jumlah simulasi Monte Carlo
rf_annual  <- 0.06           # risk-free rate tahunan (acuan SBI ~6%)
rf_daily   <- rf_annual / 252

set.seed(42)   # reproducibility untuk Monte Carlo

PENGUNDUHAN DATA DAN PERHITUNGAN LOG RETURN

cat(strrep("=", 60), "\n")
## ============================================================
cat("  BAGIAN 2: MENGUNDUH DATA DARI YAHOO FINANCE\n")
##   BAGIAN 2: MENGUNDUH DATA DARI YAHOO FINANCE
cat(strrep("=", 60), "\n")
## ============================================================
# Download harga saham
getSymbols(
  tickers,
  src         = "yahoo",
  from        = start_date,
  to          = end_date,
  auto.assign = TRUE,
  warnings    = FALSE
)
##  [1] "BBCA.JK" "BBRI.JK" "BMRI.JK" "TLKM.JK" "ASII.JK" "BBNI.JK" "ANTM.JK"
##  [8] "UNTR.JK" "ICBP.JK" "ADRO.JK"
# Ambil Adjusted Close Price dan gabungkan
price_list  <- lapply(tickers, function(t) Ad(get(t)))
prices_all  <- do.call(merge, price_list)

# Bersihkan nama kolom: "BBCA.JK.Adjusted" -> "BBCA"
colnames(prices_all) <- gsub("\\.JK\\.Adjusted", "", colnames(prices_all))
prices_all           <- na.omit(prices_all)

# Hitung log return harian: R_t = ln(P_t / P_{t-1})
returns_all <- na.omit(diff(log(prices_all)))
n_stocks    <- ncol(returns_all)
n_obs       <- nrow(returns_all)
stock_names <- colnames(returns_all)

cat("Saham          :", paste(stock_names, collapse = ", "), "\n")
## Saham          : BBCA, BBRI, BMRI, TLKM, ASII, BBNI, ANTM, UNTR, ICBP, ADRO
cat("Observasi total:", n_obs, "hari\n")
## Observasi total: 1204 hari
cat("Periode        :",
    as.character(index(returns_all)[1]), "s.d.",
    as.character(index(returns_all)[n_obs]), "\n\n")
## Periode        : 2021-01-05 s.d. 2025-12-30
# Download IHSG sebagai indeks pasar untuk SIM
cat("Mengunduh data IHSG (^JKSE) sebagai indeks pasar...\n")
## Mengunduh data IHSG (^JKSE) sebagai indeks pasar...
getSymbols("^JKSE", src = "yahoo",
           from = start_date, to = end_date,
           auto.assign = TRUE, warnings = FALSE)
## [1] "JKSE"
market_ret <- na.omit(diff(log(Ad(JKSE))))

# Sejajarkan tanggal saham dan IHSG
common_dates <- intersect(index(returns_all), index(market_ret))
R_mat        <- as.matrix(returns_all[common_dates, ])
R_mkt        <- as.numeric(market_ret[common_dates])
port_dates   <- as.Date(common_dates)

# Hitung parameter dasar
mu_vec  <- colMeans(R_mat)
Sigma   <- cov(R_mat)
var_mkt <- var(R_mkt)

cat("Observasi setelah alignment:", length(common_dates), "hari\n\n")
## Observasi setelah alignment: 1204 hari

STATISTIK DESKRIPTIF

cat(strrep("=", 60), "\n")
## ============================================================
cat("  BAGIAN 3: STATISTIK DESKRIPTIF LOG RETURN HARIAN\n")
##   BAGIAN 3: STATISTIK DESKRIPTIF LOG RETURN HARIAN
cat(strrep("=", 60), "\n")
## ============================================================
desc_stats <- data.frame(
  Saham    = stock_names,
  Mean     = round(colMeans(R_mat) * 100, 4),
  Std_Dev  = round(apply(R_mat, 2, sd) * 100, 4),
  Min      = round(apply(R_mat, 2, min) * 100, 4),
  Max      = round(apply(R_mat, 2, max) * 100, 4),
  Skewness = round(apply(R_mat, 2, skewness), 4),
  Kurtosis = round(apply(R_mat, 2, kurtosis), 4)
)
print(desc_stats, row.names = FALSE)
##  Saham    Mean Std_Dev      Min     Max Skewness Kurtosis
##   BBCA  0.0243  1.4596  -8.9153  7.3427   0.1829   5.6004
##   BBRI  0.0194  1.8179 -10.6733  8.8251   0.0472   5.5918
##   BMRI  0.0603  1.8622 -10.7500  8.4475  -0.0320   5.9465
##   TLKM  0.0215  1.8117  -6.9457 10.9434   0.3765   5.6961
##   ASII  0.0356  1.7428  -9.3685  9.4858   0.4126   5.3100
##   BBNI  0.0429  1.8975  -8.3382  8.5942   0.2314   4.8927
##   ANTM  0.0462  2.8796 -15.5171 16.3072   0.6789   7.3227
##   UNTR  0.0517  2.1334 -15.8406 14.9375   0.0289   8.5455
##   ICBP -0.0038  1.6444  -9.2932  9.5310   0.0794   6.2514
##   ADRO  0.1147  2.7340 -28.2863 17.7195  -0.1333  16.4659
cat("\nKeterangan: Mean dan Std Dev dalam % per hari\n\n")
## 
## Keterangan: Mean dan Std Dev dalam % per hari

UJI NORMALITAS (JARQUE-BERA)

cat(strrep("=", 60), "\n")
## ============================================================
cat("  BAGIAN 4: UJI NORMALITAS — JARQUE-BERA TEST\n")
##   BAGIAN 4: UJI NORMALITAS — JARQUE-BERA TEST
cat(strrep("=", 60), "\n")
## ============================================================
cat("H0: data berdistribusi normal  |  Tolak H0 jika p-value < 0.05\n\n")
## H0: data berdistribusi normal  |  Tolak H0 jika p-value < 0.05
jb_results <- t(apply(R_mat, 2, function(x) {
  jb <- jarque.test(x)
  c(
    JB_Stat  = round(jb$statistic, 4),
    P_Value  = round(jb$p.value, 6),
    Skewness = round(skewness(x), 4),
    Kurtosis = round(kurtosis(x), 4),
    Normal   = ifelse(jb$p.value > 0.05, "Ya", "Tidak")
  )
}))
print(as.data.frame(jb_results))
##      JB_Stat.JB P_Value Skewness Kurtosis Normal
## BBCA   345.9336       0   0.1829   5.6004  Tidak
## BBRI   337.4467       0   0.0472   5.5918  Tidak
## BMRI   435.7333       0   -0.032   5.9465  Tidak
## TLKM   393.0877       0   0.3765   5.6961  Tidak
## ASII   301.8704       0   0.4126     5.31  Tidak
## BBNI   190.4557       0   0.2314   4.8927  Tidak
## ANTM  1029.9143       0   0.6789   7.3227  Tidak
## UNTR  1542.9389       0   0.0289   8.5455  Tidak
## ICBP   531.6004       0   0.0794   6.2514  Tidak
## ADRO  9100.2825       0  -0.1333  16.4659  Tidak
cat("\nCatatan: Seluruh return berdistribusi tidak normal (leptokurtik),\n")
## 
## Catatan: Seluruh return berdistribusi tidak normal (leptokurtik),
cat("konsisten dengan karakteristik data keuangan (fat tail).\n\n")
## konsisten dengan karakteristik data keuangan (fat tail).

UJI STASIONERITAS (AUGMENTED DICKEY-FULLER)

cat(strrep("=", 60), "\n")
## ============================================================
cat("  BAGIAN 5: UJI STASIONERITAS — AUGMENTED DICKEY-FULLER\n")
##   BAGIAN 5: UJI STASIONERITAS — AUGMENTED DICKEY-FULLER
cat(strrep("=", 60), "\n")
## ============================================================
cat("H0: data tidak stasioner  |  Tolak H0 jika p-value < 0.05\n\n")
## H0: data tidak stasioner  |  Tolak H0 jika p-value < 0.05
adf_results <- t(sapply(stock_names, function(s) {
  pval <- suppressWarnings(adf.test(R_mat[, s])$p.value)
  c(
    ADF_PValue = round(pval, 4),
    Stasioner  = ifelse(pval < 0.05, "Ya", "Tidak")
  )
}))
print(as.data.frame(adf_results))
##      ADF_PValue Stasioner
## BBCA       0.01        Ya
## BBRI       0.01        Ya
## BMRI       0.01        Ya
## TLKM       0.01        Ya
## ASII       0.01        Ya
## BBNI       0.01        Ya
## ANTM       0.01        Ya
## UNTR       0.01        Ya
## ICBP       0.01        Ya
## ADRO       0.01        Ya
cat("\nCatatan: p-value '0.01' menunjukkan p-value <= 0.01 (batas bawah tseries).\n")
## 
## Catatan: p-value '0.01' menunjukkan p-value <= 0.01 (batas bawah tseries).
cat("Seluruh log return bersifat stasioner — memenuhi asumsi analisis.\n\n")
## Seluruh log return bersifat stasioner — memenuhi asumsi analisis.

OPTIMASI PORTOFOLIO: MODEL MARKOWITZ

cat(strrep("=", 60), "\n")
## ============================================================
cat("  BAGIAN 6: OPTIMASI PORTOFOLIO — MODEL MARKOWITZ\n")
##   BAGIAN 6: OPTIMASI PORTOFOLIO — MODEL MARKOWITZ
cat(strrep("=", 60), "\n")
## ============================================================
cat("Metode: Mean-Variance Optimization (Tangency Portfolio / Max Sharpe)\n")
## Metode: Mean-Variance Optimization (Tangency Portfolio / Max Sharpe)
cat("Constraints: bobot non-negatif, jumlah bobot = 1\n\n")
## Constraints: bobot non-negatif, jumlah bobot = 1
markowitz_optimize <- function(mu, Sigma, rf = rf_daily) {
  n    <- length(mu)
  Dmat <- 2 * Sigma
  dvec <- rep(0, n)

  # Constraints: sum(w) = 1 (equality) dan w >= 0 (inequality)
  Amat <- cbind(rep(1, n), diag(n))
  bvec <- c(1, rep(0, n))

  # Buat 100 target return untuk menyusun efficient frontier
  targets      <- seq(min(mu), max(mu), length.out = 100)
  frontier     <- data.frame(
    port_return = numeric(0),
    port_sd     = numeric(0),
    sharpe      = numeric(0)
  )
  weights_list <- list()

  for (tr in targets) {
    tryCatch({
      sol <- solve.QP(Dmat, dvec,
                      cbind(Amat, mu), c(bvec, tr), meq = 1)
      w   <- pmax(sol$solution, 0)
      w   <- w / sum(w)
      ret <- sum(w * mu)
      sig <- as.numeric(sqrt(t(w) %*% Sigma %*% w))
      sr  <- (ret - rf) / sig

      frontier <- rbind(frontier, data.frame(
        port_return = ret, port_sd = sig, sharpe = sr
      ))
      weights_list[[length(weights_list) + 1]] <- w
    }, error = function(e) NULL)
  }

  best_idx <- which.max(frontier$sharpe)
  list(
    frontier = frontier,
    weights  = weights_list[[best_idx]],
    ret      = frontier$port_return[best_idx],
    sd       = frontier$port_sd[best_idx],
    sharpe   = frontier$sharpe[best_idx]
  )
}

mw_result <- markowitz_optimize(mu_vec, Sigma)
w_mw      <- mw_result$weights
ret_mw    <- mw_result$ret
sd_mw     <- mw_result$sd
sr_mw     <- mw_result$sharpe

# Tampilkan bobot (semua saham, urut dari terbesar)
cat("--- Bobot Portofolio Markowitz (Tangency Portfolio) ---\n")
## --- Bobot Portofolio Markowitz (Tangency Portfolio) ---
mw_w_df <- data.frame(
  Saham     = stock_names,
  Bobot_pct = round(w_mw * 100, 4)
)
print(mw_w_df[order(-mw_w_df$Bobot_pct), ], row.names = FALSE)
##  Saham Bobot_pct
##   ADRO   62.5988
##   BMRI   37.4012
##   BBCA    0.0000
##   BBRI    0.0000
##   TLKM    0.0000
##   ASII    0.0000
##   BBNI    0.0000
##   ANTM    0.0000
##   UNTR    0.0000
##   ICBP    0.0000
cat(sprintf("\nExpected Return Harian  : %.6f (%.4f%%)\n", ret_mw, ret_mw * 100))
## 
## Expected Return Harian  : 0.000944 (0.0944%)
cat(sprintf("Expected Return Tahunan : %.4f%%\n", ret_mw * 252 * 100))
## Expected Return Tahunan : 23.7795%
cat(sprintf("Std Deviation Harian    : %.6f (%.4f%%)\n", sd_mw, sd_mw * 100))
## Std Deviation Harian    : 0.019961 (1.9961%)
cat(sprintf("Std Deviation Tahunan   : %.4f%%\n", sd_mw * sqrt(252) * 100))
## Std Deviation Tahunan   : 31.6868%
cat(sprintf("Sharpe Ratio (Harian)   : %.4f\n", sr_mw))
## Sharpe Ratio (Harian)   : 0.0353
cat(sprintf("Sharpe Ratio (Tahunan)  : %.4f\n\n", sr_mw * sqrt(252)))
## Sharpe Ratio (Tahunan)  : 0.5611

OPTIMASI PORTOFOLIO: SINGLE INDEX MODEL (SIM)

cat(strrep("=", 60), "\n")
## ============================================================
cat("  BAGIAN 7: OPTIMASI PORTOFOLIO — SINGLE INDEX MODEL (SIM)\n")
##   BAGIAN 7: OPTIMASI PORTOFOLIO — SINGLE INDEX MODEL (SIM)
cat(strrep("=", 60), "\n")
## ============================================================
cat("Metode: Elton-Gruber-Padberg (Excess Return to Beta / Cut-off Rate)\n")
## Metode: Elton-Gruber-Padberg (Excess Return to Beta / Cut-off Rate)
cat("Indeks Pasar: IHSG (^JKSE)\n\n")
## Indeks Pasar: IHSG (^JKSE)
# -- Langkah 7.1: Estimasi parameter SIM per saham --
cat("--- Langkah 7.1: Parameter SIM per Saham ---\n")
## --- Langkah 7.1: Parameter SIM per Saham ---
sim_par <- data.frame(
  Saham     = stock_names,
  Alpha_i   = NA_real_,
  Beta_i    = NA_real_,
  Sigma_ei2 = NA_real_,  # varians residual (risiko unik)
  Mu_i      = NA_real_,  # expected return harian
  ERB       = NA_real_   # Excess Return to Beta
)

for (i in seq_len(n_stocks)) {
  fit                  <- lm(R_mat[, i] ~ R_mkt)
  sim_par$Alpha_i[i]   <- coef(fit)[1]
  sim_par$Beta_i[i]    <- coef(fit)[2]
  sim_par$Sigma_ei2[i] <- var(residuals(fit))
  sim_par$Mu_i[i]      <- mean(R_mat[, i])
  sim_par$ERB[i]       <- (sim_par$Mu_i[i] - rf_daily) / sim_par$Beta_i[i]
}

print(
  cbind(
    Saham = sim_par$Saham,
    round(sim_par[, c("Alpha_i", "Beta_i", "Sigma_ei2", "ERB")], 6)
  ),
  row.names = FALSE
)
##  Saham   Alpha_i   Beta_i Sigma_ei2       ERB
##   BBCA -0.000045 0.994418  0.000135  0.000005
##   BBRI -0.000190 1.329659  0.000190 -0.000033
##   BMRI  0.000213 1.350011  0.000202  0.000270
##   TLKM -0.000044 0.893509  0.000265 -0.000026
##   ASII  0.000092 0.915015  0.000237  0.000129
##   BBNI  0.000051 1.307372  0.000224  0.000146
##   ANTM  0.000128 1.154791  0.000723  0.000194
##   UNTR  0.000266 0.865526  0.000396  0.000322
##   ICBP -0.000171 0.461007  0.000254 -0.000598
##   ADRO  0.000820 1.129963  0.000646  0.000804
# -- Langkah 7.2: Hitung Cut-off Rate (C*) --
cat("\n--- Langkah 7.2: Penentuan Cut-off Rate (C*) ---\n")
## 
## --- Langkah 7.2: Penentuan Cut-off Rate (C*) ---
sim_s  <- sim_par[order(sim_par$ERB, decreasing = TRUE), ]
C_num  <- cumsum((sim_s$ERB - rf_daily) * sim_s$Beta_i / sim_s$Sigma_ei2)
C_den  <- 1 + var_mkt * cumsum(sim_s$Beta_i^2 / sim_s$Sigma_ei2)
C_vec  <- var_mkt * C_num / C_den

cutoff_idx <- max(which(sim_s$ERB > C_vec))
C_star     <- C_vec[cutoff_idx]
selected   <- sim_s[seq_len(cutoff_idx), ]

cat(sprintf("Cut-off Rate (C*)   : %.6f\n", C_star))
## Cut-off Rate (C*)   : -0.000071
cat(sprintf("Jumlah saham terpilih: %d saham\n", nrow(selected)))
## Jumlah saham terpilih: 9 saham
cat(sprintf("Saham terpilih SIM  : %s\n",
            paste(selected$Saham, collapse = ", ")))
## Saham terpilih SIM  : ADRO, UNTR, BMRI, ANTM, BBNI, ASII, BBCA, TLKM, BBRI
# -- Langkah 7.3: Hitung bobot SIM --
cat("\n--- Langkah 7.3: Bobot Portofolio SIM ---\n")
## 
## --- Langkah 7.3: Bobot Portofolio SIM ---
Z_i   <- (selected$Beta_i / selected$Sigma_ei2) * (selected$ERB - C_star)
w_sel <- Z_i / sum(Z_i)

w_sim <- rep(0, n_stocks)
for (i in seq_len(nrow(selected))) {
  idx       <- which(stock_names == selected$Saham[i])
  w_sim[idx] <- w_sel[i]
}

ret_sim <- sum(w_sim * mu_vec)
sd_sim  <- as.numeric(sqrt(t(w_sim) %*% Sigma %*% w_sim))
sr_sim  <- (ret_sim - rf_daily) / sd_sim

sim_w_df <- data.frame(
  Saham     = stock_names,
  Bobot_pct = round(w_sim * 100, 4)
)
print(sim_w_df[order(-sim_w_df$Bobot_pct), ], row.names = FALSE)
##  Saham Bobot_pct
##   BMRI   28.1162
##   ADRO   18.8613
##   BBNI   15.5852
##   UNTR   10.5897
##   ASII    9.5355
##   BBCA    6.9264
##   ANTM    5.2098
##   BBRI    3.2935
##   TLKM    1.8825
##   ICBP    0.0000
cat(sprintf("\nExpected Return Harian  : %.6f (%.4f%%)\n", ret_sim, ret_sim * 100))
## 
## Expected Return Harian  : 0.000593 (0.0593%)
cat(sprintf("Expected Return Tahunan : %.4f%%\n", ret_sim * 252 * 100))
## Expected Return Tahunan : 14.9372%
cat(sprintf("Std Deviation Harian    : %.6f (%.4f%%)\n", sd_sim, sd_sim * 100))
## Std Deviation Harian    : 0.013397 (1.3397%)
cat(sprintf("Std Deviation Tahunan   : %.4f%%\n", sd_sim * sqrt(252) * 100))
## Std Deviation Tahunan   : 21.2666%
cat(sprintf("Sharpe Ratio (Harian)   : %.4f\n", sr_sim))
## Sharpe Ratio (Harian)   : 0.0265
cat(sprintf("Sharpe Ratio (Tahunan)  : %.4f\n\n", sr_sim * sqrt(252)))
## Sharpe Ratio (Tahunan)  : 0.4202

PERHITUNGAN RETURN PORTOFOLIO HARIAN

port_ret_mw  <- as.numeric(R_mat %*% w_mw)
port_ret_sim <- as.numeric(R_mat %*% w_sim)

FUNGSI VALUE AT RISK (VaR) 95%

# -- VaR Historical Simulation --
# Metode non-parametrik: ambil persentil ke-alpha dari return historis
var_hs_fn <- function(ret) {
  as.numeric(-quantile(ret, alpha))
}

# -- VaR Monte Carlo --
# Bangkitkan n_sim skenario return dari distribusi Normal(mean, sd)
# set.seed(42) sudah dipanggil di Bagian 1 untuk reproducibility
var_mc_fn <- function(ret, n = n_sim) {
  sim <- rnorm(n, mean = mean(ret), sd = sd(ret))
  as.numeric(-quantile(sim, alpha))
}

ROLLING VaR (WINDOW = 250 HARI)

cat(strrep("=", 60), "\n")
## ============================================================
cat(sprintf("  BAGIAN 10: ROLLING VaR 95%% (WINDOW = %d HARI)\n", window))
##   BAGIAN 10: ROLLING VaR 95% (WINDOW = 250 HARI)
cat(strrep("=", 60), "\n\n")
## ============================================================
compute_rolling_var <- function(port_ret) {
  n_total <- length(port_ret)
  n_roll  <- n_total - window
  hs      <- numeric(n_roll)
  mc      <- numeric(n_roll)

  for (i in seq_len(n_roll)) {
    w_data <- port_ret[i:(i + window - 1)]
    hs[i]  <- var_hs_fn(w_data)
    mc[i]  <- var_mc_fn(w_data)
  }

  data.frame(
    date   = port_dates[(window + 1):n_total],
    actual = port_ret[(window + 1):n_total],
    HS     = hs,
    MC     = mc
  )
}

cat("Menghitung rolling VaR portofolio Markowitz...\n")
## Menghitung rolling VaR portofolio Markowitz...
var_mw  <- compute_rolling_var(port_ret_mw)

cat("Menghitung rolling VaR portofolio SIM...\n")
## Menghitung rolling VaR portofolio SIM...
var_sim <- compute_rolling_var(port_ret_sim)

cat(sprintf("Periode out-of-sample : %d hari\n\n", nrow(var_mw)))
## Periode out-of-sample : 954 hari

BACKTESTING VaR

cat(strrep("=", 60), "\n")
## ============================================================
cat("  BAGIAN 11: BACKTESTING VaR\n")
##   BAGIAN 11: BACKTESTING VaR
cat(strrep("=", 60), "\n\n")
## ============================================================
# Fungsi penentu exception (pelanggaran VaR)
# Exception = 1 jika return aktual < -VaR (kerugian melebihi estimasi)
get_hits <- function(actual, var_est) {
  as.integer(actual < -var_est)
}

# ── Uji Kupiec (POF Test) ──────────────────────────────────
# H0: proporsi exception = alpha (5%)
# Statistik: LR_POF ~ chi-sq(df=1), nilai kritis = 3.841
kupiec_test <- function(hits) {
  T_n   <- length(hits)
  x     <- sum(hits)
  p     <- alpha

  if (x == 0 || x == T_n) {
    return(data.frame(
      Exceptions = x,
      Exp_Exc    = round(T_n * p, 1),
      Pct_Exc    = round(x / T_n * 100, 2),
      LR_Stat    = NA,
      Critical   = 3.841,
      P_Value    = NA,
      Keputusan  = "Data tidak cukup"
    ))
  }

  p_hat <- x / T_n
  LR    <- -2 * (
    x       * log(p / p_hat) +
    (T_n - x) * log((1 - p) / (1 - p_hat))
  )
  pval  <- 1 - pchisq(LR, df = 1)

  data.frame(
    Exceptions = x,
    Exp_Exc    = round(T_n * p, 1),
    Pct_Exc    = round(p_hat * 100, 2),
    LR_Stat    = round(LR, 4),
    Critical   = 3.841,
    P_Value    = round(pval, 4),
    Keputusan  = ifelse(pval > 0.05, "Diterima (Valid)", "Ditolak (Tidak Valid)")
  )
}

# ── Uji Christoffersen (Conditional Coverage Test) ────────
# H0: exception bersifat independen DAN proporsinya = alpha
# Statistik: LR_CC = LR_UC + LR_Ind ~ chi-sq(df=2), nilai kritis = 5.991
christoffersen_test <- function(hits) {
  T_n   <- length(hits)
  x     <- sum(hits)
  p     <- alpha

  if (x == 0 || x == T_n) {
    return(data.frame(
      Exceptions = x,
      Exp_Exc    = round(T_n * p, 1),
      Pct_Exc    = round(x / T_n * 100, 2),
      LR_UC      = NA,
      LR_Ind     = NA,
      LR_CC      = NA,
      Critical   = 5.991,
      P_Value    = NA,
      Keputusan  = "Data tidak cukup"
    ))
  }

  # Unconditional Coverage (identik dengan Kupiec)
  p_hat <- x / T_n
  LR_uc <- -2 * (
    x       * log(p / p_hat) +
    (T_n - x) * log((1 - p) / (1 - p_hat))
  )

  # Matriks transisi: n_ij = transisi dari kondisi i ke j hari berikutnya
  h_lag  <- hits[-T_n]
  h_curr <- hits[-1]
  n00    <- sum(h_lag == 0 & h_curr == 0)
  n01    <- sum(h_lag == 0 & h_curr == 1)
  n10    <- sum(h_lag == 1 & h_curr == 0)
  n11    <- sum(h_lag == 1 & h_curr == 1)

  # Probabilitas transisi (dengan clamp untuk hindari log(0))
  p01 <- ifelse((n00 + n01) > 0, n01 / (n00 + n01), 1e-10)
  p11 <- ifelse((n10 + n11) > 0, n11 / (n10 + n11), 1e-10)
  pi2 <- (n01 + n11) / (n00 + n01 + n10 + n11)

  p01 <- pmax(pmin(p01, 1 - 1e-10), 1e-10)
  p11 <- pmax(pmin(p11, 1 - 1e-10), 1e-10)
  pi2 <- pmax(pmin(pi2, 1 - 1e-10), 1e-10)

  # Statistik Independence
  LR_ind <- -2 * (
    (n00 + n10) * log(1 - pi2) + (n01 + n11) * log(pi2) -
     n00 * log(1 - p01) - n01 * log(p01) -
     n10 * log(1 - p11) - n11 * log(p11)
  )

  LR_cc <- LR_uc + LR_ind
  pval  <- 1 - pchisq(LR_cc, df = 2)

  data.frame(
    Exceptions = x,
    Exp_Exc    = round(T_n * p, 1),
    Pct_Exc    = round(p_hat * 100, 2),
    LR_UC      = round(LR_uc, 4),
    LR_Ind     = round(LR_ind, 4),
    LR_CC      = round(LR_cc, 4),
    Critical   = 5.991,
    P_Value    = round(pval, 4),
    Keputusan  = ifelse(pval > 0.05, "Diterima (Valid)", "Ditolak (Tidak Valid)")
  )
}

# ── Fungsi wrapper: jalankan backtesting untuk satu portofolio ──
run_backtesting <- function(var_df, portfolio_name) {
  methods <- list(
    list(col = "HS", label = "Historical Simulation 95%"),
    list(col = "MC", label = "Monte Carlo 95%")
  )

  cat(strrep("-", 60), "\n")
  cat(sprintf("  Portofolio: %s\n", portfolio_name))
  cat(strrep("-", 60), "\n")

  kupiec_rows <- list()
  christ_rows <- list()

  for (m in methods) {
    hits                   <- get_hits(var_df$actual, var_df[[m$col]])
    kupiec_rows[[m$label]] <- cbind(Metode = m$label, kupiec_test(hits))
    christ_rows[[m$label]] <- cbind(Metode = m$label, christoffersen_test(hits))
  }

  k_tbl <- do.call(rbind, kupiec_rows)
  c_tbl <- do.call(rbind, christ_rows)
  rownames(k_tbl) <- rownames(c_tbl) <- NULL

  cat("\n[Kupiec POF Test]  chi-sq df=1 | Nilai kritis = 3.841\n")
  print(k_tbl, row.names = FALSE)

  cat("\n[Christoffersen Conditional Coverage]  chi-sq df=2 | Nilai kritis = 5.991\n")
  print(c_tbl, row.names = FALSE)

  cat("\nKeterangan: Diterima = model VaR valid (p-value > 0.05)\n\n")
  invisible(list(kupiec = k_tbl, christoffersen = c_tbl))
}

# Jalankan backtesting untuk masing-masing portofolio
bt_mw  <- run_backtesting(var_mw,  "MARKOWITZ")
## ------------------------------------------------------------ 
##   Portofolio: MARKOWITZ
## ------------------------------------------------------------ 
## 
## [Kupiec POF Test]  chi-sq df=1 | Nilai kritis = 3.841
##                     Metode Exceptions Exp_Exc Pct_Exc LR_Stat Critical P_Value
##  Historical Simulation 95%         44    47.7    4.61  0.3098    3.841  0.5778
##            Monte Carlo 95%         38    47.7    3.98  2.2252    3.841  0.1358
##         Keputusan
##  Diterima (Valid)
##  Diterima (Valid)
## 
## [Christoffersen Conditional Coverage]  chi-sq df=2 | Nilai kritis = 5.991
##                     Metode Exceptions Exp_Exc Pct_Exc  LR_UC LR_Ind  LR_CC
##  Historical Simulation 95%         44    47.7    4.61 0.3098 3.5101 3.8199
##            Monte Carlo 95%         38    47.7    3.98 2.2252 0.1540 2.3792
##  Critical P_Value        Keputusan
##     5.991  0.1481 Diterima (Valid)
##     5.991  0.3043 Diterima (Valid)
## 
## Keterangan: Diterima = model VaR valid (p-value > 0.05)
bt_sim <- run_backtesting(var_sim, "SINGLE INDEX MODEL (SIM)")
## ------------------------------------------------------------ 
##   Portofolio: SINGLE INDEX MODEL (SIM)
## ------------------------------------------------------------ 
## 
## [Kupiec POF Test]  chi-sq df=1 | Nilai kritis = 3.841
##                     Metode Exceptions Exp_Exc Pct_Exc LR_Stat Critical P_Value
##  Historical Simulation 95%         54    47.7    5.66  0.8416    3.841  0.3589
##            Monte Carlo 95%         48    47.7    5.03  0.0020    3.841  0.9645
##         Keputusan
##  Diterima (Valid)
##  Diterima (Valid)
## 
## [Christoffersen Conditional Coverage]  chi-sq df=2 | Nilai kritis = 5.991
##                     Metode Exceptions Exp_Exc Pct_Exc  LR_UC LR_Ind  LR_CC
##  Historical Simulation 95%         54    47.7    5.66 0.8416 0.2991 1.1407
##            Monte Carlo 95%         48    47.7    5.03 0.0020 0.1455 0.1475
##  Critical P_Value        Keputusan
##     5.991  0.5653 Diterima (Valid)
##     5.991  0.9289 Diterima (Valid)
## 
## Keterangan: Diterima = model VaR valid (p-value > 0.05)
# ── Tabel ringkasan backtesting semua kombinasi ──
cat(strrep("=", 60), "\n")
## ============================================================
cat("  RINGKASAN BACKTESTING — SEMUA KOMBINASI\n")
##   RINGKASAN BACKTESTING — SEMUA KOMBINASI
cat(strrep("=", 60), "\n")
## ============================================================
ringkasan_bt <- data.frame(
  Portofolio = c("Markowitz", "Markowitz", "SIM", "SIM"),
  Metode_VaR = c("Historical Simulation", "Monte Carlo",
                 "Historical Simulation", "Monte Carlo"),
  Kupiec     = c(
    bt_mw$kupiec[1, "Keputusan"],
    bt_mw$kupiec[2, "Keputusan"],
    bt_sim$kupiec[1, "Keputusan"],
    bt_sim$kupiec[2, "Keputusan"]
  ),
  Christoffersen = c(
    bt_mw$christoffersen[1, "Keputusan"],
    bt_mw$christoffersen[2, "Keputusan"],
    bt_sim$christoffersen[1, "Keputusan"],
    bt_sim$christoffersen[2, "Keputusan"]
  )
)
ringkasan_bt$Kesimpulan <- ifelse(
  ringkasan_bt$Kupiec        == "Diterima (Valid)" &
  ringkasan_bt$Christoffersen == "Diterima (Valid)",
  "Valid", "Tidak Valid"
)
print(ringkasan_bt, row.names = FALSE)
##  Portofolio            Metode_VaR           Kupiec   Christoffersen Kesimpulan
##   Markowitz Historical Simulation Diterima (Valid) Diterima (Valid)      Valid
##   Markowitz           Monte Carlo Diterima (Valid) Diterima (Valid)      Valid
##         SIM Historical Simulation Diterima (Valid) Diterima (Valid)      Valid
##         SIM           Monte Carlo Diterima (Valid) Diterima (Valid)      Valid
cat("\n")

RINGKASAN PERBANDINGAN PORTOFOLIO

cat(strrep("=", 60), "\n")
## ============================================================
cat("  BAGIAN 12: RINGKASAN PERBANDINGAN PORTOFOLIO\n")
##   BAGIAN 12: RINGKASAN PERBANDINGAN PORTOFOLIO
cat(strrep("=", 60), "\n\n")
## ============================================================
summary_df <- data.frame(
  Model              = c("Markowitz", "SIM"),
  Ret_Harian_pct     = round(c(ret_mw, ret_sim) * 100, 4),
  Ret_Tahunan_pct    = round(c(ret_mw, ret_sim) * 252 * 100, 4),
  SD_Harian_pct      = round(c(sd_mw, sd_sim) * 100, 4),
  SD_Tahunan_pct     = round(c(sd_mw, sd_sim) * sqrt(252) * 100, 4),
  Sharpe_Harian      = round(c(sr_mw, sr_sim), 4),
  Sharpe_Tahunan     = round(c(sr_mw, sr_sim) * sqrt(252), 4),
  VaR_HS_avg_pct     = round(c(mean(var_mw$HS), mean(var_sim$HS)) * 100, 4),
  VaR_MC_avg_pct     = round(c(mean(var_mw$MC), mean(var_sim$MC)) * 100, 4)
)
print(summary_df, row.names = FALSE)
##      Model Ret_Harian_pct Ret_Tahunan_pct SD_Harian_pct SD_Tahunan_pct
##  Markowitz         0.0944         23.7795        1.9961        31.6868
##        SIM         0.0593         14.9372        1.3397        21.2666
##  Sharpe_Harian Sharpe_Tahunan VaR_HS_avg_pct VaR_MC_avg_pct
##         0.0353         0.5611         2.8776         3.1283
##         0.0265         0.4202         1.9331         2.0454
cat("\nKeterangan: Ret = Expected Return  |  SD = Std Deviation\n")
## 
## Keterangan: Ret = Expected Return  |  SD = Std Deviation
cat("VaR dihitung pada confidence level 95%\n\n")
## VaR dihitung pada confidence level 95%

VISUALISASI

cat(strrep("=", 60), "\n")
## ============================================================
cat("  BAGIAN 13: MEMBUAT VISUALISASI\n")
##   BAGIAN 13: MEMBUAT VISUALISASI
cat(strrep("=", 60), "\n\n")
## ============================================================
# -- Tema jurnal: bersih, hitam-putih friendly --
theme_jurnal <- theme_classic(base_size = 11) +
  theme(
    plot.title       = element_text(face = "bold", size = 12),
    plot.subtitle    = element_text(color = "gray40", size = 9),
    plot.caption     = element_text(color = "gray50", size = 8, hjust = 0),
    legend.position  = "top",
    legend.title     = element_text(face = "bold", size = 9),
    legend.key.size  = unit(0.4, "cm"),
    axis.title       = element_text(size = 10),
    panel.grid.major = element_line(color = "gray90", linewidth = 0.4),
    panel.grid.minor = element_blank()
  )

warna_model <- c("Markowitz" = "#1565C0", "SIM" = "#C62828")
warna_var   <- c("Historical Simulation" = "#1565C0",
                 "Monte Carlo"           = "#2E7D32")

# ── Plot 1: Perkembangan Harga Saham (Normalized, Base = 100) ──
cat("Membuat Plot 1: Harga saham normalized...\n")
## Membuat Plot 1: Harga saham normalized...
prices_norm <- sweep(as.matrix(prices_all), 2,
                     as.numeric(prices_all[1, ]), "/") * 100
prices_df        <- as.data.frame(prices_norm)
prices_df$date   <- as.Date(index(prices_all))
prices_long      <- pivot_longer(prices_df, -date,
                                 names_to = "Saham", values_to = "Harga")

p1 <- ggplot(prices_long, aes(x = date, y = Harga, color = Saham)) +
  geom_line(linewidth = 0.5, alpha = 0.85) +
  geom_hline(yintercept = 100, linetype = "dashed",
             color = "gray40", linewidth = 0.4) +
  scale_x_date(date_breaks = "6 months", date_labels = "%b %Y") +
  scale_color_brewer(palette = "Paired") +
  labs(
    title    = "Perkembangan Harga Saham LQ45 (Normalized, Base = 100)",
    subtitle = paste0("Periode: ", start_date, " s.d. ", end_date),
    x        = NULL,
    y        = "Indeks Harga (Base = 100)",
    color    = "Saham",
    caption  = "Sumber: Yahoo Finance"
  ) +
  theme_jurnal +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 8))

ggsave("01_harga_saham_normalized.png",
       plot = p1, width = 11, height = 6, dpi = 150, bg = "white")

# ── Plot 2: Heatmap Korelasi Return ──
cat("Membuat Plot 2: Heatmap korelasi...\n")
## Membuat Plot 2: Heatmap korelasi...
p2 <- ggcorrplot(
  cor(R_mat),
  method        = "square",
  type          = "lower",
  lab           = TRUE,
  lab_size      = 3,
  colors        = c("#C62828", "white", "#1565C0"),
  tl.cex        = 9,
  tl.col        = "black",
  outline.color = "white"
) +
  labs(
    title    = "Heatmap Korelasi Log Return Harian — 10 Saham LQ45",
    subtitle = paste0("Periode: ", start_date, " s.d. ", end_date),
    caption  = "Biru = korelasi positif kuat  |  Merah = korelasi negatif"
  ) +
  theme(
    plot.title    = element_text(face = "bold", size = 12),
    plot.subtitle = element_text(color = "gray40", size = 9),
    plot.caption  = element_text(color = "gray50", size = 8, hjust = 0)
  )

ggsave("02_heatmap_korelasi.png",
       plot = p2, width = 9, height = 8, dpi = 150, bg = "white")

# ── Plot 3: Efficient Frontier Markowitz ──
cat("Membuat Plot 3: Efficient frontier...\n")
## Membuat Plot 3: Efficient frontier...
saham_pts <- data.frame(
  sd_pct  = apply(R_mat, 2, sd) * 100,
  ret_pct = colMeans(R_mat) * 100,
  Saham   = stock_names
)

p3 <- ggplot(mw_result$frontier,
             aes(x = port_sd * 100, y = port_return * 100)) +
  geom_path(color = "#1565C0", linewidth = 1.3) +
  geom_point(data = saham_pts,
             aes(x = sd_pct, y = ret_pct),
             color = "gray50", size = 2, shape = 1) +
  geom_text(data = saham_pts,
            aes(x = sd_pct, y = ret_pct, label = Saham),
            size = 2.5, hjust = -0.15, color = "gray40") +
  geom_point(data = data.frame(x = sd_mw * 100, y = ret_mw * 100),
             aes(x = x, y = y),
             color = "#C62828", size = 5, shape = 18) +
  annotate("text",
           x = sd_mw * 100 + 0.03, y = ret_mw * 100,
           label = "Tangency\nPortfolio",
           color = "#C62828", size = 3, hjust = 0, fontface = "bold") +
  labs(
    title    = "Efficient Frontier — Model Markowitz",
    subtitle = "◆ = Tangency Portfolio (Max Sharpe Ratio)  |  ○ = Saham Individual",
    x        = "Risiko: Std Deviation Harian (%)",
    y        = "Expected Return Harian (%)",
    caption  = paste0("Periode estimasi: ", start_date, " s.d. ", end_date)
  ) +
  theme_jurnal

ggsave("03_efficient_frontier.png",
       plot = p3, width = 9, height = 6, dpi = 150, bg = "white")

# ── Plot 4: Perbandingan Bobot Portofolio ──
cat("Membuat Plot 4: Perbandingan bobot...\n")
## Membuat Plot 4: Perbandingan bobot...
weights_df <- data.frame(
  Saham = rep(stock_names, 2),
  Bobot = c(w_mw * 100, w_sim * 100),
  Model = rep(c("Markowitz", "SIM"), each = n_stocks)
)

p4 <- ggplot(weights_df,
             aes(x = reorder(Saham, Bobot), y = Bobot, fill = Model)) +
  geom_col(position = "dodge", width = 0.7,
           color = "white", linewidth = 0.3) +
  geom_text(
    aes(label = ifelse(Bobot > 0.5, paste0(round(Bobot, 1), "%"), "")),
    position = position_dodge(0.7),
    hjust = -0.1, size = 2.8
  ) +
  coord_flip() +
  scale_fill_manual(values = warna_model) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(
    title    = "Perbandingan Bobot Portofolio Optimal",
    subtitle = "Markowitz (Tangency Portfolio) vs Single Index Model (SIM)",
    x        = NULL,
    y        = "Bobot (%)",
    fill     = "Model Portofolio"
  ) +
  theme_jurnal

ggsave("04_perbandingan_bobot.png",
       plot = p4, width = 9, height = 5.5, dpi = 150, bg = "white")

# ── Plot 5: Beta Saham — SIM ──
cat("Membuat Plot 5: Beta saham SIM...\n")
## Membuat Plot 5: Beta saham SIM...
beta_df <- sim_par[, c("Saham", "Beta_i", "ERB")]
beta_df$Terpilih <- ifelse(
  beta_df$Saham %in% selected$Saham, "Terpilih", "Tidak Terpilih"
)

p5 <- ggplot(beta_df,
             aes(x = reorder(Saham, Beta_i), y = Beta_i, fill = Terpilih)) +
  geom_col(width = 0.65, color = "white") +
  geom_hline(yintercept = 1, linetype = "dashed",
             color = "#C62828", linewidth = 0.6) +
  geom_text(aes(label = round(Beta_i, 3)), hjust = -0.15, size = 3) +
  coord_flip() +
  scale_fill_manual(values = c("Terpilih"       = "#1565C0",
                               "Tidak Terpilih" = "gray70")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  annotate("text", x = 0.6, y = 1.05,
           label = "β = 1", color = "#C62828", size = 3, hjust = 0) +
  labs(
    title    = "Nilai Beta (β) Saham — Single Index Model",
    subtitle = "Biru = saham terpilih  |  Garis merah putus-putus = β pasar (β = 1)",
    x        = NULL,
    y        = "Beta (β)",
    fill     = NULL,
    caption  = "β > 1: lebih agresif dari pasar  |  β < 1: lebih defensif dari pasar"
  ) +
  theme_jurnal

ggsave("05_beta_saham_sim.png",
       plot = p5, width = 9, height = 5.5, dpi = 150, bg = "white")

# ── Plot 6: Return Harian Portofolio Markowitz vs SIM ──
cat("Membuat Plot 6: Return harian portofolio...\n")
## Membuat Plot 6: Return harian portofolio...
ret_compare_df <- data.frame(
  date      = port_dates,
  Markowitz = port_ret_mw * 100,
  SIM       = port_ret_sim * 100
) |> pivot_longer(-date, names_to = "Model", values_to = "Return")

p6 <- ggplot(ret_compare_df, aes(x = date, y = Return, color = Model)) +
  geom_line(linewidth = 0.35, alpha = 0.85) +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.4) +
  scale_color_manual(values = warna_model) +
  scale_x_date(date_breaks = "6 months", date_labels = "%b %Y") +
  labs(
    title    = "Return Harian Portofolio: Markowitz vs SIM",
    subtitle = paste0("Log Return Harian (%)  |  Periode: ",
                      start_date, " s.d. ", end_date),
    x        = NULL,
    y        = "Log Return Harian (%)",
    color    = "Model Portofolio"
  ) +
  theme_jurnal +
  theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 8))

ggsave("06_return_harian_portofolio.png",
       plot = p6, width = 11, height = 5, dpi = 150, bg = "white")

# ── Plot 7: Distribusi Return Portofolio ──
cat("Membuat Plot 7: Distribusi return portofolio...\n")
## Membuat Plot 7: Distribusi return portofolio...
ret_dist_df <- data.frame(
  Return = c(port_ret_mw * 100, port_ret_sim * 100),
  Model  = rep(c("Markowitz", "SIM"), each = length(port_ret_mw))
)

p7 <- ggplot(ret_dist_df, aes(x = Return, fill = Model, color = Model)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 60, alpha = 0.4,
                 position = "identity", linewidth = 0.2) +
  # Kurva normal Markowitz
  geom_function(
    fun   = dnorm,
    args  = list(mean = mean(port_ret_mw * 100),
                 sd   = sd(port_ret_mw * 100)),
    color = "#1565C0", linewidth = 0.8, linetype = "solid",
    inherit.aes = FALSE
  ) +
  # Kurva normal SIM
  geom_function(
    fun   = dnorm,
    args  = list(mean = mean(port_ret_sim * 100),
                 sd   = sd(port_ret_sim * 100)),
    color = "#C62828", linewidth = 0.8, linetype = "dashed",
    inherit.aes = FALSE
  ) +
  # Garis VaR HS rata-rata
  geom_vline(xintercept = -mean(var_mw$HS)  * 100,
             color = "#1565C0", linewidth = 0.8, linetype = "dotted") +
  geom_vline(xintercept = -mean(var_sim$HS) * 100,
             color = "#C62828", linewidth = 0.8, linetype = "dotted") +
  scale_fill_manual(values  = warna_model) +
  scale_color_manual(values = warna_model) +
  labs(
    title    = "Distribusi Return Harian Portofolio",
    subtitle = "Histogram + kurva normal  |  Garis titik-titik = rata-rata VaR HS 95%",
    x        = "Log Return Harian (%)",
    y        = "Densitas",
    fill     = "Model",
    color    = "Model"
  ) +
  theme_jurnal

ggsave("07_distribusi_return.png",
       plot = p7, width = 9, height = 5.5, dpi = 150, bg = "white")

# ── Plot 8: Rolling VaR 95% — Panel Markowitz & SIM ──
cat("Membuat Plot 8: Rolling VaR backtesting...\n")
## Membuat Plot 8: Rolling VaR backtesting...
make_var_plot <- function(var_df, port_name) {
  df_long <- var_df |>
    pivot_longer(cols     = c(HS, MC),
                 names_to = "Metode", values_to = "VaR") |>
    mutate(Metode = recode(Metode,
                           "HS" = "Historical Simulation",
                           "MC" = "Monte Carlo"))

  exc_pts <- var_df[var_df$actual < -var_df$HS, ]

  ggplot(df_long, aes(x = date)) +
    geom_line(aes(y = actual * 100),
              color = "gray60", linewidth = 0.3, alpha = 0.8) +
    geom_line(aes(y = -VaR * 100, color = Metode), linewidth = 0.7) +
    geom_point(data = exc_pts,
               aes(x = date, y = actual * 100),
               color = "#C62828", size = 1.2, alpha = 0.8,
               shape = 16, inherit.aes = FALSE) +
    scale_color_manual(values = warna_var) +
    scale_x_date(date_breaks = "6 months", date_labels = "%b %Y") +
    labs(
      title    = paste0("Rolling VaR 95% — Portofolio ", port_name),
      subtitle = paste0("Garis berwarna = batas VaR  |  ",
                        "Titik merah = Exception (return aktual < −VaR HS)  |  ",
                        "Window = ", window, " hari"),
      x        = NULL,
      y        = "Log Return Harian (%)",
      color    = "Metode VaR"
    ) +
    theme_jurnal +
    theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 8))
}

p8a <- make_var_plot(var_mw,  "MARKOWITZ")
p8b <- make_var_plot(var_sim, "SIM")
p8  <- grid.arrange(p8a, p8b, ncol = 1)

ggsave("08_rolling_var_backtesting.png",
       plot = p8, width = 12, height = 9, dpi = 150, bg = "white")

# ── Plot 9: Perbandingan Rata-rata VaR ──
cat("Membuat Plot 9: Perbandingan rata-rata VaR...\n")
## Membuat Plot 9: Perbandingan rata-rata VaR...
var_cmp <- data.frame(
  Model  = rep(c("Markowitz", "SIM"), each = 2),
  Metode = rep(c("Historical Simulation", "Monte Carlo"), 2),
  VaR    = c(mean(var_mw$HS),  mean(var_mw$MC),
             mean(var_sim$HS), mean(var_sim$MC)) * 100
)

p9 <- ggplot(var_cmp, aes(x = Metode, y = VaR, fill = Model)) +
  geom_col(position = "dodge", width = 0.6,
           color = "white", linewidth = 0.3) +
  geom_text(aes(label = paste0(round(VaR, 4), "%")),
            position = position_dodge(0.6),
            vjust = -0.5, size = 3.2) +
  scale_fill_manual(values = warna_model) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(
    title    = "Perbandingan Rata-rata VaR 95% Harian",
    subtitle = "Markowitz vs SIM  |  Historical Simulation & Monte Carlo",
    x        = "Metode Estimasi VaR",
    y        = "Rata-rata VaR (%)",
    fill     = "Model Portofolio",
    caption  = paste0("Rolling window = ", window,
                      " hari  |  Confidence level = 95%")
  ) +
  theme_jurnal

ggsave("09_perbandingan_ratarata_var.png",
       plot = p9, width = 8, height = 6, dpi = 150, bg = "white")

# ── Plot 10: Perbandingan Kinerja (Return, Risiko, Sharpe) ──
cat("Membuat Plot 10: Perbandingan kinerja portofolio...\n")
## Membuat Plot 10: Perbandingan kinerja portofolio...
kinerja_df <- data.frame(
  Model     = rep(c("Markowitz", "SIM"), each = 3),
  Indikator = rep(c("Expected Return\n(% per hari)",
                    "Std Deviation\n(% per hari)",
                    "Sharpe Ratio"), 2),
  Nilai     = c(ret_mw * 100, sd_mw * 100, sr_mw,
                ret_sim * 100, sd_sim * 100, sr_sim)
)

p10 <- ggplot(kinerja_df, aes(x = Indikator, y = Nilai, fill = Model)) +
  geom_col(position = "dodge", width = 0.6,
           color = "white", linewidth = 0.3) +
  geom_text(aes(label = round(Nilai, 4)),
            position = position_dodge(0.6),
            vjust = -0.5, size = 3) +
  scale_fill_manual(values = warna_model) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.2))) +
  labs(
    title    = "Perbandingan Kinerja Portofolio: Markowitz vs SIM",
    subtitle = "Expected Return, Risiko (Std Deviation), dan Sharpe Ratio",
    x        = NULL,
    y        = "Nilai",
    fill     = "Model Portofolio"
  ) +
  theme_jurnal

ggsave("10_perbandingan_kinerja_portofolio.png",
       plot = p10, width = 9, height = 6, dpi = 150, bg = "white")