## ===========================================
## 0) PACKAGES, THEME, SEED
## ===========================================
# Blok suppressPackageStartupMessages untuk menekan pesan loading package
suppressPackageStartupMessages({
# Mendefinisikan vektor berisi nama-nama package yang dibutuhkan
pkgs <- c(
"quantmod", # Untuk download data finansial (getSymbols)
"dplyr", # Manipulasi data (filter, mutate, select, dll)
"tidyr", # Mengatur struktur data (pivot, gather, dll)
"moments", # Menghitung skewness dan kurtosis (untuk VaR Ekspansi Cornish-Fisher)
"ggplot2", # Visualisasi data (scatter plot, line chart, dll)
"ggrepel", # Label yang tidak overlap pada plot
"plotly", # Visualisasi interaktif
"heatmaply", # Heatmap interaktif (untuk correlation matrix)
"gt", # Membuat tabel yang cantik untuk laporan
"scales", # Format angka dan skala untuk visualisasi
"viridis", # Palet warna untuk visualisasi
"RColorBrewer", # Palet warna untuk visualisasi
"zoo", # Time series object dan fungsi
"xts", # Extensible time series (digunakan quantmod)
"stats", # Fungsi statistik dasar (cov, cor, dll)
"openxlsx", # Export hasil ke excel
"DT" # Interactive DataTable untuk R Markdown
)
# Mengecek package mana yang belum terinstall
# sapply: apply function ke setiap elemen pkgs
# requireNamespace: cek apakah package tersedia (TRUE/FALSE)
# !sapply: negasi, ambil yang FALSE (belum terinstall)
inst <- pkgs[!sapply(pkgs, requireNamespace, quietly = TRUE)]
# Jika ada package yang belum terinstall, install otomatis
if (length(inst)) install.packages(inst)
# Load semua package ke R session
# invisible: suppress output dari lapply
# lapply: apply library() ke setiap package
invisible(lapply(pkgs, library, character.only = TRUE))
})
## Warning: package 'quantmod' was built under R version 4.5.1
## Warning: package 'plotly' was built under R version 4.5.1
## Warning: package 'heatmaply' was built under R version 4.5.1
## Warning: package 'gt' was built under R version 4.5.1
# Set seed untuk reproducibility
# Semua proses random akan menghasilkan hasil yang sama setiap kali dijalankan
set.seed(42)
# Set theme default untuk semua plot ggplot2
# theme_minimal: theme bersih tanpa grid background berlebihan
# base_size = 12: ukuran font dasar 12pt
theme_set(theme_minimal(base_size = 12))
## =============================================
## 1) KONFIGURASI DATA & PARAMETER GLOBAL
## =============================================
## ---------------------------------------------
## A. PARAMETER DATA & PERIODE
## ---------------------------------------------
# Jumlah hari perdagangan dalam setahun untuk annualisasi
trading_days <- 252
# Daftar 17 saham yang akan dianalisis (Drop saham PGEO karena banyak missing data)
# Berasal dari indeks SMInfra18 evaluasi mayor April 2025
tickers <- c("AKRA","BBNI","BBRI","BMRI","EXCL","INTP",
"ISAT","JSMR","MEDC","MTEL","PGAS",
"PTPP","SMGR","SSIA","TLKM","TOWR","UNTR")
# Menambahkan suffix ".JK" untuk Jakarta Stock Exchange
# Diperlukan untuk download data dari Yahoo Finance via quantmod
yahoo_tickers <- paste0(tickers, ".JK")
# Periode data historis untuk analisis
from_date <- as.Date("2022-01-01") # Tanggal mulai: 1 Januari 2022
to_date <- as.Date("2025-09-30") # Tanggal akhir: 30 September 2025
## ---------------------------------------------
## B. CONSTRAINT PORTOFOLIO
## ---------------------------------------------
# CONSTRAINT 1: Batas bobot individual
alpha_w <- 0.10 # Bobot minimum: 10% (jika saham masuk portofolio)
beta_w <- 0.40 # Bobot maksimum: 40% (mencegah konsentrasi berlebihan)
# wi = 0 atau α ≤ wi ≤ β
# CONSTRAINT 2: Cardinality constraint
d_min <- 3 # Minimum 3 saham dalam portofolio (diversifikasi minimum)
d_max <- 5 # Maksimum 5 saham dalam portofolio (tidak terlalu kompleks)
# dmin ≤ dnz ≤ dmax, dimana dnz = jumlah bobot non-zero
# CONSTRAINT 3: Must-have asset
must_have_name <- "TLKM" # TLKM wajib ada di setiap portofolio
# IMPLICIT CONSTRAINT:
# Σwi = 1 (total bobot = 100%)
# wi ≥ 0 (no short selling)
## ---------------------------------------------
## C. PARAMETER ALGORITMA MOCv-ABC
## ---------------------------------------------
# PARAMETER 1: Scout Bee Limit
limit_factor <- 10 # Multiplier untuk dimensi
# limit = limit_factor * d
# Jika solusi tidak improve dalam 'limit' iterasi,
# solusi akan di-reinisialisasi (Scout Bee Phase)
# PARAMETER 2: Sigma untuk Onlooker Bee
sigma_abc <- 0.5 # Scaling factor untuk covariance-guided exploration
# Zi = m + σBDr
# Nilai lebih besar = eksplorasi lebih luas
# PARAMETER 3: Maximum Iterations
max_iter <- 200 # Jumlah maksimum iterasi
# PARAMETER 4: Debug flag
debug_opt <- TRUE # TRUE = print progress setiap iterasi
# FALSE = hanya hasil akhir
# Berguna untuk monitoring konvergensi
## ---------------------------------------------
## D. PARAMETER EVALUASI KINERJA PORTOFOLIO
## ---------------------------------------------
# PARAMETER 1: Risk-Free Rate untuk Sharpe Ratio
rf <- 0.0475 # 4.75% per tahun (BI Rate)
# Sharpe Ratio = (Rp - Rf) / σp
# Digunakan untuk memilih portofolio terbaik
# PARAMETER 2: VaR Cornish-Fisher
cl <- 0.95 # Confidence Level = 95%
# Artinya kita 95% yakin kerugian yang dialami tidak akan melebihi VaR
alpha <- 1 - cl # Alpha = 5% (tail risk yang diukur)
V0 <- 1 # Nilai investasi awal
hp <- 1 # Holding period = 1 hari
## =============================================
## 2) AMBIL DATA SAHAM & BENTUK PANEL
## =============================================
## -------------------------------------------
## A. DOWNLOAD DATA HARGA SAHAM
## -------------------------------------------
# Loop setiap ticker untuk download data dari Yahoo Finance
prices_list <- lapply(yahoo_tickers, function(sym) {
# suppressWarnings: menekan warning dari getSymbols (misal: data tidak lengkap)
# getSymbols: fungsi dari quantmod untuk download data OHLCV
# - sym: ticker symbol (misal: "AKRA.JK")
# - from/to: periode data
# Output: xts object dengan kolom OHLCV (Open, High, Low, Close, Volume, Adjusted)
suppressWarnings(
getSymbols(sym, from = from_date, to = to_date, auto.assign = FALSE)
)
})
## -------------------------------------------
## B. EKSTRAKSI ADJUSTED CLOSE PRICE
## -------------------------------------------
# Ad() adalah fungsi dari quantmod untuk extract kolom "Adjusted" price
# Adjusted price sudah memperhitungkan stock split dan dividend
adj_list <- lapply(prices_list, Ad)
## -------------------------------------------
## C. MERGE SEMUA SAHAM MENJADI 1 XTS OBJECT
## -------------------------------------------
# do.call(merge, list): gabungkan semua xts dalam list menjadi 1 xts
# Hasilnya xts dengan 18 kolom (1 untuk setiap saham)
prices_xts <- do.call(merge, adj_list)
# Ganti nama kolom
colnames(prices_xts) <- tickers
# na.omit: buang baris yang ada missing value (NA)
prices_xts <- na.omit(prices_xts)
## -------------------------------------------
## D. HITUNG LOG-RETURN HARIAN
## -------------------------------------------
# Log-return: r_t = ln(P_t / P_{t-1}) = ln(P_t) - ln(P_{t-1})
rets_xts <- na.omit(diff(log(prices_xts)))
## -------------------------------------------
## E. DATA QUALITY CHECK: HAPUS SAHAM TANPA DATA
## -------------------------------------------
# Cek apakah ada kolom (saham) yang semua nilai-nya NA
# Bisa terjadi jika data download gagal atau saham delisting
na_all <- sapply(rets_xts, function(x) all(is.na(x)))
# Jika ada saham dengan all-NA, hapus dari dataset
if (any(na_all)) {
keep <- which(!na_all) # Index kolom yang valid
rets_xts <- rets_xts[, keep, drop = FALSE] # Keep kolom valid
tickers <- tickers[keep] # Update daftar ticker
# Print pesan notifikasi
message("Menghapus aset tanpa data lengkap. Sisa aset:", length(tickers), " -> ", paste(tickers, collapse = ", "))
}
## -------------------------------------------
## F. TRANSFORMASI KE LONG FORMAT (untuk EDA & Plotting)
## -------------------------------------------
# 📊 CONVERT RETURNS KE LONG FORMAT
# Wide format (xts): date | AKRA | BBNI | BBRI | ...
# Long format (df): date | ticker | ret
rets_df <- rets_xts |>
as.data.frame() |> # xts -> data.frame
tibble::rownames_to_column("date") |> # Date jadi kolom (dari rownames)
mutate(date = as.Date(date)) |> # Convert ke Date object
pivot_longer( # Wide -> Long transformation
-date, # Keep kolom date
names_to = "ticker", # Nama kolom -> "ticker"
values_to = "ret" # Value -> "ret"
)
# Hasil: tibble dengan kolom: date | ticker | ret
# Total rows = n_dates × n_tickers
# 📊 CONVERT PRICES KE LONG FORMAT (analog dengan returns)
prices_df <- prices_xts |>
as.data.frame() |>
tibble::rownames_to_column("date") |>
mutate(date = as.Date(date)) |>
pivot_longer(-date, names_to = "ticker", values_to = "price")
# Hasil: tibble dengan kolom: date | ticker | price
## -------------------------------------------
## G. EXPORT DATA KE EXCEL (untuk Dokumentasi & Backup)
## -------------------------------------------
# Nama file output
output_file_prices_returns <- "data_harga_dan_return.xlsx"
# Buat workbook Excel baru (objek untuk menulis Excel)
wb <- createWorkbook()
## ======================
## 📋 SHEET 1: HARGA SAHAM (Wide Format)
## ======================
# Tambah worksheet bernama "Harga_Saham"
addWorksheet(wb, "Harga_Saham")
# Transform ke wide format untuk Excel (lebih mudah dibaca)
# Long: date | ticker | price → Wide: date | AKRA | BBNI | BBRI | ...
prices_export <- prices_df |>
tidyr::pivot_wider(
names_from = ticker, # Kolom ticker jadi nama kolom baru
values_from = price # Value price jadi isi cell
) |>
dplyr::arrange(date) # Sort berdasarkan tanggal
# Tulis data ke Excel dengan format table
writeDataTable(
wb,
sheet = "Harga_Saham", # Nama sheet
x = prices_export, # Data yang ditulis
tableStyle = "TableStyleLight9", # Style Excel table (auto-filter, stripes)
withFilter = TRUE # Enable filter di header
)
## ======================
## 📋 SHEET 2: RETURN SAHAM (Wide Format)
## ======================
addWorksheet(wb, "Return_Saham")
# Transform return ke wide format
rets_export <- rets_df |>
tidyr::pivot_wider(names_from = ticker, values_from = ret) |>
dplyr::arrange(date)
# Tulis data return ke Excel
writeDataTable(
wb,
sheet = "Return_Saham",
x = rets_export,
tableStyle = "TableStyleLight9",
withFilter = TRUE
)
## ======================
## 💾 SIMPAN WORKBOOK KE FILE
## ======================
saveWorkbook(wb, file = output_file_prices_returns, overwrite = TRUE)
# Print konfirmasi dengan full path
cat("✅ File Excel berhasil disimpan ke:\n",
normalizePath(output_file_prices_returns), "\n")
## ✅ File Excel berhasil disimpan ke:
## D:\KULIAH\BISMILLAH SKRIPSI\data_harga_dan_return.xlsx
## ======================
## 👁️ PREVIEW DATA (Opsional)
## ======================
# Tampilkan 5 baris pertama dari masing-masing dataset
# Untuk verifikasi bahwa data sudah benar
print(head(prices_export, 5)) # Preview harga
## # A tibble: 5 × 18
## date AKRA BBNI BBRI BMRI EXCL INTP ISAT JSMR MEDC MTEL PGAS
## <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2022-01-03 656. 2772. 3285. 2767. 2706. 10293. 1403. 3720. 389. 728. 961.
## 2 2022-01-04 632. 2896. 3269. 2816. 2689. 10095. 1365. 3682. 395. 728. 961.
## 3 2022-01-05 619. 2865. 3308. 2757. 2646. 9941. 1343. 3645. 407. 728. 937.
## 4 2022-01-06 626. 2865. 3269. 2757. 2543. 9809. 1376. 3636. 411. 724. 927.
## 5 2022-01-07 645. 2917. 3293. 2767. 2586. 9633. 1376. 3692. 414. 724. 930.
## # ℹ 6 more variables: PTPP <dbl>, SMGR <dbl>, SSIA <dbl>, TLKM <dbl>,
## # TOWR <dbl>, UNTR <dbl>
print(head(rets_export, 5)) # Preview return
## # A tibble: 5 × 18
## date AKRA BBNI BBRI BMRI EXCL INTP ISAT
## <date> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2022-01-04 -0.0374 0.0436 -0.00480 0.0176 -0.00635 -0.0194 -0.0276
## 2 2022-01-05 -0.0192 -0.0107 0.0119 -0.0211 -0.0161 -0.0154 -0.0161
## 3 2022-01-06 0.00966 0 -0.0119 0 -0.0396 -0.0134 0.0241
## 4 2022-01-07 0.0308 0.0178 0.00719 0.00355 0.0167 -0.0181 0
## 5 2022-01-10 0.00233 -0.00354 -0.00239 0 -0.00664 -0.00229 -0.0160
## # ℹ 10 more variables: JSMR <dbl>, MEDC <dbl>, MTEL <dbl>, PGAS <dbl>,
## # PTPP <dbl>, SMGR <dbl>, SSIA <dbl>, TLKM <dbl>, TOWR <dbl>, UNTR <dbl>
# Cek periode data yang benar-benar tersedia
actual_start <- start(prices_xts)
actual_end <- end(prices_xts)
n_obs <- nrow(prices_xts)
cat("\n📅 RINGKASAN DATA HARGA:\n")
##
## 📅 RINGKASAN DATA HARGA:
cat(" Periode diminta : ", from_date, "s.d.", to_date, "\n")
## Periode diminta : 18993 s.d. 20361
cat(" Periode tersedia: ", actual_start, "s.d.", actual_end, "\n")
## Periode tersedia: 18995 s.d. 20360
cat(" Jumlah observasi: ", n_obs, "hari\n")
## Jumlah observasi: 894 hari
cat(" Jumlah saham : ", ncol(prices_xts), "saham\n")
## Jumlah saham : 17 saham
cat(" Missing data : ", sum(is.na(prices_xts)), "cells\n\n")
## Missing data : 0 cells
# Cek periode data yang benar-benar tersedia
rets_start <- start(rets_xts)
rets_end <- end(rets_xts)
rets_obs <- nrow(rets_xts)
cat("\n📅 RINGKASAN DATA RETURN:\n")
##
## 📅 RINGKASAN DATA RETURN:
cat(" Periode diminta : ", from_date, "s.d.", to_date, "\n")
## Periode diminta : 18993 s.d. 20361
cat(" Periode tersedia: ", rets_start, "s.d.", rets_end, "\n")
## Periode tersedia: 18996 s.d. 20360
cat(" Jumlah observasi: ", rets_obs, "hari\n")
## Jumlah observasi: 893 hari
cat(" Jumlah saham : ", ncol(rets_xts), "saham\n")
## Jumlah saham : 17 saham
cat(" Missing data : ", sum(is.na(rets_xts)), "cells\n\n")
## Missing data : 0 cells
## =======================================================
## 3) EDA INTERAKTIF (Harga, Distribusi, Deskriptif, Korelasi, Rolling Vol)
## =======================================================
## -------------------------------------------
## Helper Function: Membagi Saham ke Dalam Grup Visualisasi
## -------------------------------------------
# Fungsi pembantu untuk membagi daftar ticker saham menjadi grup-grup kecil
# Tujuannya agar visualisasi tidak terlalu padat dan tetap mudah dibaca
make_groups <- function(x, size = 9) {
# x : vektor nama ticker saham (misal: c("AKRA", "BBNI", ...))
# size : jumlah maksimal saham per grup (default 9)
# Membuat indeks grup untuk setiap ticker:
# - seq_len(ceiling(length(x) / size)): membuat urutan grup (1, 2, 3, ...)
# - each = size: setiap grup diulang sebanyak 'size' kali
# - length.out = length(x): memastikan panjang indeks sama dengan jumlah ticker
group_indices <- rep(seq_len(ceiling(length(x) / size)), each = size, length.out = length(x))
# split(x, group_indices): membagi vektor ticker menjadi list per grup
split(x, group_indices)
}
# Membagi ticker menjadi grup-grup berisi maksimal 6 saham
groups <- make_groups(tickers, size = 6)
## -------------------------------------------
## 3a). Visualisasi Harga Saham
## -------------------------------------------
# Fungsi utama untuk membuat plot harga saham interaktif per grup
plot_prices_by_group <- function(g_ids = 1:length(groups)) {
# g_ids: vektor indeks grup yang ingin divisualisasikan (default: semua grup)
# Inisialisasi list kosong untuk menyimpan plot hasil
pls <- vector("list", length(g_ids))
# Loop untuk setiap grup yang dipilih
for (ii in seq_along(g_ids)) {
# Mengambil daftar ticker untuk grup ke-ii
g <- groups[[ g_ids[ii] ]]
# Filter data harga hanya untuk ticker dalam grup ini
df_sub <- prices_df |> dplyr::filter(ticker %in% g)
# Buat palet viridis sebanyak jumlah ticker di grup
pal <- viridis(length(g), option = "D")
names(pal) <- g
# Membuat plot time series harga saham dengan ggplot2
p <- ggplot(df_sub, aes(date, price, color = ticker)) +
# Garis utama harga saham per ticker
geom_line(linewidth = 0.6, alpha = 0.9, show.legend = FALSE) +
# Membagi plot menjadi panel per ticker (facet)
facet_wrap(~ ticker, ncol = 3, scales = "free_y") +
# Menambahkan garis trend untuk tiap saham
geom_smooth(method = "loess", se = FALSE, linewidth = 0.2, color = "black") +
scale_color_manual(values = pal) +
# Judul dan label sumbu
labs(title = "Time Series Harga Saham", x = NULL, y = "Harga") +
# Kustomisasi tema:
theme(axis.text.y = element_text(size = 7),
axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
plot.title = element_text(size = 16, face = "bold", hjust = 0.5))
# Mengubah plot ggplot2 menjadi plotly
# Fitur: hover, zoom, pan, export, dsb
pls[[ii]] <- plotly::ggplotly(p)
}
# Mengembalikan list plot interaktif untuk setiap grup
pls
}
# Membuat plot interaktif untuk grup 1-3
price_plots <- plot_prices_by_group(1:3)
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
# Menampilkan plot untuk grup pertama
price_plots[[1]]
# Menampilkan plot untuk grup kedua
price_plots[[2]]
price_plots[[3]]
## -------------------------------------------
## 3b). Visualisasi Distribusi Return
## -------------------------------------------
# Fungsi untuk membuat plot distribusi return (density) per grup saham
plot_density_by_group <- function(g_ids = 1:length(groups)) {
pls <- vector("list", length(g_ids)) # list kosong untuk plot per grup
# Hitung statistik deskriptif per ticker (mean, sd, quantile)
summ_by_tkr_all <- rets_df |>
group_by(ticker) |>
summarise(
mean_d = mean(ret, na.rm = TRUE), # Rata-rata return harian
sd_d = sd(ret, na.rm = TRUE), # Standar deviasi return harian
x_min = quantile(ret, 0.005, na.rm = TRUE), # Batas bawah (hindari outlier ekstrem)
x_max = quantile(ret, 0.995, na.rm = TRUE), # Batas atas
.groups = "drop"
)
# Buat data kurva normal teoritis per ticker (untuk overlay)
norm_df_all <- summ_by_tkr_all |>
rowwise() |>
mutate(x = list(seq(x_min, x_max, length.out = 400)), # Titik x untuk kurva normal
y = list(dnorm(x, mean = mean_d, sd = sd_d)) # Density normal
) |>
tidyr::unnest(c(x, y)) |> ungroup()
for (ii in seq_along(g_ids)) {
tick_sub <- groups[[ g_ids[ii] ]] # Ticker dalam grup ke-ii
rets_sub <- rets_df |> filter(ticker %in% tick_sub) # Data return grup ini
summ_sub <- summ_by_tkr_all |> filter(ticker %in% tick_sub) # Statistik grup
norm_sub <- norm_df_all |> filter(ticker %in% tick_sub) # Kurva normal grup
# Buat palet viridis
pal <- viridis(length(tick_sub), option = "D")
names(pal) <- tick_sub
# Plot density empiris dan overlay kurva normal
p <- ggplot() +
geom_density(data = rets_sub, aes(ret, after_stat(density), color = ticker),
linewidth = 0.8, alpha = 0.9, show.legend = FALSE) + # Density empiris
geom_line(data = norm_sub, aes(x, y), color = "black", linewidth = 0.4) + # Kurva normal
geom_vline(xintercept = 0, color = "grey40", linewidth = 0.4) + # Garis di nol
geom_vline(data = summ_sub, aes(xintercept = mean_d),
color = "black", linewidth = 0.5, linetype = "dashed") + # Garis mean
facet_wrap(~ ticker, ncol = 3, scales = "free_y") +
scale_color_manual(values = pal) +
labs(title = "Distribusi Return Harian", x = "Return", y = "Kepadatan") +
theme(plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
strip.text = element_text(face = "bold"),
panel.grid.minor = element_blank(),
legend.position = "none")
pls[[ii]] <- plotly::ggplotly(p)
}
pls
}
density_plots <- plot_density_by_group(1:3)
density_plots[[1]]
density_plots[[2]]
density_plots[[3]]
## -------------------------------------------
## 3c). Statistik Deskriptif & Uji Normalitas
## -------------------------------------------
# -----------------------------------------------------------------------------
# Fungsi custom untuk uji normalitas Jarque-Bera pada vektor return
jb_test <- function(x) {
x <- x[is.finite(x)] # Hapus nilai NA, Inf, -Inf agar perhitungan valid
n <- length(x) # Hitung jumlah data setelah pembersihan
if (n < 8) return(NA_real_) # Jika data kurang dari 8, return NA (JB tidak valid)
s <- moments::skewness(x) # Hitung skewness
k <- moments::kurtosis(x) # Hitung kurtosis
JB <- n/6 * (s^2 + ((k-3)^2 / 4)) # Hitung statistik JB sesuai rumus standar
1 - pchisq(JB, df = 2) # Hitung p-value dari distribusi chi-square (df = 2)
}
# -----------------------------------------------------------------------------
# Hitung statistik deskriptif per saham
stats_df <- rets_df |>
group_by(ticker) |>
summarise(
mean_d = mean(ret, na.rm = TRUE), # Rata-rata return harian
var_d = var(ret, na.rm = TRUE), # Variansi harian
sd_d = sd(ret, na.rm = TRUE), # Standar deviasi harian (volatilitas)
skew = moments::skewness(ret, na.rm = TRUE), # Skewness harian
kurt = moments::kurtosis(ret, na.rm = TRUE), # Kurtosis harian
jb_p = jb_test(ret), # p-value uji normalitas JB
.groups = "drop" # Drop grouping untuk operasi selanjutnya
) |>
mutate(
mean_ann = mean_d * trading_days, # Annualisasi mean
sd_ann = sd_d * sqrt(trading_days), # Annualisasi SD
var_ann = sd_ann^2, # Varians tahunan
rank_mean = rank(-mean_ann, ties.method = "min"), # Ranking return (semakin tinggi semakin baik)
rank_sd = rank(sd_ann, ties.method = "min") # Ranking risiko (semakin rendah semakin baik)
) |>
arrange(rank_sd, rank_mean) # Urutkan dari risk terendah, lalu return tinggi
# -----------------------------------------------------------------------------
# Tampilkan hasil statistik deskriptif
stats_gt <- stats_df |>
select(ticker, mean_ann, sd_ann, skew, kurt, jb_p, rank_mean, rank_sd) |> # Pilih kolom penting
gt(rowname_col = "ticker") |> # Ticker sebagai rowname
fmt_number(columns = c(mean_ann, sd_ann), decimals = 4) |> # Format mean & SD 4 desimal
fmt_number(columns = c(skew, kurt), decimals = 2) |> # Format skew & kurt 2 desimal
fmt_number(columns = c(jb_p, rank_mean, rank_sd), decimals = 2) |> # Format p-value & ranking 2 desimal
tab_header(title = md("**Statistik Deskriptif Return**")) |> # Judul tabel
cols_label(mean_ann = "Mean (ann.)", sd_ann="SD (ann.)", skew="Skew", kurt="Kurt",
jb_p="JB p-val", rank_mean="Rank Mean", rank_sd="Rank SD") |> # Label kolom
data_color(columns = c(rank_sd), # color coding untuk risk ranking
colors = scales::col_numeric(c("#2DC937","#E7B416","#CC3232"), domain = NULL)) |>
tab_source_note(md("Warna hijau = lebih baik (SD rendah). JB p-val kecil → indikasi non-normal."))
## Warning: Since gt v0.9.0, the `colors` argument has been deprecated.
## • Please use the `fn` argument instead.
## This warning is displayed once every 8 hours.
stats_gt
| Statistik Deskriptif Return |
|
Mean (ann.) |
SD (ann.) |
Skew |
Kurt |
JB p-val |
Rank Mean |
Rank SD |
| TLKM |
−0.0263 |
0.2756 |
−0.08 |
4.65 |
0.00 |
12.00 |
1.00 |
|---|
| MTEL |
−0.0641 |
0.2778 |
−0.32 |
10.13 |
0.00 |
13.00 |
2.00 |
|---|
| BBRI |
0.0542 |
0.2871 |
−0.03 |
6.15 |
0.00 |
9.00 |
3.00 |
|---|
| BBNI |
0.1152 |
0.2986 |
0.09 |
5.28 |
0.00 |
7.00 |
4.00 |
|---|
| BMRI |
0.1309 |
0.2995 |
−0.13 |
6.41 |
0.00 |
6.00 |
5.00 |
|---|
| JSMR |
−0.0124 |
0.3076 |
0.37 |
4.90 |
0.00 |
10.00 |
6.00 |
|---|
| PGAS |
0.1601 |
0.3080 |
0.31 |
5.80 |
0.00 |
5.00 |
7.00 |
|---|
| UNTR |
0.1901 |
0.3281 |
−0.08 |
11.02 |
0.00 |
3.00 |
8.00 |
|---|
| INTP |
−0.1212 |
0.3294 |
0.18 |
7.85 |
0.00 |
14.00 |
9.00 |
|---|
| EXCL |
−0.0157 |
0.3368 |
0.32 |
5.91 |
0.00 |
11.00 |
10.00 |
|---|
| TOWR |
−0.1568 |
0.3484 |
0.75 |
8.35 |
0.00 |
15.00 |
11.00 |
|---|
| SMGR |
−0.2255 |
0.3744 |
0.16 |
8.30 |
0.00 |
16.00 |
12.00 |
|---|
| AKRA |
0.1741 |
0.4009 |
0.54 |
9.01 |
0.00 |
4.00 |
13.00 |
|---|
| ISAT |
0.0672 |
0.4188 |
−0.11 |
8.89 |
0.00 |
8.00 |
14.00 |
|---|
| MEDC |
0.3524 |
0.4963 |
0.63 |
6.08 |
0.00 |
2.00 |
15.00 |
|---|
| PTPP |
−0.2600 |
0.5159 |
1.50 |
11.56 |
0.00 |
17.00 |
16.00 |
|---|
| SSIA |
0.3688 |
0.5521 |
0.93 |
10.45 |
0.00 |
1.00 |
17.00 |
|---|
| Warna hijau = lebih baik (SD rendah). JB p-val kecil → indikasi non-normal. |
## 3c.1) Simpan Expected Return & Risiko (SD) tiap saham ke Excel
# Pilih kolom ticker, mean_ann (expected return tahunan), dan sd_ann (risiko tahunan) dari stats_df
risk_return_df <- stats_df |>
dplyr::select(ticker, mean_ann, sd_ann) |>
dplyr::rename(Expected_Return = mean_ann, Risk_SD = sd_ann)
# - stats_df: Data frame berisi statistik deskriptif tiap saham
# - select(): Ambil hanya kolom yang dibutuhkan
# - rename(): Ubah nama kolom
# Tentukan nama file output Excel
output_file <- "saham_risk_return.xlsx"
# Ekspor data risk-return ke file Excel
openxlsx::write.xlsx(risk_return_df, file = output_file,
sheetName = "RiskReturn", overwrite = TRUE)
# - openxlsx::write.xlsx: Fungsi ekspor ke Excel
# - overwrite = TRUE: File lama akan ditimpa jika sudah ada
# Tampilkan konfirmasi lokasi file yang berhasil disimpan
cat("✅ File disimpan ke:", normalizePath(output_file), "\n")
## ✅ File disimpan ke: D:\KULIAH\BISMILLAH SKRIPSI\saham_risk_return.xlsx
## 3c.2) Top-5 Expected Return & Top-5 Risiko Terendah
# Ambil 5 saham dengan expected return tahunan tertinggi
top5_mean_ann <- stats_df |>
slice_max(order_by = mean_ann, n = 5) |>
arrange(mean_ann)
# - slice_max(): Pilih 5 baris dengan mean_ann terbesar
# - arrange(): Urutkan dari return terendah ke tertinggi
# Buat bar chart horizontal untuk top-5 expected return
p_top_mean_ann <- ggplot(top5_mean_ann, aes(reorder(ticker, mean_ann), mean_ann)) +
geom_col(fill = "#2E86DE") + coord_flip() +
geom_text(aes(label = round(mean_ann, 3)),
hjust = -0.1, size = 3.2) + # Label angka di ujung bar
expand_limits(y = max(top5_mean_ann$mean_ann) * 1.15) + # Tambah ruang di atas bar
labs(title = "Top-5 Expected Return Tertinggi", x = NULL, y = "Mean")
plotly::ggplotly(p_top_mean_ann)
# Ambil 5 saham dengan risiko tahunan (SD) terendah
top5_low_sd_ann <- stats_df |>
slice_min(order_by = sd_ann, n = 5) |>
arrange(sd_ann)
# Buat bar char horizontal untuk top-5 risiko terendah
p_low_sd_ann <- ggplot(top5_low_sd_ann, aes(reorder(ticker, -sd_ann), sd_ann)) +
geom_col(fill = "#27AE60") + coord_flip() +
geom_text(aes(label = round(sd_ann, 4)), hjust = -0.1, size = 3.2) +
expand_limits(y = max(top5_low_sd_ann$sd_ann) * 1.15) +
labs(title = "Top-5 Risiko Terendah", x = NULL, y = "SD")
plotly::ggplotly(p_low_sd_ann)
## -------------------------------------------
## 3d). Heatmap Matriks Varian-Kovarian & Korelasi
## -------------------------------------------
## 3d.1) Heatmap Matriks Varian-Kovarian
# 1. Transformasi data return ke format wide
wide_ret <- rets_df |>
tidyr::pivot_wider(names_from = ticker, values_from = ret) |> # Ubah long ke wide
dplyr::select(-date) # Hilangkan kolom tanggal
# 2. Hitung matriks kovarians harian antar saham
cov_daily <- stats::cov(wide_ret, use = "pairwise.complete.obs", method = "pearson")
# - 'pariwise.complete.obs': gunakan semua data yang tersedia untuk setiap pasangan saham
# - 'pearson': metode standar untuk kovarians
# 3. Hitung matriks korelasi harian antar saham
cor_mat <- stats::cor(wide_ret, use = "pairwise.complete.obs", method = "pearson")
# 4. Annualisasi matriks kovarians
cov_ann <- cov_daily * trading_days
# 5. Tampilkan matriks kovarians tahunan dalam tabel interaktif
datatable(
round(cov_ann, 4), # Matriks kovarians dibulatkan 4 desimal
extensions = "Buttons", # Fitur export
options = list(
dom = "Bfrtip", # Layout tombol export
buttons = c("copy", "csv", "excel"),
scrollX = TRUE, # Scroll horizontal jika tabel lebar
pageLength = 10 # Jumlah baris per halaman
),
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left;",
"Tabel Varian–Kovarian (Annualized)"
)
)
datatable(
round(cor_mat, 4), # Matriks kovarians dibulatkan 4 desimal
extensions = "Buttons", # Fitur export
options = list(
dom = "Bfrtip", # Layout tombol export
buttons = c("copy", "csv", "excel"),
scrollX = TRUE, # Scroll horizontal jika tabel lebar
pageLength = 10 # Jumlah baris per halaman
),
caption = htmltools::tags$caption(
style = "caption-side: top; text-align: left;",
"Tabel Korelasi Saham"
)
)
## 3d.2) Heatmap Korelasi + Dendrogram
# -----------------------------------------------------------------------------
# Fungsi untuk menghitung matriks p-value korelasi antar saham
pval_mat <- function(M) {
k <- ncol(M)
P <- matrix(NA_real_, k, k, dimnames = list(colnames(M), colnames(M)))
for (i in 1:k) for (j in i:k) {
ct <- suppressWarnings(cor.test(M[, i], M[, j], use = "pairwise.complete.obs"))
P[i, j] <- P[j, i] <- ct$p.value
}
P
}
# -----------------------------------------------------------------------------
# Hitung matriks korelasi Pearson antar saham (simetris, diagonal = 1)
cor_Pear <- cor(wide_ret, use = "pairwise.complete.obs", method = "pearson")
# Hitung matriks p-value untuk korelasi (opsional, untuk masking signifikansi)
p_Pear <- pval_mat(as.data.frame(wide_ret))
# Pilih matriks korelasi yang akan divisualisasikan (bisa masking p-value jika mau)
cor_sig <- cor_Pear
# -----------------------------------------------------------------------------
# Lookup sektor untuk anotasi warna pada heatmap
sector_lookup <- c(
AKRA="Energi", BBNI="Keuangan", BBRI="Keuangan", BMRI="Keuangan",
EXCL="Infrastruktur", INTP="Barang Baku", ISAT="Infrastruktur",
JSMR="Infrastruktur", MEDC="Energi", MTEL="Infrastruktur",
PGAS="Energi", PGEO="Infrastruktur", PTPP="Infrastruktur",
SMGR="Barang Baku", SSIA="Infrastruktur", TLKM="Infrastruktur",
TOWR="Infrastruktur", UNTR="Perindustrian"
)
side_anno <- data.frame(Sector = sector_lookup[colnames(cor_sig)])
row.names(side_anno) <- colnames(cor_sig)
custom_text <- cor_sig
custom_text[lower.tri(custom_text, diag = FALSE)] <- NA
text_matrix <- matrix("", nrow = nrow(cor_sig), ncol = ncol(cor_sig))
for (i in 1:nrow(cor_sig)) {
for (j in i:ncol(cor_sig)) { # j mulai dari i
text_matrix[i, j] <- sprintf("%.2f", cor_sig[i, j])
}
}
sector_colors <- c(
"Energi" = "#E74C3C", # Merah
"Keuangan" = "#3498DB", # Biru
"Infrastruktur" = "#2ECC71", # Hijau
"Barang Baku" = "#F39C12", # Oranye
"Perindustrian" = "#9B59B6" # Ungu
)
side_anno$SectorColor <- sector_colors[side_anno$Sector]
# -----------------------------------------------------------------------------
# HEATMAPLY: Visualisasi heatmap korelasi dengan dendrogram simetris
cor_heatmap <- heatmaply::heatmaply(
cor_sig,
dendrogram = "both", # KOREKSI: Cluster baris & kolom agar simetris
seriate = "OLO", # Optimal Leaf Ordering untuk urutan cluster
symm = TRUE, # KOREKSI: Enforce simetri matriks
row_dend_right = TRUE, # Dendrogram di kiri (baris)
column_dend_top = TRUE, # Dendrogram di atas (kolom)
limits = c(-1, 1), # Skala warna -1 s/d 1
colors = colorRampPalette(rev(RColorBrewer::brewer.pal(11, "RdBu")))(256),
grid_color = "black", grid_width = 0.0001,
label_names = c("Saham i", "Saham j", "Korelasi"),
main = "Heatmap Korelasi Saham dengan Clustering Hierarkis",
showticklabels = c(TRUE, TRUE), na.value = "grey90",
row_side_colors = side_anno[, "Sector", drop = FALSE], # Anotasi sektor di baris
#col_side_colors = side_anno[, "Sector", drop = FALSE], # Anotasi sektor di kolom
custom_hovertext = text_matrix,
cellnote = text_matrix,
cellnote_textposition = "middle center",
cellnote_size = 6,
cellnote_color = "black",
fontsize_row = 8,
fontsize_col = 8,
plot_method = "plotly"
)
cor_heatmap
## 3e) Rolling volatility 30 hari
roll_sd_xts <- zoo::rollapply(rets_xts, 30, sd, by.column = TRUE, align = "right")
roll_sd_df <- roll_sd_xts |>
as.data.frame() |> tibble::rownames_to_column("date") |>
mutate(date = as.Date(date)) |>
pivot_longer(-date, names_to = "ticker", values_to = "roll_sd")
plot_rollvol_by_group <- function(g_ids = 1:length(groups)) {
pls <- vector("list", length(g_ids))
for (ii in seq_along(g_ids)) {
tick_sub <- groups[[ g_ids[ii] ]]
roll_sub <- roll_sd_df |> filter(ticker %in% tick_sub)
p <- ggplot(roll_sub, aes(date, roll_sd, color = ticker)) +
geom_line(linewidth = 0.4, alpha = 0.9, show.legend = FALSE) +
facet_wrap(~ ticker, ncol = 3, scales = "free_y") +
labs(title = "Rolling 30-hari Volatilitas", x = NULL, y = "SD") +
theme(axis.text.y = element_text(size = 7),
axis.text.x = element_text(angle = 45, hjust = 1, size = 8),
plot.title = element_text(size = 16, face = "bold", hjust = 0.5))
pls[[ii]] <- plotly::ggplotly(p)
}
pls
}
roll_plots <- plot_rollvol_by_group(1:3)
roll_plots[[1]]
roll_plots[[2]]
roll_plots[[3]]
## =============================================
## 4) PARAMETER OPTIMASI (annualized μ & Σ)
## =============================================
mu <- colMeans(rets_xts, na.rm = TRUE) * trading_days
# Ambil rata-rata return harian tiap saham, lalu annualisasi
# mu: vektor expected return tahunan untuk setiap saham
S <- cov(rets_xts, use = "pairwise.complete.obs") * trading_days
# Hitung matriks kovarians harian antar saham, lalu annualisasi
# S: matriks kovarians tahunan (jumlah saham x jumlah saham)
stopifnot(is.numeric(mu), is.matrix(S), length(mu) == ncol(S), ncol(S) == nrow(S))
# Validasi: mu harus vektor numerik, S matriks persegi, dimensi konsisten
d <- length(mu)
# d: jumlah saham (dimensi portofolio)
## =============================================
## 5) UTILITAS CCMV + ABC + NSGA-II
## =============================================
# Proyeksi bobot ke simplex terbatas (sum=1, L<=w<=U)
project_box_simplex <- function(w, L, U, s = 1) {
k <- length(w) # Jumlah aset dalam portofolio
L <- rep_len(L, k); U <- rep_len(U, k) # Pastikan L dan U sepanjang k
w <- pmin(pmax(w, L), U) # Batasi setiap bobot di antara L dan U (box constraint)
S <- sum(w) # Hitung jumlah bobot saat ini
if (abs(S - s) < 1e-12) return(w) # Jika sum(w) sudah sama dengan s (default 1), return w
free <- rep(TRUE, k); iter <- 0 # Inisialisasi aset "bebas" dan counter iterasi
while (abs(S - s) > 1e-12 && iter < 1000) { # Iterasi hingga sum(w) ≈ s atau max iterasi
iter <- iter + 1
idx <- which(free); if (!length(idx)) break # Indeks aset yang masih "bebas"
delta <- (S - s) / length(idx) # Koreksi rata-rata untuk aset bebas
w[idx] <- w[idx] - delta # Koreksi bobot agar sum(w) mendekati s
hitL <- which(w < L + 1e-15); hitU <- which(w > U - 1e-15) # Cek yang melewati batas
if (length(hitL)) { w[hitL] <- L[hitL]; free[hitL] <- FALSE } # Kunci di L jika < L
if (length(hitU)) { w[hitU] <- U[hitU]; free[hitU] <- FALSE } # Kunci di U jika > U
S <- sum(w); if (!any(free)) break # Update sum dan cek jika semua sudah fix
}
# Koreksi minor jika sum(w) masih belum pas
if (abs(sum(w) - s) > 1e-8) {
idx <- which((w > L + 1e-12) & (w < U - 1e-12))
if (length(idx) > 0) {
adj <- (sum(w) - s) / length(idx)
w[idx] <- pmin(pmax(w[idx] - adj, L[idx]), U[idx])
}
}
w # Return bobot yang sudah feasible
}
# Bersihkan solusi: tegakkan kardinalitas & aset wajib
clean_portfolio <- function(w, d, alpha, beta, d_min, d_max, must_have_idx_or_NULL) {
stopifnot(length(w) == d) # Validasi panjang bobot = jumlah aset
w[!is.finite(w)] <- 0; w[w < 1e-12] <- 0 # Set bobot tidak valid atau sangat kecil ke nol
active <- which(w > 0) # Indeks aset aktif (bobot > 0)
if (!is.null(must_have_idx_or_NULL)) { # Jika ada aset wajib
if (!(must_have_idx_or_NULL %in% active)) active <- c(must_have_idx_or_NULL, active)
}
nz <- length(active) # Jumlah aset aktif
if (nz < d_min) { # Jika kurang dari minimum
cand <- setdiff(seq_len(d), active)
add <- if (length(cand) >= (d_min - nz)) sample(cand, d_min - nz) else cand
active <- c(active, add)
} else if (nz > d_max) { # Jika lebih dari maksimum
drop_cand <- if (is.null(must_have_idx_or_NULL)) active else setdiff(active, must_have_idx_or_NULL)
if (length(drop_cand) > 0) {
drop <- sample(drop_cand, nz - d_max) # acak agar ragam pola meningkat
active <- setdiff(active, drop)
}
}
w[-active] <- 0 # Set bobot non-aktif ke nol
if (length(active) > 0) {
w_active <- w[active]; w_active[!is.finite(w_active)] <- 0
if (sum(w_active) <= 0) w_active <- rep(1/length(active), length(active))
w_active <- project_box_simplex(w_active, L = alpha, U = beta, s = 1)
w[active] <- w_active
}
w[abs(w) < 1e-12] <- 0
w
}
# Tujuan multi: f1=min SD, f2=min(-Return) → max Return
obj_vec <- function(w, mu, S) {
risk_sd <- sqrt(max(0, as.numeric(t(w) %*% S %*% w))) # Hitung risiko portofolio (SD)
ret <- sum(w * mu) # Hitung expected return portofolio
c(f1 = risk_sd, f2 = -ret) # f1: minimasi risiko; f2; minimasi negatif return
}
# Pareto rules + sorting + crowding
dominates <- function(a, b)
all(a <= b) && any(a < b)
# Cek apakah solusi a mendominasi solusi b:
# - Semua tujuan a <= b (tidak lebih buruk di semua tujuan)
# - Ada minimal satu tujuan a < b (lebih baik di setidaknya satu tujuan)
# Ini adalah definisi dominasi Pareto standar pada optimasi multi-objektif
fast_non_dominated_sort <- function(F) {
F <- as.matrix(F) # Pastikan F adalah matriks (baris: solusi, kolom: tujuan)
n <- nrow(F) # Jumlah solusi
if (n == 0) return(list(fronts = list(), rank = integer(0)))
Slist <- vector("list", n) # Slist[p]: solusi yang didominasi oleh solusi p
n_dom <- integer(n) # n_dom[p]: berapa solusi yang mendominasi p
ranks <- integer(n) # ranks[p]: ranking Pareto solusi p (front ke berapa)
fronts <- list() # List front Pareto
front1 <- integer(0) # Front 1 (solusi yang tidak didominasi siapapun)
for (p in 1:n) {
Slist[[p]] <- integer(0); n_dom[p] <- 0
for (q in 1:n) if (p != q) { # Pengujian untuk p tidak sama dengan q (tidak membandingkan dengan solusi itu sendiri)
if (dominates(F[p,], F[q,])) Slist[[p]] <- c(Slist[[p]], q) # Jika p mendominasi q, masukkan q ke Slist[[p]]
else if (dominates(F[q,], F[p,])) n_dom[p] <- n_dom[p] + 1 # Jika q mendominasi p, n_dom[p] + 1 (artinya bertambah satu solusi yang mendominasi p)
}
if (n_dom[p] == 0) { ranks[p] <- 1; front1 <- c(front1, p) } # Jika n_dom[p] = 0 (artinya tidak ada solusi yang mendominasi p) tetapkan rank p = 1 dan masukkan ke dalam front 1
}
fronts[[1]] <- front1
# Lanjutkan untuk front berikutnya
i <- 1
while (length(fronts[[i]]) > 0) {
Q <- integer(0)
for (p in fronts[[i]]) for (q in Slist[[p]]) { # Contoh: untuk setiap p pada front 1 dan q pada Slist[p]
n_dom[q] <- n_dom[q] - 1
if (n_dom[q] == 0) { ranks[q] <- i + 1; Q <- c(Q, q) } # Jika n_dom[q] = 0 maka tetapkan rank q = 2 (i + 1)
}
i <- i + 1; fronts[[i]] <- Q # tetapkan i = 2 dan masukkan Q ke dalam front 2 (karena n_dom[q] sudah 0) dan kembali ke iterasi awal
}
fronts[length(fronts)] <- NULL # Hapus front kosong terakhir
list(fronts = fronts, rank = ranks)
# Output: list front Pareto (indeks solusi per front) dan ranking Pareto
}
# Fungsi ini adalah implementasi klasik fast non-dominated sorting pada NSGA-II
crowding_distance <- function(Fvals, idxs) {
if (is.null(idxs) || length(idxs) == 0) return(numeric(0)) # Jika tidak ada solusi, return vektor kosong
if (length(idxs) == 1) return(Inf) # Jika hanya satu solusi, crowding distance = Inf
Fvals <- as.matrix(Fvals)
k <- length(idxs) # Jumlah solusi di front
m <- ncol(Fvals) # Jumlah tujuan/objektif
dist <- rep(0, k) # Inisialisasi crowding distance
for (j in 1:m) { # Untuk setiap tujuan
ord <- order(Fvals[idxs, j]) # Urutkan solusi berdasarkan tujuan j
dist[ord[1]] <- Inf; dist[ord[k]] <- Inf # Solusi ekstrem diberi Inf (selalu dipilih)
v <- Fvals[idxs[ord], j]
denom <- max(v) - min(v)
if (!is.finite(denom) || denom == 0) denom <- 1e-12 # Hindari pembagian nol
for (t in 2:(k-1))
dist[ord[t]] <- dist[ord[t]] + (v[t+1] - v[t-1]) / denom
# Tambahkan crowding distance normalisasi antar tetangga
}
dist
}
.active_pattern <- function(W, eps = 1e-8) {
if (is.null(W) || nrow(W) == 0) return(character(0))
apply(W, 1, function(w) paste(as.integer(w > eps), collapse = ""))
}
# Untuk setiap solusi (baris W), buat string biner yang menandakan aset aktif (bobot > eps)
# Contoh: w = [0.2, 0, 0.3, 0.5] → "1011"
# Ini memungkinkan identifikasi pola portofolio unik di populasi.
.pattern_div_score <- function(cands_W, chosen_patterns) {
if (length(chosen_patterns) == 0) return(rep(1, nrow(cands_W)))
pat_c <- .active_pattern(cands_W) # Pola aktif kandidat solusi
tab <- table(chosen_patterns) # Hitung frekuensi pola yang sudah terpilih
sapply(pat_c, function(p) 1 / (1 + ifelse(is.na(tab[p]), 0, tab[p])))
# Skor tinggi untuk pola yang jarang/unik, skor rendah untuk pola yang sering muncul
# Mendorong diversitas portofolio di ruang keputusan (decision space diversity).
}
# Seleksi NSGA-II + bonus diversitas di ruang keputusan
nsga2_select <- function(P, Fvals, target_size,
cd_threshold = 0.01, use_pattern_tiebreak = TRUE) {
stopifnot(is.matrix(P), is.matrix(Fvals), nrow(P) == nrow(Fvals), target_size > 0)
P <- as.matrix(P); Fvals <- as.matrix(Fvals)
fs <- fast_non_dominated_sort(Fvals); fronts <- fs$fronts
nextP <- matrix(numeric(0), 0, ncol(P)) # Populasi terpilih berikutnya
nextF <- matrix(numeric(0), 0, ncol(Fvals)) # Nilai objektif solusi terpilih
chosen_patterns <- character(0) # Pola aktif yang sudah terpilih
for (f in fronts) {
if (is.null(f) || length(f) == 0) next
remain <- target_size - nrow(nextP) # Sisa slot populasi
if (remain <= 0) break
if (length(f) <= remain) {
valid_indices <- f[f > 0 & f <= nrow(P)]
if (length(valid_indices) == 0) next
# Jika solusi di front ini cukup/tidak melebihi kapasitas, ambil semua
nextP <- rbind(nextP, P[valid_indices, , drop = FALSE])
nextF <- rbind(nextF, Fvals[valid_indices, , drop = FALSE])
pattern_subset <- P[valid_indices, , drop = FALSE]
if (nrow(pattern_subset) > 0) {
chosen_patterns <- c(chosen_patterns, .active_pattern(pattern_subset))
}
} else {
# Step 1: Hitung crowding distance (primary)
dist_obj <- crowding_distance(Fvals, f)
if (is.null(dist_obj) || length(dist_obj) != length(f)) {
warning("crowding_distance error, pakai uniform distance")
dist_obj <- rep(1, length(f))
}
# Step 2: Hitung pattern diversity (SECONDARY)
if (use_pattern_tiebreak && length(f) > 0) {
candW <- P[f, , drop = FALSE]
pat_cands <- .active_pattern(candW)
if (length(pat_cands) > 0) {
div_dec <- sapply(pat_cands, function(p) {
if (length(chosen_patterns) == 0) return(1)
1 / (1 + sum(chosen_patterns == p, na.rm = TRUE))
})
} else {
div_dec <- rep(0, length(f))
}
} else {
div_dec <- rep(0, length(f))
}
# Validasi panjang vektor sebelum data.frame
if (length(f) != length(dist_obj) || length(f) != length(div_dec)) {
min_len <- min(length(f), length(dist_obj), length(div_dec))
f <- f[1:min_len]; dist_obj <- dist_obj[1:min_len]; div_dec <- div_dec[1:min_len]
}
# Step 3: Buat data frame untuk sorting
selection_df <- data.frame(
idx = f,
cd = dist_obj,
pds = div_dec,
stringsAsFactors = FALSE
)
# Step 4: Lexicographic sort
# Primary: CD (descending)
# Secondary: PDS (descending) - hanya untuk CD yang sama
if (nrow(selection_df) == 0) next
selection_df <- selection_df[order(-selection_df$cd, -selection_df$pds), ]
# Step 5: Ambil top-k
take_idx <- head(selection_df$idx, remain)
valid_take_idx <- take_idx[take_idx > 0 & take_idx <= nrow(P)]
if (length(valid_take_idx) == 0) next
nextP <- rbind(nextP, P[valid_take_idx, , drop = FALSE])
nextF <- rbind(nextF, Fvals[valid_take_idx, , drop = FALSE])
pattern_subset <- P[valid_take_idx, , drop = FALSE]
if (nrow(pattern_subset) > 0) {
chosen_patterns <- c(chosen_patterns, .active_pattern(pattern_subset))
}
break
}
}
# Jika populasi terpilih melebihi kapasitas, potong
if (nrow(nextP) > target_size) {
keep <- seq_len(target_size)
nextP <- nextP[keep, , drop = FALSE]
nextF <- nextF[keep, , drop = FALSE]
}
rownames(nextP) <- NULL; rownames(nextF) <- NULL
stopifnot(nrow(nextP) == nrow(nextF), nrow(nextP) <= target_size)
list(P = nextP, Fval = nextF, fronts = fronts)
}
# Seleksi NSGA-II:
# - Ambil solusi dari front Pareto terbaik, jika perlu seleksi dengan crowding distance + pattern diversity
# - Menjaga diversity baik di ruang objektif (risk-return) maupun ruang keputusan (komposisi aset)
## =============================================
## 6) LOOP UTAMA MOCv-ABC
## =============================================
run_mcabc <- function(mu, S, tickers,
alpha = 0.1, beta = 0.4, d_min = 3, d_max = 5,
must_have_name = NULL, must_have_idx = NULL,
sigma = 0.5, limit_factor = 10, max_iter = 700,
seed = 123, debug = FALSE) {
set.seed(seed) # Set seed random agar hasil eksperimen reproducible
d <- length(mu) # Jumlah aset (dimensi portofolio)
stopifnot(ncol(S) == d, nrow(S) == d) # Validasi dimensi matriks kovarians S
# ==========================================
# DEBUGGING INFRASTRUCTURE
# ==========================================
if (debug) {
debug_log <- list(
iteration = integer(),
phase = character(),
metric = character(),
value = numeric(),
detail = character()
)
log_event <- function(iter, phase, metric, value, detail = "") {
debug_log$iteration <<- c(debug_log$iteration, iter)
debug_log$phase <<- c(debug_log$phase, phase)
debug_log$metric <<- c(debug_log$metric, metric)
debug_log$value <<- c(debug_log$value, value)
debug_log$detail <<- c(debug_log$detail, detail)
}
cat("═══════════════════════════════════════════════════════════════\n")
cat(" M-CABC EXECUTION WITH DEBUG MODE \n")
cat("═══════════════════════════════════════════════════════════════\n")
}
# ---------- assert dim ----------
assert_dims <- function(P, Fvals, stage, n_employed, d) {
rP <- NROW(P); cP <- NCOL(P); rF <- NROW(Fvals); cF <- NCOL(Fvals)
if (!(rP == n_employed && cP == d))
stop(sprintf("[%s] dim(P)=%s, expected %d x %d", stage, paste(dim(P), collapse="x"), n_employed, d))
if (!(rF == n_employed && cF == 2))
stop(sprintf("[%s] dim(Fvals)=%s, expected %d x 2", stage, paste(dim(Fvals), collapse="x"), n_employed))
}
# Fungsi utilitas untuk memastikan dimensi populasi dan nilai objektif benar
# Jika tidak sesuai, akan error dan menghentikan eksekusi
# ---------- indeks aset wajib ----------
m_idx <- NULL # Inisialisasi indeks aset wajib
if (!is.null(must_have_name) && nzchar(must_have_name)) {
pos <- match(must_have_name, tickers) # Cari posisi aset wajib berdasarkan nama
if (!is.na(pos)) m_idx <- pos
} else if (!is.null(must_have_idx)) {
if (must_have_idx >= 1 && must_have_idx <= d) m_idx <- must_have_idx # Atau berdasarkan indeks langsung
}
if (debug && !is.null(m_idx)) {
cat(sprintf("✓ Must-have asset: %s (index %d)\n",
tickers[m_idx], m_idx))
}
k_missing_count <- c(`3` = 0, `4` = 0, `5` = 0)
# ---------- ukuran koloni ----------
hive_size <- d * 10L # Ukuran koloni (standar: 10x jumlah aset)
n_employed <- as.integer(floor(hive_size / 2)) # Jumlah employed bees
n_onlooker <- hive_size - n_employed # Jumlah onlooker bees
limit <- as.integer(limit_factor * d) # Batas stagnasi (untuk scout bee)
if (debug) {
cat("───────────────────────────────────────────────────────────────\n")
cat(sprintf("✓ Problem dimension : %d assets\n", d))
cat(sprintf("✓ Hive size : %d\n", hive_size))
cat(sprintf("✓ Employed bees : %d\n", n_employed))
cat(sprintf("✓ Onlooker bees : %d\n", n_onlooker))
cat(sprintf("✓ Scout limit : %d iterations\n", limit))
cat(sprintf("✓ Sigma : %.2f\n", sigma))
cat(sprintf("✓ Cardinality range : [%d, %d]\n", d_min, d_max))
cat(sprintf("✓ Weight bounds : [%.2f, %.2f]\n", alpha, beta))
cat("───────────────────────────────────────────────────────────────\n\n")
}
# ---------- ELITE PER-K & HELPER ----------
elite_by_k <- list(`3` = NULL, `4` = NULL, `5` = NULL)
# Menyimpan portofolio elite (Sharpe tertinggi) untuk setiap jumlah aset aktif (k = 3,4,5)
.metrics_from_W <- function(W) {
ER <- as.numeric(W %*% mu) # Hitung expected return
VAR <- apply(W, 1, function(w) as.numeric(t(w) %*% S %*% w)) # Hitung varians portofolio
SD <- sqrt(VAR) # Hitung standar deviasi (risiko)
Sharpe <- (ER - rf) / pmax(SD, 1e-12) # Hitung sharpe ratio
data.frame(ER, SD, Sharpe) # Return data frame metrik
}
.update_elite <- function(W, elite) {
if (is.null(W) || nrow(W) == 0) return(elite) # Jika tidak ada kandidat, return elite lama
kW <- rowSums(W > 1e-8) # Hitung jumlah aset aktif tiap portofolio
for (kk in c(3L,4L,5L)) { # Untuk setiap k yang di-cover
idx <- which(kW == kk) # Cari portofolio dengan k aset aktif
if (length(idx)) {
M <- .metrics_from_W(W[idx,,drop=FALSE]) # Hitung metriknya
id <- idx[which.max(M$Sharpe)] # Ambil yang sharpe ratio-nya tertinggi
cand <- W[id,,drop=FALSE]
if (is.null(elite[[as.character(kk)]])) {
elite[[as.character(kk)]] <- cand # Jika belum ada, simpan
} else {
old <- .metrics_from_W(elite[[as.character(kk)]])$Sharpe
new <- .metrics_from_W(cand)$Sharpe
if (new > old + 1e-12) elite[[as.character(kk)]] <- cand # Update jika lebih baik dari elite sebelumnya
}
}
}
elite # Return elite yang sudah diupdate
}
# ---------- inisialisasi populasi ----------
init_population <- function(size) {
Pclean <- matrix(NA_real_, nrow = size, ncol = d)
for (i in seq_len(size)) {
# Step 1: Generate random weights dalam bounds [alpha, beta]
w <- alpha + runif(d) * (beta - alpha)
# Step 2: Pre-normalize untuk mendekatkan ke feasible region
w <- w / sum(w) # Enforce sum = 1
# Step 3: Clean untuk enforce semua constraint
Pclean[i, ] <- clean_portfolio(w, d, alpha, beta, d_min, d_max, m_idx)
}
Pclean
}
# ==========================================
# INITIALIZATION
# ==========================================
if (debug) cat("⏳ Initializing population...\n")
P <- init_population(n_employed) # Populasi employed bees awal
Fvals <- t(apply(P, 1, function(w) obj_vec(w, mu, S))) # Hitung nilai objektif awal (risk, -return)
pareto_result <- fast_non_dominated_sort(Fvals)
best_rank_seen <- pareto_result$rank
stagn <- integer(n_employed) # Counter stagnasi untuk setiap employed bee
hist_min_risk <- numeric(max_iter) # History min risiko per iterasi
hist_max_ret <- numeric(max_iter) # History max return per iterasi
# Debug: Initial population stats
if (debug) {
k_init <- table(rowSums(P > 1e-8))
n_fronts <- length(pareto_result$fronts)
f1_size <- length(pareto_result$fronts[[1]])
cat("✓ Population initialized\n")
cat(sprintf(" ├─ Cardinality distribution: %s\n",
paste(sprintf("k=%s:%d", names(k_init), as.integer(k_init)), collapse=", ")))
cat(sprintf(" ├─ Pareto fronts: %d\n", n_fronts))
cat(sprintf(" ├─ Front 1 size: %d (%.1f%%)\n", f1_size, 100*f1_size/n_employed))
cat(sprintf(" ├─ Min risk: %.6f\n", min(Fvals[,1])))
cat(sprintf(" └─ Max return: %.6f\n\n", max(-Fvals[,2])))
log_event(0, "INIT", "n_fronts", n_fronts)
log_event(0, "INIT", "f1_size", f1_size)
log_event(0, "INIT", "min_risk", min(Fvals[,1]))
log_event(0, "INIT", "max_return", max(-Fvals[,2]))
}
# ==========================================
# MAIN LOOP
# ==========================================
if (debug) {
cat("═══════════════════════════════════════════════════════════════\n")
cat(" STARTING ITERATIONS \n")
cat("═══════════════════════════════════════════════════════════════\n\n")
}
for (g in seq_len(max_iter)) { # Loop utama iterasi algoritma
if (debug && (g == 1 || g %% 50 == 0 || g == max_iter)) {
cat(sprintf("\n╔═══ ITERATION %d/%d ═══╗\n", g, max_iter))
}
## =======================================
## EMPLOYED BEES PHASE
## =======================================
phi <- runif(1, -1, 1) # Parameter random eksplorasi
varphi <- runif(1, 0, 1.5) # Parameter random eksploitasi
# Pareto sorting
tmp <- fast_non_dominated_sort(Fvals) # Pareto sorting populasi saat ini
F1 <- tmp$fronts[[1]]
if (length(F1) == 0) F1 <- 1L # Ambil front 1 (Pareto optimal)
BEST <- as.numeric(P[sample(F1, 1), , drop = FALSE])
if (debug && (g == 1 || g %% 50 == 0)) {
best_idx <- sample(F1, 1)
best_metrics <- .metrics_from_W(matrix(P[best_idx,], nrow=1))
cat(sprintf("├─ EMPLOYED PHASE\n"))
cat(sprintf("│ ├─ phi = %.3f, varphi = %.3f\n", phi, varphi))
cat(sprintf("│ ├─ Front 1 size: %d\n", length(F1)))
cat(sprintf("│ └─ BEST solution: idx=%d, SR=%.3f, risk=%.4f, ret=%.4f\n",
best_idx, best_metrics$Sharpe, best_metrics$SD, best_metrics$ER))
}
P_emp <- P # Salin populasi employed bees
for (i in seq_len(n_employed)) { # Untuk setiap employed bee (solusi)
j <- sample(setdiff(seq_len(n_employed), i), 1) # Pilih solusi lain secara acak (bukan diri sendiri)
dm <- sample.int(d, 1) # Pilih satu dimensi (aset) secara acak
Xi <- P[i, ]; Xj <- P[j, ]
Xi_new <- Xi
Xi_new[dm] <- Xi[dm] + phi*(Xi[dm]-Xj[dm]) + varphi*(BEST[dm]-Xi[dm]) # Update posisi bee: eksplorasi + eksploitasi (buat solusi tetangga)
Xi_new <- clean_portfolio(Xi_new, d, alpha, beta, d_min, d_max, m_idx) # Bersihkan agar feasible
P_emp[i, ] <- Xi_new # Simpan solusi baru
}
if (debug && (g == 1 || g %% 50 == 0)) {
log_event(g, "EMPLOYED", "phi", phi)
log_event(g, "EMPLOYED", "varphi", varphi)
log_event(g, "EMPLOYED", "f1_size", length(F1))
}
## =======================================
## ONLOOKER BEES PHASE
## =======================================
fronts <- tmp$fronts
# Step 1: Tentukan jumlah front yang akan digunakan
max_fronts_to_use <- ceiling(log2(n_onlooker))
n_available_fronts <- length(fronts)
n_fronts_to_use <- min(max_fronts_to_use, n_available_fronts)
# Step 2: Hitung alokasi untuk setiap front
# Formula: nof_j = [o / 2^j]
alloc_counts <- integer(n_fronts_to_use)
for (j in 1:n_fronts_to_use) {
alloc_counts[j] <- ceiling(n_onlooker / (2^j))
}
# Step 3: Normalisasi agar total = n_onlooker
total_alloc <- sum(alloc_counts)
if (total_alloc > n_onlooker) {
# Kurangi dari front terakhir
excess <- total_alloc - n_onlooker
alloc_counts[n_fronts_to_use] <- max(1, alloc_counts[n_fronts_to_use] - excess)
} else if (total_alloc < n_onlooker) {
# Tambahkan ke front pertama (yang paling banyak)
deficit <- n_onlooker - total_alloc
alloc_counts[1] <- alloc_counts[1] + deficit
}
# Step 4: Generate allocation vector
alloc <- rep(seq_len(n_fronts_to_use), times = alloc_counts)
if (debug && (g == 1 || g %% 50 == 0)) {
cat(sprintf("├─ ONLOOKER PHASE\n"))
cat(sprintf("│ ├─ Fronts used: %d (max: %d)\n", n_fronts_to_use, max_fronts_to_use))
cat(sprintf("│ └─ Allocation: %s\n",
paste(sprintf("F%d:%d", 1:n_fronts_to_use, alloc_counts), collapse=", ")))
log_event(g, "ONLOOKER", "n_fronts_used", n_fronts_to_use)
log_event(g, "ONLOOKER", "f1_onlookers", alloc_counts[1])
}
# Step 5: Generate onlooker bees
P_onl <- matrix(NA_real_, nrow = n_onlooker, ncol = d)
for (idx in seq_len(n_onlooker)) {
grp <- alloc[idx]
ids <- fronts[[grp]]
if (length(ids) <= 1) {
m <- colMeans(P)
C <- stats::cov(P)
} else {
m <- colMeans(P[ids, , drop = FALSE])
C <- stats::cov(P[ids, , drop = FALSE])
}
eig <- eigen(C, symmetric = TRUE)
Z <- as.numeric(m + sigma * (eig$vectors %*% diag(sqrt(pmax(0, eig$values))) %*% runif(d)))
P_onl[idx, ] <- clean_portfolio(Z, d, alpha, beta, d_min, d_max, m_idx)
}
## =======================================
## SCOUT BEES PHASE
## =======================================
F_emp <- t(apply(P_emp, 1, function(w) obj_vec(w, mu, S))) # Hitung nilai objektif solusi tetangga employed bees
# Step 1: Gabungkan populasi lama dan baru untuk comparison
P_for_rank <- rbind(P, P_emp)
F_for_rank <- rbind(Fvals, F_emp)
# Step2 2: Hitung Pareto rank
rank_result <- fast_non_dominated_sort(F_for_rank)
all_ranks <- rank_result$rank
# Step 3: Pisahkan rank
current_ranks <- all_ranks[(n_employed + 1):(2 * n_employed)]
old_ranks <- all_ranks[1:n_employed]
# Step 4: Check improvement (rank membaik atau sama)
improved <- current_ranks <= best_rank_seen
stagn[improved] <- 0L
stagn[!improved] <- stagn[!improved] + 1L
# Step 5: Update best rank seen
best_rank_seen <- pmin(best_rank_seen, current_ranks)
# Step 6: Scout bee (re-initialize stagnant)
idx_stag <- which(stagn >= limit)
if (debug && (g == 1 || g %% 50 == 0)) {
n_improved <- sum(improved)
max_stag <- max(stagn)
cat(sprintf("├─ SCOUT PHASE\n"))
cat(sprintf("│ ├─ Improved solutions: %d/%d (%.1f%%)\n",
n_improved, n_employed, 100*n_improved/n_employed))
cat(sprintf("│ ├─ Max stagnation: %d (limit: %d)\n", max_stag, limit))
cat(sprintf("│ └─ Re-initialized: %d bees\n", length(idx_stag)))
log_event(g, "SCOUT", "n_improved", n_improved)
log_event(g, "SCOUT", "n_reinit", length(idx_stag))
log_event(g, "SCOUT", "max_stagnation", max_stag)
}
if (length(idx_stag)) {
for (ii in idx_stag) {
# Gunakan fungsi inisialisasi
w <- alpha + runif(d) * (beta - alpha)
w <- w / sum(w)
P_emp[ii, ] <- clean_portfolio(w, d, alpha, beta, d_min, d_max, m_idx)
}
stagn[idx_stag] <- 0L
# Re-evaluate
F_emp <- t(apply(P_emp, 1, function(w) obj_vec(w, mu, S)))
# Update ranks setelah re-inisialisasi
P_for_rank <- rbind(P, P_emp)
F_for_rank <- rbind(Fvals, F_emp)
rank_result <- fast_non_dominated_sort(F_for_rank)
current_ranks <- rank_result$rank[(n_employed + 1):(2 * n_employed)]
best_rank_seen <- pmin(best_rank_seen, current_ranks)
}
## =======================================
## TARGETED SCOUT (k-coverage)
## =======================================
ensure_k <- c(3L, 4L, 5L) # Daftar jumlah aset aktif yang wajib ada di populasi
k_now <- rowSums(P_emp > 1e-8) # Hitung jumlah aset aktif tiap solusi
missing_k <- setdiff(ensure_k, unique(k_now)) # Cari k yang belum ada di populasi
if (length(missing_k) > 0) {
# Update counter
for (kk in missing_k) {
k_missing_count[as.character(kk)] <- k_missing_count[as.character(kk)] + 1
}
# Adaptive Rate
n_missing <- length(missing_k)
base_rate <- 0.03 * n_missing
max_missing_count <- max(k_missing_count[as.character(missing_k)])
adaptive_multiplier <- min(3.0, 1.0 + (max_missing_count / 10.0))
injection_rate <- min(0.15, base_rate * adaptive_multiplier)
n_inject <- max(2L, floor(injection_rate * n_employed))
if (debug && (g == 1 || g %% 50 == 0)) {
cat(sprintf("├─ TARGETED SCOUT (ADAPTIVE)\n"))
cat(sprintf("│ ├─ Missing k=%s (count: %s)\n",
paste(missing_k, collapse=","),
paste(k_missing_count[as.character(missing_k)], collapse=",")))
cat(sprintf("│ ├─ Multiplier: %.2fx → Rate: %.1f%%\n",
adaptive_multiplier, injection_rate*100))
cat(sprintf("│ └─ Injecting: %d solutions\n", n_inject))
log_event(g, "TARGETED", "n_missing_k", n_missing)
log_event(g, "TARGETED", "n_injected", n_inject)
}
for (kk in missing_k) {
idx_reseed <- sample(seq_len(n_employed), n_inject)
for (ii in idx_reseed) {
w <- alpha + runif(d) * (beta - alpha)
w <- w / sum(w)
P_emp[ii, ] <- clean_portfolio(w, d, alpha, beta, kk, kk, m_idx)
}
}
F_emp <- t(apply(P_emp, 1, function(w) obj_vec(w, mu, S)))
}
elite_by_k <- .update_elite(P_emp, elite_by_k)
## =======================================
## NSGA-II SELECTION
## =======================================
F_onl <- t(apply(P_onl, 1, function(w) obj_vec(w, mu, S))) # Hitung nilai objektif solusi onlooker bees
elite_stack <- do.call(rbind, elite_by_k[!sapply(elite_by_k, is.null)]) # Gabungkan semua solusi
if (!is.null(elite_stack) && nrow(elite_stack) > 0) { # Jika ada solusi elite
F_elite <- t(apply(elite_stack, 1, function(w) obj_vec(w, mu, S))) # Hitung nilai objektifnya
P_comb <- rbind(P, P_emp, P_onl, elite_stack) # Gabungkan semua solusi (lama, employed, onlooker, elite)
F_comb <- rbind(Fvals, F_emp, F_onl, F_elite)
n_elite <- nrow(elite_stack)
} else {
P_comb <- rbind(P, P_emp, P_onl) # Gabungkan tanpa elite
F_comb <- rbind(Fvals, F_emp, F_onl)
n_elite <- 0
}
if (debug && (g == 1 || g %% 50 == 0)) {
cat(sprintf("├─ SELECTION\n"))
cat(sprintf("│ ├─ Combined pool: %d solutions\n", nrow(P_comb)))
cat(sprintf("│ │ ├─ Old generation: %d\n", n_employed))
cat(sprintf("│ │ ├─ Employed: %d\n", n_employed))
cat(sprintf("│ │ ├─ Onlooker: %d\n", n_onlooker))
cat(sprintf("│ │ └─ Elite: %d\n", n_elite))
log_event(g, "SELECTION", "pool_size", nrow(P_comb))
log_event(g, "SELECTION", "n_elite", n_elite)
}
sel <- nsga2_select(P_comb, F_comb, n_employed) # Seleksi generasi berikutnya dengan NSGA-II (Pareto + crowding + pattern_diversity tie breaker)
P_selected <- as.matrix(sel$P)
Fvals_selected <- as.matrix(sel$Fval)
# ============================================
# FORCE CARDINALITY COVERAGE
# ============================================
k_selected <- rowSums(P_selected > 1e-8)
missing_k_final <- setdiff(ensure_k, unique(k_selected))
if (length(missing_k_final) > 0 && !is.null(elite_stack)) {
# Identifikasi elite untuk missing k
elite_indices_to_force <- which(sapply(names(elite_by_k), function(kk) {
as.integer(kk) %in% missing_k_final && !is.null(elite_by_k[[kk]])
}))
if (length(elite_indices_to_force) > 0) {
# Ambil elite yang missing
elite_to_force <- do.call(rbind, elite_by_k[elite_indices_to_force])
F_elite_to_force <- t(apply(elite_to_force, 1, function(w) obj_vec(w, mu, S)))
# Hitung rank elite
P_test <- rbind(P_selected, elite_to_force)
F_test <- rbind(Fvals_selected, F_elite_to_force)
rank_test <- fast_non_dominated_sort(F_test)$rank
elite_ranks <- rank_test[(n_employed + 1):nrow(P_test)]
# Replace worst solutions dengan elite
# Strategy: Replace solutions dengan rank tertinggi (worst)
current_ranks <- rank_test[1:n_employed]
worst_rank <- max(current_ranks)
worst_indices <- which(current_ranks == worst_rank)
n_to_replace <- min(length(elite_indices_to_force), length(worst_indices))
if (n_to_replace > 0) {
# Replace dari belakang (worst crowding distance di rank terburuk)
replace_idx <- tail(worst_indices, n_to_replace)
P_selected[replace_idx, ] <- elite_to_force[1:n_to_replace, ]
Fvals_selected[replace_idx, ] <- F_elite_to_force[1:n_to_replace, ]
}
}
}
P <- P_selected
Fvals <- Fvals_selected
# ==========================================
# ITERATION SUMMARY
# ==========================================
hist_min_risk[g] <- min(Fvals[,1]) # Simpan min risiko iterasi ini
hist_max_ret[g] <- max(-Fvals[,2]) # Simpan max return iterasi ini
if (debug && (g == 1 || g %% 50 == 0 || g == max_iter)) {
k_final <- table(rowSums(P > 1e-8))
final_fronts <- fast_non_dominated_sort(Fvals)
# Elite stats
elite_sharpe <- sapply(elite_by_k, function(e) {
if (is.null(e)) return(NA)
.metrics_from_W(matrix(e, nrow=1))$Sharpe
})
cat(sprintf("└─ ITERATION SUMMARY\n"))
cat(sprintf(" ├─ Cardinality: %s\n",
paste(sprintf("k=%s:%d", names(k_final), as.integer(k_final)), collapse=", ")))
cat(sprintf(" ├─ Fronts: %d (F1: %d solutions)\n",
length(final_fronts$fronts), length(final_fronts$fronts[[1]])))
cat(sprintf(" ├─ Min risk: %.6f\n", hist_min_risk[g]))
cat(sprintf(" ├─ Max return: %.6f\n", hist_max_ret[g]))
cat(sprintf(" └─ Elite Sharpe: k3=%.3f, k4=%.3f, k5=%.3f\n",
ifelse(is.na(elite_sharpe[1]), 0, elite_sharpe[1]),
ifelse(is.na(elite_sharpe[2]), 0, elite_sharpe[2]),
ifelse(is.na(elite_sharpe[3]), 0, elite_sharpe[3])))
if (g == 1 || g %% 50 == 0) cat("╚═══════════════════════════════════════════════════════════════╝\n")
log_event(g, "SUMMARY", "min_risk", hist_min_risk[g])
log_event(g, "SUMMARY", "max_return", hist_max_ret[g])
log_event(g, "SUMMARY", "n_fronts", length(final_fronts$fronts))
log_event(g, "SUMMARY", "f1_size", length(final_fronts$fronts[[1]]))
}
}
# ==========================================
# FINAL REPORT
# ==========================================
if (debug) {
cat("\n═══════════════════════════════════════════════════════════════\n")
cat(" EXECUTION COMPLETED \n")
cat("═══════════════════════════════════════════════════════════════\n")
final_fronts <- fast_non_dominated_sort(Fvals)
k_final <- table(rowSums(P > 1e-8))
cat(sprintf("✓ Total iterations: %d\n", max_iter))
cat(sprintf("✓ Final population:\n"))
cat(sprintf(" ├─ Pareto fronts: %d\n", length(final_fronts$fronts)))
cat(sprintf(" ├─ Front 1 size: %d (%.1f%%)\n",
length(final_fronts$fronts[[1]]),
100*length(final_fronts$fronts[[1]])/n_employed))
cat(sprintf(" ├─ Cardinality: %s\n",
paste(sprintf("k=%s:%d", names(k_final), as.integer(k_final)), collapse=", ")))
cat(sprintf(" ├─ Best risk: %.6f\n", min(Fvals[,1])))
cat(sprintf(" └─ Best return: %.6f\n", max(-Fvals[,2])))
# Elite summary
cat("\n✓ Elite portfolios (best Sharpe per cardinality):\n")
for (kk in c(3,4,5)) {
e <- elite_by_k[[as.character(kk)]]
if (!is.null(e)) {
m <- .metrics_from_W(matrix(e, nrow=1))
active_assets <- which(e > 1e-8)
cat(sprintf(" ├─ k=%d: Sharpe=%.4f, Risk=%.4f, Return=%.4f\n",
kk, m$Sharpe, m$SD, m$ER))
cat(sprintf(" │ Assets: %s\n",
paste(tickers[active_assets], collapse=", ")))
} else {
cat(sprintf(" ├─ k=%d: Not found\n", kk))
}
}
cat("═══════════════════════════════════════════════════════════════\n\n")
# Save debug log
debug_df <- as.data.frame(debug_log)
attr(P, "debug_log") <- debug_df
}
# Return
result <- list(
P = P,
Fvals = Fvals,
tickers = tickers,
elite_by_k = elite_by_k,
history = data.frame(
iter = seq_len(max_iter),
min_risk = hist_min_risk,
max_return = hist_max_ret
)
)
if (debug) {
result$debug_log <- as.data.frame(debug_log)
}
return(result)
}
## =============================================
## 7) JALANKAN OPTIMASI & VISUALISASI FRONTIER
## =============================================
mcabc_cov <- run_mcabc(
mu = mu, S = S, tickers = tickers,
alpha = alpha_w, beta = beta_w,
d_min = d_min, d_max = d_max,
must_have_name = must_have_name,
sigma = sigma_abc,
limit_factor = limit_factor,
max_iter = max_iter,
debug = TRUE
)
## ═══════════════════════════════════════════════════════════════
## M-CABC EXECUTION WITH DEBUG MODE
## ═══════════════════════════════════════════════════════════════
## ✓ Must-have asset: TLKM (index 15)
## ───────────────────────────────────────────────────────────────
## ✓ Problem dimension : 17 assets
## ✓ Hive size : 170
## ✓ Employed bees : 85
## ✓ Onlooker bees : 85
## ✓ Scout limit : 170 iterations
## ✓ Sigma : 0.50
## ✓ Cardinality range : [3, 5]
## ✓ Weight bounds : [0.10, 0.40]
## ───────────────────────────────────────────────────────────────
##
## ⏳ Initializing population...
## ✓ Population initialized
## ├─ Cardinality distribution: k=5:85
## ├─ Pareto fronts: 17
## ├─ Front 1 size: 7 (8.2%)
## ├─ Min risk: 0.175542
## └─ Max return: 0.203194
##
## ═══════════════════════════════════════════════════════════════
## STARTING ITERATIONS
## ═══════════════════════════════════════════════════════════════
##
##
## ╔═══ ITERATION 1/200 ═══╗
## ├─ EMPLOYED PHASE
## │ ├─ phi = 0.150, varphi = 0.037
## │ ├─ Front 1 size: 7
## │ └─ BEST solution: idx=34, SR=-0.001, risk=0.1790, ret=0.0473
## ├─ ONLOOKER PHASE
## │ ├─ Fronts used: 7 (max: 7)
## │ └─ Allocation: F1:43, F2:22, F3:11, F4:6, F5:3, F6:2, F7:1
## ├─ SCOUT PHASE
## │ ├─ Improved solutions: 20/85 (23.5%)
## │ ├─ Max stagnation: 1 (limit: 170)
## │ └─ Re-initialized: 0 bees
## ├─ TARGETED SCOUT (ADAPTIVE)
## │ ├─ Missing k=3,4 (count: 1,1)
## │ ├─ Multiplier: 1.10x → Rate: 6.6%
## │ └─ Injecting: 5 solutions
## ├─ SELECTION
## │ ├─ Combined pool: 258 solutions
## │ │ ├─ Old generation: 85
## │ │ ├─ Employed: 85
## │ │ ├─ Onlooker: 85
## │ │ └─ Elite: 3
## └─ ITERATION SUMMARY
## ├─ Cardinality: k=3:2, k=4:1, k=5:82
## ├─ Fronts: 7 (F1: 14 solutions)
## ├─ Min risk: 0.173130
## ├─ Max return: 0.203194
## └─ Elite Sharpe: k3=0.421, k4=-0.149, k5=0.711
## ╚═══════════════════════════════════════════════════════════════╝
##
## ╔═══ ITERATION 50/200 ═══╗
## ├─ EMPLOYED PHASE
## │ ├─ phi = 0.896, varphi = 1.474
## │ ├─ Front 1 size: 84
## │ └─ BEST solution: idx=69, SR=0.800, risk=0.2560, ret=0.2523
## ├─ ONLOOKER PHASE
## │ ├─ Fronts used: 2 (max: 7)
## │ └─ Allocation: F1:63, F2:22
## ├─ SCOUT PHASE
## │ ├─ Improved solutions: 66/85 (77.6%)
## │ ├─ Max stagnation: 7 (limit: 170)
## │ └─ Re-initialized: 0 bees
## ├─ TARGETED SCOUT (ADAPTIVE)
## │ ├─ Missing k=3 (count: 14)
## │ ├─ Multiplier: 2.40x → Rate: 7.2%
## │ └─ Injecting: 6 solutions
## ├─ SELECTION
## │ ├─ Combined pool: 258 solutions
## │ │ ├─ Old generation: 85
## │ │ ├─ Employed: 85
## │ │ ├─ Onlooker: 85
## │ │ └─ Elite: 3
## └─ ITERATION SUMMARY
## ├─ Cardinality: k=3:1, k=4:2, k=5:82
## ├─ Fronts: 2 (F1: 84 solutions)
## ├─ Min risk: 0.169115
## ├─ Max return: 0.283916
## └─ Elite Sharpe: k3=0.718, k4=0.802, k5=0.816
## ╚═══════════════════════════════════════════════════════════════╝
##
## ╔═══ ITERATION 100/200 ═══╗
## ├─ EMPLOYED PHASE
## │ ├─ phi = 0.265, varphi = 0.288
## │ ├─ Front 1 size: 84
## │ └─ BEST solution: idx=2, SR=0.803, risk=0.3104, ret=0.2967
## ├─ ONLOOKER PHASE
## │ ├─ Fronts used: 2 (max: 7)
## │ └─ Allocation: F1:63, F2:22
## ├─ SCOUT PHASE
## │ ├─ Improved solutions: 74/85 (87.1%)
## │ ├─ Max stagnation: 57 (limit: 170)
## │ └─ Re-initialized: 0 bees
## ├─ TARGETED SCOUT (ADAPTIVE)
## │ ├─ Missing k=3 (count: 24)
## │ ├─ Multiplier: 3.00x → Rate: 9.0%
## │ └─ Injecting: 7 solutions
## ├─ SELECTION
## │ ├─ Combined pool: 258 solutions
## │ │ ├─ Old generation: 85
## │ │ ├─ Employed: 85
## │ │ ├─ Onlooker: 85
## │ │ └─ Elite: 3
## └─ ITERATION SUMMARY
## ├─ Cardinality: k=3:1, k=4:7, k=5:77
## ├─ Fronts: 2 (F1: 84 solutions)
## ├─ Min risk: 0.168992
## ├─ Max return: 0.299741
## └─ Elite Sharpe: k3=0.753, k4=0.815, k5=0.817
## ╚═══════════════════════════════════════════════════════════════╝
##
## ╔═══ ITERATION 150/200 ═══╗
## ├─ EMPLOYED PHASE
## │ ├─ phi = 0.729, varphi = 0.357
## │ ├─ Front 1 size: 84
## │ └─ BEST solution: idx=18, SR=0.754, risk=0.2240, ret=0.2164
## ├─ ONLOOKER PHASE
## │ ├─ Fronts used: 2 (max: 7)
## │ └─ Allocation: F1:63, F2:22
## ├─ SCOUT PHASE
## │ ├─ Improved solutions: 65/85 (76.5%)
## │ ├─ Max stagnation: 107 (limit: 170)
## │ └─ Re-initialized: 0 bees
## ├─ SELECTION
## │ ├─ Combined pool: 258 solutions
## │ │ ├─ Old generation: 85
## │ │ ├─ Employed: 85
## │ │ ├─ Onlooker: 85
## │ │ └─ Elite: 3
## └─ ITERATION SUMMARY
## ├─ Cardinality: k=3:1, k=4:12, k=5:72
## ├─ Fronts: 2 (F1: 84 solutions)
## ├─ Min risk: 0.168957
## ├─ Max return: 0.301112
## └─ Elite Sharpe: k3=0.753, k4=0.815, k5=0.817
## ╚═══════════════════════════════════════════════════════════════╝
##
## ╔═══ ITERATION 200/200 ═══╗
## ├─ EMPLOYED PHASE
## │ ├─ phi = 0.909, varphi = 0.411
## │ ├─ Front 1 size: 84
## │ └─ BEST solution: idx=67, SR=0.810, risk=0.3147, ret=0.3025
## ├─ ONLOOKER PHASE
## │ ├─ Fronts used: 2 (max: 7)
## │ └─ Allocation: F1:63, F2:22
## ├─ SCOUT PHASE
## │ ├─ Improved solutions: 66/85 (77.6%)
## │ ├─ Max stagnation: 157 (limit: 170)
## │ └─ Re-initialized: 0 bees
## ├─ SELECTION
## │ ├─ Combined pool: 258 solutions
## │ │ ├─ Old generation: 85
## │ │ ├─ Employed: 85
## │ │ ├─ Onlooker: 85
## │ │ └─ Elite: 3
## └─ ITERATION SUMMARY
## ├─ Cardinality: k=3:1, k=4:13, k=5:71
## ├─ Fronts: 2 (F1: 84 solutions)
## ├─ Min risk: 0.168943
## ├─ Max return: 0.303047
## └─ Elite Sharpe: k3=0.753, k4=0.815, k5=0.817
## ╚═══════════════════════════════════════════════════════════════╝
##
## ═══════════════════════════════════════════════════════════════
## EXECUTION COMPLETED
## ═══════════════════════════════════════════════════════════════
## ✓ Total iterations: 200
## ✓ Final population:
## ├─ Pareto fronts: 2
## ├─ Front 1 size: 84 (98.8%)
## ├─ Cardinality: k=3:1, k=4:13, k=5:71
## ├─ Best risk: 0.168943
## └─ Best return: 0.303047
##
## ✓ Elite portfolios (best Sharpe per cardinality):
## ├─ k=3: Sharpe=0.7532, Risk=0.3130, Return=0.2832
## │ Assets: MEDC, SSIA, TLKM
## ├─ k=4: Sharpe=0.8151, Risk=0.2948, Return=0.2878
## │ Assets: MEDC, SSIA, TLKM, UNTR
## ├─ k=5: Sharpe=0.8170, Risk=0.2773, Return=0.2741
## │ Assets: MEDC, PGAS, SSIA, TLKM, UNTR
## ═══════════════════════════════════════════════════════════════
P <- mcabc_cov$P; Fvals <- mcabc_cov$Fvals
# Plot frontier (semua solusi, highlight Pareto)
fs <- fast_non_dominated_sort(Fvals)
F1 <- fs$fronts[[1]]
cat("Jumlah solusi di Front-1:", length(F1), "\n")
## Jumlah solusi di Front-1: 84
dfF <- as.data.frame(Fvals); colnames(dfF) <- c("risk_sd","neg_ret")
dfF$return <- -dfF$neg_ret; dfF$pareto <- FALSE; dfF$pareto[F1] <- TRUE
p_front <- ggplot(dfF, aes(risk_sd, return, color = pareto)) +
geom_point(alpha = .8) +
scale_color_manual(values = c("#B0B0B0","#2E86DE")) +
labs(title = "Risk-Return Frontier (SD vs Return)",
x = "Risk (SD, annualized)", y = "Return (annualized)", color = NULL)
plotly::ggplotly(p_front)
# Frekuensi keikutsertaan aset di F1
P_F1 <- P[F1, , drop = FALSE]
freq <- colSums(P_F1 > 1e-8); pct <- 100 * freq / nrow(P_F1)
freq_df <- tibble::tibble(ticker = tickers, freq = freq, pct = pct) |> dplyr::arrange(desc(pct))
print(freq_df)
## # A tibble: 17 × 3
## ticker freq pct
## <chr> <dbl> <dbl>
## 1 TLKM 84 100
## 2 UNTR 83 98.8
## 3 SSIA 67 79.8
## 4 MEDC 52 61.9
## 5 PGAS 46 54.8
## 6 BBNI 26 31.0
## 7 BMRI 26 31.0
## 8 MTEL 17 20.2
## 9 JSMR 4 4.76
## 10 AKRA 2 2.38
## 11 BBRI 0 0
## 12 EXCL 0 0
## 13 INTP 0 0
## 14 ISAT 0 0
## 15 PTPP 0 0
## 16 SMGR 0 0
## 17 TOWR 0 0
top5_freq <- head(freq_df, 5)
p_freq <- ggplot(top5_freq, aes(reorder(ticker, pct), pct)) +
geom_col(fill = "#2E86DE") + coord_flip() +
labs(title = "Top-5 Saham Paling Sering di Front-1", x = NULL, y = "Frekuensi (%)")
plotly::ggplotly(p_freq)
# Ringkasan jumlah solusi per-k (populasi & Pareto-1)
k_each <- rowSums(P > 1e-8)
count_pop <- as.data.frame(table(k = k_each)); names(count_pop) <- c("k","count_in_population"); count_pop$k <- as.integer(as.character(count_pop$k))
count_f1 <- as.data.frame(table(k = k_each[F1])); names(count_f1) <- c("k","count_in_F1"); count_f1$k <- as.integer(as.character(count_f1$k))
counts_summary <- dplyr::full_join(count_pop, count_f1, by = "k") |>
dplyr::arrange(k) |>
dplyr::mutate(dplyr::across(dplyr::starts_with("count"), ~replace(., is.na(.), 0L)))
print(counts_summary)
## k count_in_population count_in_F1
## 1 3 1 0
## 2 4 13 13
## 3 5 71 71
debug_df <- mcabc_cov$debug_log
ggplot(debug_df %>% filter(metric == "min_risk"),
aes(x = iteration, y = value)) +
geom_line() +
labs(title = "Risk Convergence",
x = "Iteration", y = "Minimum Risk") +
theme_minimal()
ggplot(debug_df %>% filter(metric == "max_return"),
aes(x = iteration, y = value)) +
geom_line() +
labs(title = "Return Convergence",
x = "Iteration", y = "Maximum Return") +
theme_minimal()
mean_risk_history <- mcabc_cov$history
mean_risk_history
## iter min_risk max_return
## 1 1 0.1731297 0.2031939
## 2 2 0.1714163 0.2109224
## 3 3 0.1714163 0.2109224
## 4 4 0.1701895 0.2109224
## 5 5 0.1701895 0.2118447
## 6 6 0.1701895 0.2118447
## 7 7 0.1701895 0.2136776
## 8 8 0.1701895 0.2136776
## 9 9 0.1701895 0.2208829
## 10 10 0.1701895 0.2214201
## 11 11 0.1701895 0.2214201
## 12 12 0.1701895 0.2214201
## 13 13 0.1698612 0.2214201
## 14 14 0.1698612 0.2236294
## 15 15 0.1698612 0.2236294
## 16 16 0.1698612 0.2236294
## 17 17 0.1698612 0.2239790
## 18 18 0.1697364 0.2408740
## 19 19 0.1693840 0.2422052
## 20 20 0.1693840 0.2424305
## 21 21 0.1693840 0.2424305
## 22 22 0.1693840 0.2424305
## 23 23 0.1693840 0.2431988
## 24 24 0.1693840 0.2431988
## 25 25 0.1693840 0.2464948
## 26 26 0.1691642 0.2643594
## 27 27 0.1691642 0.2643594
## 28 28 0.1691642 0.2686975
## 29 29 0.1691642 0.2686975
## 30 30 0.1691642 0.2690507
## 31 31 0.1691642 0.2690507
## 32 32 0.1691338 0.2696328
## 33 33 0.1691338 0.2696328
## 34 34 0.1691338 0.2706578
## 35 35 0.1691338 0.2783484
## 36 36 0.1691338 0.2783484
## 37 37 0.1691338 0.2783484
## 38 38 0.1691338 0.2783484
## 39 39 0.1691338 0.2783484
## 40 40 0.1691151 0.2821111
## 41 41 0.1691151 0.2821111
## 42 42 0.1691151 0.2821111
## 43 43 0.1691151 0.2821111
## 44 44 0.1691151 0.2821111
## 45 45 0.1691151 0.2839163
## 46 46 0.1691151 0.2839163
## 47 47 0.1691151 0.2839163
## 48 48 0.1691151 0.2839163
## 49 49 0.1691151 0.2839163
## 50 50 0.1691151 0.2839163
## 51 51 0.1691151 0.2839163
## 52 52 0.1691151 0.2839163
## 53 53 0.1691151 0.2911587
## 54 54 0.1691151 0.2911587
## 55 55 0.1691151 0.2911587
## 56 56 0.1691151 0.2911587
## 57 57 0.1691151 0.2911587
## 58 58 0.1691151 0.2911587
## 59 59 0.1691151 0.2911587
## 60 60 0.1691151 0.2911587
## 61 61 0.1691148 0.2911587
## 62 62 0.1691148 0.2911587
## 63 63 0.1691148 0.2911587
## 64 64 0.1691148 0.2911587
## 65 65 0.1691148 0.2911587
## 66 66 0.1691148 0.2911587
## 67 67 0.1691148 0.2911587
## 68 68 0.1691148 0.2911587
## 69 69 0.1691148 0.2911587
## 70 70 0.1691148 0.2911587
## 71 71 0.1691148 0.2919384
## 72 72 0.1691148 0.2919384
## 73 73 0.1691148 0.2919384
## 74 74 0.1691148 0.2943328
## 75 75 0.1691148 0.2943328
## 76 76 0.1691148 0.2943328
## 77 77 0.1691148 0.2943328
## 78 78 0.1691148 0.2943328
## 79 79 0.1691148 0.2943328
## 80 80 0.1691148 0.2943328
## 81 81 0.1691148 0.2943328
## 82 82 0.1690296 0.2943328
## 83 83 0.1690296 0.2943328
## 84 84 0.1689966 0.2943328
## 85 85 0.1689966 0.2943328
## 86 86 0.1689966 0.2943328
## 87 87 0.1689966 0.2943328
## 88 88 0.1689966 0.2943328
## 89 89 0.1689966 0.2943328
## 90 90 0.1689966 0.2943328
## 91 91 0.1689966 0.2943328
## 92 92 0.1689966 0.2943328
## 93 93 0.1689966 0.2959032
## 94 94 0.1689966 0.2959032
## 95 95 0.1689966 0.2959032
## 96 96 0.1689966 0.2959032
## 97 97 0.1689966 0.2959032
## 98 98 0.1689966 0.2966774
## 99 99 0.1689966 0.2966774
## 100 100 0.1689922 0.2997415
## 101 101 0.1689922 0.2997415
## 102 102 0.1689922 0.2997415
## 103 103 0.1689922 0.2997415
## 104 104 0.1689922 0.2997415
## 105 105 0.1689922 0.2997415
## 106 106 0.1689922 0.2997415
## 107 107 0.1689922 0.2997415
## 108 108 0.1689922 0.2997415
## 109 109 0.1689922 0.2997415
## 110 110 0.1689922 0.2997415
## 111 111 0.1689922 0.2997415
## 112 112 0.1689922 0.2997415
## 113 113 0.1689922 0.2997415
## 114 114 0.1689922 0.2997415
## 115 115 0.1689922 0.2997415
## 116 116 0.1689922 0.2997415
## 117 117 0.1689922 0.2997415
## 118 118 0.1689922 0.2997415
## 119 119 0.1689922 0.2997415
## 120 120 0.1689922 0.2997415
## 121 121 0.1689922 0.2997415
## 122 122 0.1689922 0.2997415
## 123 123 0.1689922 0.2997415
## 124 124 0.1689922 0.2997415
## 125 125 0.1689922 0.2997415
## 126 126 0.1689922 0.2997415
## 127 127 0.1689922 0.2997415
## 128 128 0.1689922 0.3006859
## 129 129 0.1689653 0.3006859
## 130 130 0.1689653 0.3006859
## 131 131 0.1689653 0.3006859
## 132 132 0.1689653 0.3006859
## 133 133 0.1689653 0.3006859
## 134 134 0.1689653 0.3006859
## 135 135 0.1689653 0.3006859
## 136 136 0.1689653 0.3006859
## 137 137 0.1689653 0.3007455
## 138 138 0.1689653 0.3007455
## 139 139 0.1689566 0.3007455
## 140 140 0.1689566 0.3007455
## 141 141 0.1689566 0.3007455
## 142 142 0.1689566 0.3010206
## 143 143 0.1689566 0.3010206
## 144 144 0.1689566 0.3010206
## 145 145 0.1689566 0.3010206
## 146 146 0.1689566 0.3010206
## 147 147 0.1689566 0.3010206
## 148 148 0.1689566 0.3010206
## 149 149 0.1689566 0.3010206
## 150 150 0.1689566 0.3011122
## 151 151 0.1689566 0.3011122
## 152 152 0.1689566 0.3011122
## 153 153 0.1689566 0.3011122
## 154 154 0.1689566 0.3011122
## 155 155 0.1689566 0.3011122
## 156 156 0.1689566 0.3022828
## 157 157 0.1689566 0.3022828
## 158 158 0.1689566 0.3022828
## 159 159 0.1689566 0.3022828
## 160 160 0.1689566 0.3022828
## 161 161 0.1689566 0.3022828
## 162 162 0.1689566 0.3022828
## 163 163 0.1689566 0.3022828
## 164 164 0.1689566 0.3022828
## 165 165 0.1689566 0.3022828
## 166 166 0.1689566 0.3022828
## 167 167 0.1689566 0.3022828
## 168 168 0.1689566 0.3022828
## 169 169 0.1689566 0.3022828
## 170 170 0.1689566 0.3022828
## 171 171 0.1689430 0.3022828
## 172 172 0.1689430 0.3022828
## 173 173 0.1689430 0.3022828
## 174 174 0.1689430 0.3022828
## 175 175 0.1689430 0.3024687
## 176 176 0.1689430 0.3024687
## 177 177 0.1689430 0.3024687
## 178 178 0.1689430 0.3024687
## 179 179 0.1689430 0.3024687
## 180 180 0.1689430 0.3024687
## 181 181 0.1689430 0.3024687
## 182 182 0.1689430 0.3024687
## 183 183 0.1689430 0.3024687
## 184 184 0.1689430 0.3024687
## 185 185 0.1689430 0.3024687
## 186 186 0.1689430 0.3024687
## 187 187 0.1689430 0.3024687
## 188 188 0.1689430 0.3024687
## 189 189 0.1689430 0.3024687
## 190 190 0.1689430 0.3024687
## 191 191 0.1689430 0.3024687
## 192 192 0.1689430 0.3024687
## 193 193 0.1689430 0.3024687
## 194 194 0.1689430 0.3030473
## 195 195 0.1689430 0.3030473
## 196 196 0.1689430 0.3030473
## 197 197 0.1689430 0.3030473
## 198 198 0.1689430 0.3030473
## 199 199 0.1689430 0.3030473
## 200 200 0.1689430 0.3030473
## =============================================
## 8) EMPAT PORTOFOLIO TERPILIH:
## (1) Cluster manual + (2-4) MOCv-ABC k=3/4/5 (Sharpe tertinggi)
## =============================================
# Helper metrik & pemilihan berdasarkan Sharpe
portfolio_metrics <- function(W, mu, S, rf) {
ER <- as.numeric(W %*% mu)
VAR <- apply(W, 1, function(w) as.numeric(t(w) %*% S %*% w))
SD <- sqrt(VAR)
Sharpe <- (ER - rf) / pmax(SD, 1e-12)
data.frame(ER = ER, SD = SD, Sharpe = Sharpe, row.names = NULL)
}
pick_best_by_k <- function(P, mu, S, rf, idx_rows) {
if (length(idx_rows) == 0) return(NULL)
Wk <- P[idx_rows, , drop = FALSE]; M <- portfolio_metrics(Wk, mu, S, rf)
i_local <- which.max(M$Sharpe)
list(w = Wk[i_local, ], ER = M$ER[i_local], SD = M$SD[i_local], Sharpe = M$Sharpe[i_local],
idx_global = idx_rows[i_local])
}
# (A) Portofolio #1 dari hasil cluster/dendrogram (input manual saham)
cluster_manual <- c("PGAS", "EXCL", "MTEL", "TLKM", "BBRI") # EDIT sesuai hasil dendrogrammu
stopifnot(all(cluster_manual %in% tickers))
sub_idx <- match(cluster_manual, tickers)
mu_sub <- mu[sub_idx]; S_sub <- S[sub_idx, sub_idx, drop = FALSE]
tkr_sub <- tickers[sub_idx]; k_sub <- length(tkr_sub)
stopifnot(alpha_w * k_sub <= 1 + 1e-12, beta_w * k_sub >= 1 - 1e-12)
res_cluster <- run_mcabc(mu = mu_sub, S = S_sub, tickers = tkr_sub,
alpha = alpha_w, beta = beta_w,
d_min = k_sub, d_max = k_sub,
must_have_name = "", sigma = sigma_abc,
limit_factor = limit_factor, max_iter = max_iter, debug = FALSE)
P_cluster <- res_cluster$P
met_c <- portfolio_metrics(P_cluster, mu_sub, S_sub, rf)
best_c_id <- which.max(met_c$Sharpe)
w_best_sub <- as.numeric(P_cluster[best_c_id, ])
# KUNCI: paksa proyeksi ulang tepat di simplex [α,β] dan sum=1
w_best_sub <- project_box_simplex(
w = w_best_sub, L = alpha_w, U = beta_w, s = 1
)
# Jika masih ada nol karena numerik, aktifkan semuanya lalu proyeksi lagi
if (sum(w_best_sub > 1e-12) < length(w_best_sub)) {
w_best_sub[] <- pmax(w_best_sub, alpha_w)
w_best_sub <- project_box_simplex(w_best_sub, L = alpha_w, U = beta_w, s = 1)
}
stopifnot(all(w_best_sub > 0))
# Ekspansi ke universe penuh
w_cluster <- numeric(length(tickers))
names(w_cluster) <- tickers
w_cluster[tkr_sub] <- w_best_sub
# (B) Portofolio #2-#4 dari sekali run MOCv-ABC (k=3,4,5; Sharpe terbesar)
P_all <- mcabc_cov$P; k_count <- rowSums(P_all > 1e-8)
idx_k3 <- which(k_count == 3L); idx_k4 <- which(k_count == 4L); idx_k5 <- which(k_count == 5L)
best_k3 <- pick_best_by_k(P_all, mu, S, rf, idx_k3)
best_k4 <- pick_best_by_k(P_all, mu, S, rf, idx_k4)
best_k5 <- pick_best_by_k(P_all, mu, S, rf, idx_k5)
if (is.null(best_k3)) message("Tidak ada k=3 pada run ini—pertimbangkan menambah iterasi/seed.")
if (is.null(best_k4)) message("Tidak ada k=4 pada run ini—pertimbangkan menambah iterasi/seed.")
if (is.null(best_k5)) message("Tidak ada k=5 pada run ini—pertimbangkan menambah iterasi/seed.")
W_list <- list(
`Cluster` = w_cluster,
`ABC_k3` = if (!is.null(best_k3)) best_k3$w else rep(NA_real_, length(tickers)),
`ABC_k4` = if (!is.null(best_k4)) best_k4$w else rep(NA_real_, length(tickers)),
`ABC_k5` = if (!is.null(best_k5)) best_k5$w else rep(NA_real_, length(tickers))
)
W4 <- do.call(rbind, W_list); colnames(W4) <- tickers
M4 <- portfolio_metrics(W4, mu, S, rf)
opt4_summary <- data.frame(Portfolio = rownames(W4), M4, row.names = NULL)
print(opt4_summary)
## Portfolio ER SD Sharpe
## 1 Cluster 0.0592954 0.1862985 0.06331452
## 2 ABC_k3 0.2832434 0.3129769 0.75322950
## 3 ABC_k4 0.2877790 0.2947821 0.81510733
## 4 ABC_k5 0.2740621 0.2773238 0.81695868
# Return harian 4 portofolio (untuk VaR & backtest)
common <- intersect(colnames(rets_xts), colnames(W4))
rets_use <- rets_xts[, common, drop = FALSE]
R4_mat <- as.matrix(rets_use) %*% as.matrix(t(W4[, common, drop = FALSE]))
opt4_daily_log_returns <- xts::xts(R4_mat, order.by = index(rets_use))
colnames(opt4_daily_log_returns) <- rownames(W4)
# Scatterplot Risiko–Return (Plotly, label rapi)
text_positions <- c("top center","bottom center","top right","bottom right")
pos_vec <- rep(text_positions, length.out = nrow(opt4_summary))
y_offset <- (max(opt4_summary$ER) - min(opt4_summary$ER)) * 0.03
opt4_summary$label_y <- opt4_summary$ER + y_offset * c(1,-1,1,-1)
opt4_summary$label_hover <- paste0(
"<b>", opt4_summary$Portfolio, "</b><br>",
"Expected Return: ", sprintf("%.4f", opt4_summary$ER), "<br>",
"Risk (SD): ", sprintf("%.4f", opt4_summary$SD), "<br>",
"Sharpe Ratio: ", sprintf("%.4f", opt4_summary$Sharpe)
)
p_opt4_plotly_clean <- plot_ly(
data = opt4_summary, x = ~SD, y = ~ER, type = "scatter", mode = "markers+text",
text = ~Portfolio, textposition = pos_vec,
textfont = list(size = 12, color = "black", family = "Arial"),
marker = list(size = 14, color = ~Sharpe, colorscale = list(c(0,1), c("#d7191c","#1a9641")),
line = list(color = "black", width = 1)),
hovertemplate = opt4_summary$label_hover, showlegend = FALSE
) %>%
layout(
title = list(text = "Scatterplot Risiko–Return", x = 0.03, y = 0.97,
xanchor = "left", font = list(size = 18)),
xaxis = list(title = "Risiko (SD, annualized)", tickformat = ".2f",
gridcolor = "rgba(200,200,200,0.3)", zeroline = TRUE),
yaxis = list(title = "Expected Return (annualized)", tickformat = ".2f",
gridcolor = "rgba(200,200,200,0.3)", zeroline = TRUE),
margin = list(t = 90, r = 40, b = 60, l = 70),
paper_bgcolor = "white", plot_bgcolor = "white"
)
p_opt4_plotly_clean
## =======================================================
## TABEL SAHAM TERPILIH & BOBOT OPTIMAL (4 PORTOFOLIO)
## =======================================================
# Pastikan paket gt aktif
if (!requireNamespace("gt", quietly = TRUE)) install.packages("gt")
library(gt)
# Konversi W4 (4xN matrix) menjadi bentuk long agar mudah ditampilkan
w_long <- W4 |>
as.data.frame() |>
tibble::rownames_to_column("Portfolio") |>
tidyr::pivot_longer(-Portfolio, names_to = "Saham", values_to = "Bobot") |>
dplyr::filter(Bobot > 1e-6) |> # tampilkan hanya saham aktif
dplyr::arrange(Portfolio, desc(Bobot)) |>
dplyr::mutate(Bobot = round(Bobot, 4))
# Gabungkan dengan return & SD individual untuk informasi tambahan (opsional)
w_long <- w_long |>
left_join(stats_df |> select(ticker, mean_ann, sd_ann),
by = c("Saham" = "ticker")) |>
rename(`Return (ann.)` = mean_ann, `SD (ann.)` = sd_ann)
# Buat tabel gt
gt_port_weights <- w_long |>
gt(groupname_col = "Portfolio") |>
fmt_number(columns = c(Bobot, `Return (ann.)`, `SD (ann.)`), decimals = 4) |>
cols_label(
Saham = "Kode Saham",
Bobot = html("Bobot<br/>(Optimal)"),
`Return (ann.)` = html("Return<br/>(Annualized)"),
`SD (ann.)` = html("Risiko (SD<br/>Annualized)")
) |>
data_color(
columns = Bobot,
colors = scales::col_numeric(
palette = c("#f7fbff", "#6baed6", "#08306b"),
domain = NULL
)
) |>
tab_header(
title = md("**Saham Terpilih dan Bobot Optimal per Portofolio**"),
subtitle = md("Menampilkan kombinasi saham aktif dan proporsi bobot hasil optimasi MOCv-ABC & Cluster Manual")
) |>
tab_options(
table.font.size = px(13),
data_row.padding = px(4)
)
gt_port_weights
| Saham Terpilih dan Bobot Optimal per Portofolio |
| Menampilkan kombinasi saham aktif dan proporsi bobot hasil optimasi MOCv-ABC & Cluster Manual |
| Kode Saham |
Bobot
(Optimal) |
Return
(Annualized) |
Risiko (SD
Annualized) |
| ABC_k3 |
| MEDC |
0.4000 |
0.3524 |
0.4963 |
| SSIA |
0.4000 |
0.3688 |
0.5521 |
| TLKM |
0.2000 |
−0.0263 |
0.2756 |
| ABC_k4 |
| MEDC |
0.3737 |
0.3524 |
0.4963 |
| SSIA |
0.3282 |
0.3688 |
0.5521 |
| UNTR |
0.1981 |
0.1901 |
0.3281 |
| TLKM |
0.1000 |
−0.0263 |
0.2756 |
| ABC_k5 |
| MEDC |
0.3279 |
0.3524 |
0.4963 |
| SSIA |
0.3099 |
0.3688 |
0.5521 |
| UNTR |
0.1623 |
0.1901 |
0.3281 |
| PGAS |
0.1000 |
0.1601 |
0.3080 |
| TLKM |
0.1000 |
−0.0263 |
0.2756 |
| Cluster |
| PGAS |
0.3333 |
0.1601 |
0.3080 |
| BBRI |
0.3255 |
0.0542 |
0.2871 |
| TLKM |
0.1411 |
−0.0263 |
0.2756 |
| EXCL |
0.1000 |
−0.0157 |
0.3368 |
| MTEL |
0.1000 |
−0.0641 |
0.2778 |
## =============================================
## 9) EFFICIENT FRONTIER — overlay 4 portofolio terpilih
## Dua varian (A: dengan bound & wajib; B: hanya kardinalitas)
## Tiga plot per k (k=3,4,5) → total 6 plot
## =============================================
# Util metrik & kardinalitas
metrics_from_W <- function(W, mu, S, rf) {
ER <- as.numeric(W %*% mu)
VAR <- apply(W, 1, function(w) as.numeric(t(w) %*% S %*% w))
SD <- sqrt(VAR)
Sharpe <- (ER - rf) / pmax(SD, 1e-12)
data.frame(ER = ER, SD = SD, Sharpe = Sharpe)
}
k_card <- function(W, eps = 1e-8) rowSums(W > eps)
# Layer overlay untuk 4 portofolio (highlight k yang sama)
overlay_layer <- function(p, overlay_df, k_focus) {
overlay_df$alpha_pt <- ifelse(overlay_df$k == k_focus, 1, 0.35)
overlay_df$size_pt <- ifelse(overlay_df$k == k_focus, 3.8, 3.0)
overlay_df$stroke <- ifelse(overlay_df$k == k_focus, 1.2, 0.6)
p +
geom_point(data = overlay_df, aes(SD, ER), shape = 21, fill = "skyblue", color = "black",
size = overlay_df$size_pt, stroke = overlay_df$stroke, alpha = overlay_df$alpha_pt) +
ggrepel::geom_text_repel(
data = overlay_df, aes(SD, ER, label = Portfolio),
size = 3.0, max.overlaps = 50, seed = 42, box.padding = 0.25, point.padding = 0.3
)
}
pal_grad <- colorRampPalette(c("#d7191c","#fdae61","lightgreen","#1a9641"))(256)
k_overlay <- k_card(W4)
overlay_df <- cbind(opt4_summary, k = k_overlay) # kolom: Portfolio, ER, SD, Sharpe, k
# Varian A: batas bobot [alpha_w, beta_w], sum=1, PGAS wajib
set.seed(123)
d_all <- length(tickers)
m_idx <- match("PGAS", tickers); if (is.na(m_idx)) m_idx <- NULL
rand_weights_A <- function(k, n = 600L) {
stopifnot(alpha_w * k <= 1 + 1e-12, beta_w * k >= 1 - 1e-12)
W <- matrix(NA_real_, nrow = n, ncol = d_all); colnames(W) <- tickers
for (i in 1:n) {
w0 <- runif(d_all)
W[i, ] <- clean_portfolio(w0, d_all, alpha_w, beta_w, d_min = k, d_max = k, must_have_idx_or_NULL = m_idx)
}
W
}
make_frontier_plot_A <- function(k, xlim, ylim) {
W_rand <- rand_weights_A(k, n = 600)
df_rand <- metrics_from_W(W_rand, mu, S, rf)
p <- ggplot(df_rand, aes(SD, ER, color = Sharpe)) +
coord_cartesian(xlim = xlim, ylim = ylim) +
geom_point(alpha = 0.85, size = 1.8) +
scale_color_gradientn(colors = pal_grad) +
labs(title = paste0("Efficient Frontier — Varian A (k = ", k, ")"),
subtitle = "Random portfolios dengan batas bobot [α,β] & PGAS wajib",
x = "Risiko (SD, annualized)", y = "Expected Return (annualized)", color = "Sharpe") +
theme_minimal(base_size = 12)
overlay_layer(p, overlay_df, k_focus = k) |> plotly::ggplotly()
}
pA_k3 <- make_frontier_plot_A(3, xlim = c(0.15, 0.51), ylim = c(-0.3, 0.4))
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
pA_k4 <- make_frontier_plot_A(4, xlim = c(0.15, 0.4), ylim = c(-0.3, 0.35))
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
pA_k5 <- make_frontier_plot_A(5, xlim = c(0.15, 0.38), ylim = c(-0.2, 0.35))
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
# Varian B: hanya kardinalitas = k (tanpa batas bobot & tanpa wajib)
set.seed(123)
rand_weights_B <- function(k, n = 600L) {
W <- matrix(0, nrow = n, ncol = d_all); colnames(W) <- tickers
for (i in 1:n) {
idx <- sample.int(d_all, k, replace = FALSE)
w <- runif(k); w <- w / sum(w)
W[i, idx] <- w
}
W
}
make_frontier_plot_B <- function(k, xlim, ylim) {
W_rand <- rand_weights_B(k, n = 600)
df_rand <- metrics_from_W(W_rand, mu, S, rf)
p <- ggplot(df_rand, aes(SD, ER, color = Sharpe)) +
coord_cartesian(xlim = xlim, ylim = ylim) +
geom_point(alpha = 0.85, size = 1.8) +
scale_color_gradientn(colors = pal_grad) +
labs(title = paste0("Efficient Frontier — Varian B (k = ", k, ")"),
subtitle = "Random portfolios kardinalitas saja (tanpa batas bobot & wajib)",
x = "Risiko (SD, annualized)", y = "Expected Return (annualized)", color = "Sharpe") +
theme_minimal(base_size = 12)
overlay_layer(p, overlay_df, k_focus = k) |> plotly::ggplotly()
}
pB_k3 <- make_frontier_plot_B(3, xlim = c(0.15, 0.51), ylim = c(-0.3, 0.4))
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
pB_k4 <- make_frontier_plot_B(4, xlim = c(0.15, 0.4), ylim = c(-0.3, 0.35))
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
pB_k5 <- make_frontier_plot_B(5, xlim = c(0.15, 0.38), ylim = c(-0.2, 0.35))
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomTextRepel() has yet to be implemented in plotly.
## If you'd like to see this geom implemented,
## Please open an issue with your example code at
## https://github.com/ropensci/plotly/issues
# Tampilkan (opsional)
pA_k3; pB_k3
pA_k4; pB_k4
pA_k5; pB_k5
## =============================================
## 10) VAR CORNISH–FISHER (harian) UNTUK 4 PORTOFOLIO
## =============================================
# Statistik harian
Rport_daily <- as.matrix(opt4_daily_log_returns) # [T x 4]
port_names <- colnames(Rport_daily)
mu_d <- colMeans(Rport_daily, na.rm = TRUE)
sd_d <- apply(Rport_daily, 2, sd, na.rm = TRUE)
sk_d <- apply(Rport_daily, 2, function(x) moments::skewness(x, na.rm = TRUE))
kurt_d <- apply(Rport_daily, 2, function(x) moments::kurtosis(x, na.rm = TRUE))
ex_kurt_d <- kurt_d - 3
stats_top4_daily <- data.frame(
Portfolio = port_names, mu_d = mu_d, sd_d = sd_d, skew = sk_d, kurt = kurt_d, ex_kurt = ex_kurt_d
)
# CF quantile & VaR
cf_transform <- function(z, skew, ex_kurt) {
z + (1/6)*(z^2 - 1)*skew + (1/24)*(z^3 - 3*z)*ex_kurt + (1/36)*(2*z^3 - 5*z)*(skew^2)
}
cf_quantile <- function(alpha, skew, ex_kurt) {
a <- qnorm(alpha)
a + (skew/6)*(a^2 - 1) + (ex_kurt/24)*(a^3 - 3*a) - ((skew^2)/36)*(2*a^3 - 5*a)
}
var_cf_vector <- function(mu, sd, skew, ex_kurt, alpha, V0 = 1, hp = 1) {
a_cf <- cf_quantile(alpha, skew, ex_kurt)
VaR_CF <- -V0 * (mu - a_cf * sd) * hp
list(VaR_CF = VaR_CF, a_cf = a_cf)
}
res_list <- mapply(
var_cf_vector, mu = mu_d, sd = sd_d, skew = sk_d, ex_kurt = ex_kurt_d,
MoreArgs = list(alpha = alpha, V0 = V0, hp = hp), SIMPLIFY = FALSE
)
VaR_CF_vals <- sapply(res_list, function(x) x$VaR_CF)
a_cf_vals <- sapply(res_list, function(x) x$a_cf)
VaR_CF_summary <- data.frame(
Portfolio = port_names, CL = cl, alpha = alpha,
mu_d = mu_d, sd_d = sd_d, skew = sk_d, kurt = kurt_d, excess_k = ex_kurt_d,
a_cf = a_cf_vals, VaR_CF = VaR_CF_vals
)
VaR_CF_summary
## Portfolio CL alpha mu_d sd_d skew kurt
## Cluster Cluster 0.95 0.05 0.0002352992 0.01173571 -0.356812016 7.063192
## ABC_k3 ABC_k3 0.95 0.05 0.0011239817 0.01971569 0.304637208 5.251243
## ABC_k4 ABC_k4 0.95 0.05 0.0011419803 0.01856953 0.039122898 5.719642
## ABC_k5 ABC_k5 0.95 0.05 0.0010875478 0.01746975 -0.001683576 5.806608
## excess_k a_cf VaR_CF
## Cluster 4.063192 -1.661890 -0.01973876
## ABC_k3 2.251243 -1.511083 -0.03091603
## ABC_k4 2.719642 -1.578820 -0.03045991
## ABC_k5 2.806608 -1.588693 -0.02884162
## =============================================
# 11) BACKTEST VAR: KUPIEC + PLOT INTERAKTIF & TABEL GT
## =============================================
# Kupiec (unconditional coverage) — breach jika return < VaR (VaR negatif)
rets_mat <- as.matrix(opt4_daily_log_returns)
port_names <- colnames(rets_mat)
VaR_compare_vec <- VaR_CF_summary$VaR_CF[seq_along(port_names)]
kupiec_test <- function(returns, VaR, alpha = 0.05) {
x <- as.numeric(returns)
v <- as.numeric(VaR); if (length(v) == 1) v <- rep(v, length(x))
I <- as.integer(x < v) # breach jika return < VaR (ingat VaR negatif)
I <- I[is.finite(I)]
n <- length(I); xB <- sum(I)
if (n == 0) return(data.frame(n = 0, x = NA, prop = NA, LR_uc = NA, p_value = NA))
p_hat <- pmax(pmin(xB / n, 1 - 1e-12), 1e-12)
LR_uc <- -2 * ( (n - xB)*log(1 - alpha) + xB*log(alpha) - (n - xB)*log(1 - p_hat) - xB*log(p_hat) )
p_uc <- 1 - pchisq(LR_uc, df = 1)
data.frame(n = n, x = xB, prop = xB/n, LR_uc = LR_uc, p_value = p_uc)
}
results_list <- lapply(seq_along(port_names), function(i) {
ku <- kupiec_test(rets_mat[, i], VaR_compare_vec[i], alpha = alpha)
data.frame(Portfolio = port_names[i], CL = 1 - alpha,
Violations = ku$x, ViolationRate = ku$prop,
Kupiec_LR = ku$LR_uc, Kupiec_p = ku$p_value)
})
VaR_backtest_table <- do.call(rbind, results_list)
VaR_backtest_table
## Portfolio CL Violations ViolationRate Kupiec_LR Kupiec_p
## 1 Cluster 0.95 38 0.04255319 1.095613 0.29523094
## 2 ABC_k3 0.95 34 0.03807391 2.903565 0.08838391
## 3 ABC_k4 0.95 37 0.04143337 1.461484 0.22669410
## 4 ABC_k5 0.95 36 0.04031355 1.883753 0.16990822
# Data untuk plot/gt (built from scratch agar konsisten dengan Kupiec)
if (exists("plot_df", inherits = FALSE)) rm(plot_df)
ret_df <- opt4_daily_log_returns |>
as.data.frame() |> tibble::rownames_to_column("date") |>
dplyr::mutate(date = as.Date(date)) |>
tidyr::pivot_longer(-date, names_to = "Portfolio", values_to = "ret")
var_df <- VaR_CF_summary |> dplyr::transmute(Portfolio, VaR_line = VaR_CF) # NEGATIF
plot_df <- dplyr::left_join(ret_df, var_df, by = "Portfolio") |>
dplyr::mutate(breach = ret < VaR_line, breach_n = ave(breach, Portfolio, FUN = cumsum))
# validasi jumlah breach plot_df vs hasil kupiec
check_cnt <- plot_df |>
dplyr::group_by(Portfolio) |>
dplyr::summarise(breach = sum(breach), .groups = "drop")
print(check_cnt)
## # A tibble: 4 × 2
## Portfolio breach
## <chr> <int>
## 1 ABC_k3 34
## 2 ABC_k4 37
## 3 ABC_k5 36
## 4 Cluster 38
# Plot pelanggaran per-portofolio (plotly)
plot_var_breach <- function(port_name) {
dat <- plot_df |> filter(Portfolio == port_name) |> arrange(date)
var_line <- unique(dat$VaR_line); stopifnot(length(var_line) == 1)
p <- plot_ly()
p <- add_bars(p, data = dat, x = ~date, y = ~ret, name = "Return",
legendgroup = "g1",
marker = list(line = list(width = 0)),
hovertemplate = "<b>%{x}</b><br>Return: %{y:.4f}<extra>Return</extra>")
p <- add_lines(p, data = dat, x = ~date, y = ~VaR_line, name = "VaR (ambang)",
legendgroup = "g2", line = list(width = 2), hoverinfo = "skip")
p <- add_markers(p, data = subset(dat, breach), x = ~date, y = ~ret, name = "Pelanggaran",
legendgroup = "g3", marker = list(size = 8, symbol = "circle-open-dot", color = "red"),
customdata = ~breach_n,
hovertemplate = paste0("<b>%{x}</b><br>",
"Return: %{y:.4f}<br>",
"Ambang VaR: ", sprintf("%.4f", var_line), "<br>",
"<b>PELANGGARAN</b> # %{customdata}<extra>Pelanggaran</extra>"))
layout(p,
title = list(text = paste0("Plot Pelanggaran VaR — ", port_name),
x = 0.02, y = 0.98, xanchor = "left", yanchor = "top"),
margin = list(t = 90, r = 20, b = 40, l = 55),
legend = list(orientation = "h", x = 0.02, y = 1.05, xanchor = "left", yanchor = "bottom",
bgcolor = "rgba(255,255,255,0.6)"),
xaxis = list(title = NULL, automargin = TRUE, rangeslider = list(visible = TRUE)),
yaxis = list(title = "Return harian", automargin = TRUE),
barmode = "overlay")
}
plot_var_breach("Cluster")
plot_var_breach("ABC_k3")
plot_var_breach("ABC_k4")
plot_var_breach("ABC_k5")
# Tabel gt: daftar tanggal pelanggaran & ringkasannya
var_breach_table_gt <- function(port_name, digits = 4) {
dat_port <- plot_df |> dplyr::filter(Portfolio == port_name)
stopifnot(nrow(dat_port) > 0)
n_total <- nrow(dat_port)
n_breach <- sum(dat_port$breach, na.rm = TRUE)
rate <- n_breach / n_total
df_breach <- dat_port |>
dplyr::filter(breach) |>
arrange(date) |>
transmute(
Tanggal = date,
`Return` = ret,
`VaR` = VaR_line,
`Selisih`= ret - VaR_line,
`Ke-` = breach_n
)
if (nrow(df_breach) == 0) {
df_breach <- tibble::tibble(Tanggal = as.Date(character()),
`Return` = numeric(), `VaR` = numeric(),
`Selisih` = numeric(), `Ke-` = integer())
}
df_breach |>
gt() |>
fmt_number(columns = c(`Return`,`VaR`,`Selisih`), decimals = digits) |>
data_color(columns = c(`Return`), colors = col_bin(palette = c("#CC3232","#2DC937"), bins = c(-Inf,0,Inf))) |>
tab_header(
title = md(paste0("**Pelanggaran VaR ", port_name, "**")),
subtitle = md(paste0("Total observasi: **", n_total,
"** · Jumlah pelanggaran: **", n_breach,
"** · Violation rate: **", scales::percent(rate, accuracy = 0.01), "**"))
) |>
cols_label(`Ke-` = "Pelanggaran ke-") |>
tab_source_note(md("Catatan: Pelanggaran terjadi saat **return < VaR** (ambang return minimum harian)."))
}
var_breach_table_gt("Cluster")
| Pelanggaran VaR Cluster |
| Total observasi: 893 · Jumlah pelanggaran: 38 · Violation rate: 4.26% |
| Tanggal |
Return |
VaR |
Selisih |
Pelanggaran ke- |
| 2022-03-08 |
−0.0200 |
−0.0197 |
−0.0003 |
1 |
| 2022-03-10 |
−0.0230 |
−0.0197 |
−0.0033 |
2 |
| 2022-05-09 |
−0.0395 |
−0.0197 |
−0.0198 |
3 |
| 2022-05-12 |
−0.0329 |
−0.0197 |
−0.0132 |
4 |
| 2022-07-04 |
−0.0261 |
−0.0197 |
−0.0063 |
5 |
| 2022-07-06 |
−0.0200 |
−0.0197 |
−0.0002 |
6 |
| 2022-09-07 |
−0.0216 |
−0.0197 |
−0.0019 |
7 |
| 2023-01-04 |
−0.0255 |
−0.0197 |
−0.0057 |
8 |
| 2023-01-05 |
−0.0348 |
−0.0197 |
−0.0151 |
9 |
| 2023-03-14 |
−0.0304 |
−0.0197 |
−0.0107 |
10 |
| 2023-03-20 |
−0.0215 |
−0.0197 |
−0.0017 |
11 |
| 2023-09-07 |
−0.0274 |
−0.0197 |
−0.0077 |
12 |
| 2023-10-26 |
−0.0320 |
−0.0197 |
−0.0123 |
13 |
| 2023-10-30 |
−0.0256 |
−0.0197 |
−0.0059 |
14 |
| 2023-11-01 |
−0.0234 |
−0.0197 |
−0.0037 |
15 |
| 2024-04-16 |
−0.0341 |
−0.0197 |
−0.0143 |
16 |
| 2024-04-19 |
−0.0243 |
−0.0197 |
−0.0046 |
17 |
| 2024-04-26 |
−0.0325 |
−0.0197 |
−0.0127 |
18 |
| 2024-05-27 |
−0.0225 |
−0.0197 |
−0.0028 |
19 |
| 2024-06-14 |
−0.0278 |
−0.0197 |
−0.0080 |
20 |
| 2024-08-05 |
−0.0332 |
−0.0197 |
−0.0134 |
21 |
| 2024-08-27 |
−0.0280 |
−0.0197 |
−0.0083 |
22 |
| 2024-09-30 |
−0.0222 |
−0.0197 |
−0.0025 |
23 |
| 2024-10-02 |
−0.0226 |
−0.0197 |
−0.0028 |
24 |
| 2024-11-11 |
−0.0201 |
−0.0197 |
−0.0004 |
25 |
| 2024-12-17 |
−0.0209 |
−0.0197 |
−0.0012 |
26 |
| 2025-01-07 |
−0.0249 |
−0.0197 |
−0.0051 |
27 |
| 2025-01-14 |
−0.0204 |
−0.0197 |
−0.0007 |
28 |
| 2025-02-06 |
−0.0311 |
−0.0197 |
−0.0114 |
29 |
| 2025-02-25 |
−0.0265 |
−0.0197 |
−0.0068 |
30 |
| 2025-02-27 |
−0.0355 |
−0.0197 |
−0.0158 |
31 |
| 2025-02-28 |
−0.0434 |
−0.0197 |
−0.0237 |
32 |
| 2025-04-08 |
−0.0846 |
−0.0197 |
−0.0649 |
33 |
| 2025-05-08 |
−0.0224 |
−0.0197 |
−0.0027 |
34 |
| 2025-06-02 |
−0.0201 |
−0.0197 |
−0.0004 |
35 |
| 2025-06-19 |
−0.0214 |
−0.0197 |
−0.0017 |
36 |
| 2025-08-27 |
−0.0204 |
−0.0197 |
−0.0006 |
37 |
| 2025-09-08 |
−0.0226 |
−0.0197 |
−0.0029 |
38 |
| Catatan: Pelanggaran terjadi saat return < VaR (ambang return minimum harian). |
var_breach_table_gt("ABC_k3")
| Pelanggaran VaR ABC_k3 |
| Total observasi: 893 · Jumlah pelanggaran: 34 · Violation rate: 3.81% |
| Tanggal |
Return |
VaR |
Selisih |
Pelanggaran ke- |
| 2022-01-25 |
−0.0371 |
−0.0309 |
−0.0061 |
1 |
| 2022-02-16 |
−0.0351 |
−0.0309 |
−0.0042 |
2 |
| 2022-02-21 |
−0.0377 |
−0.0309 |
−0.0068 |
3 |
| 2022-07-01 |
−0.0504 |
−0.0309 |
−0.0195 |
4 |
| 2022-09-16 |
−0.0313 |
−0.0309 |
−0.0004 |
5 |
| 2022-10-18 |
−0.0327 |
−0.0309 |
−0.0018 |
6 |
| 2022-11-18 |
−0.0360 |
−0.0309 |
−0.0051 |
7 |
| 2023-01-05 |
−0.0375 |
−0.0309 |
−0.0066 |
8 |
| 2023-02-02 |
−0.0380 |
−0.0309 |
−0.0071 |
9 |
| 2023-02-09 |
−0.0391 |
−0.0309 |
−0.0082 |
10 |
| 2023-02-10 |
−0.0387 |
−0.0309 |
−0.0078 |
11 |
| 2023-02-17 |
−0.0415 |
−0.0309 |
−0.0106 |
12 |
| 2023-03-06 |
−0.0370 |
−0.0309 |
−0.0060 |
13 |
| 2023-03-14 |
−0.0474 |
−0.0309 |
−0.0165 |
14 |
| 2023-03-16 |
−0.0355 |
−0.0309 |
−0.0045 |
15 |
| 2023-03-20 |
−0.0412 |
−0.0309 |
−0.0103 |
16 |
| 2023-05-02 |
−0.0406 |
−0.0309 |
−0.0097 |
17 |
| 2023-08-01 |
−0.0416 |
−0.0309 |
−0.0106 |
18 |
| 2023-10-19 |
−0.0337 |
−0.0309 |
−0.0028 |
19 |
| 2023-10-30 |
−0.0544 |
−0.0309 |
−0.0235 |
20 |
| 2023-11-08 |
−0.0408 |
−0.0309 |
−0.0099 |
21 |
| 2024-04-19 |
−0.0319 |
−0.0309 |
−0.0010 |
22 |
| 2024-06-14 |
−0.0377 |
−0.0309 |
−0.0068 |
23 |
| 2024-08-05 |
−0.0763 |
−0.0309 |
−0.0454 |
24 |
| 2024-09-30 |
−0.0366 |
−0.0309 |
−0.0057 |
25 |
| 2024-11-18 |
−0.0335 |
−0.0309 |
−0.0026 |
26 |
| 2025-01-02 |
−0.0668 |
−0.0309 |
−0.0358 |
27 |
| 2025-01-30 |
−0.0398 |
−0.0309 |
−0.0089 |
28 |
| 2025-02-28 |
−0.0327 |
−0.0309 |
−0.0018 |
29 |
| 2025-03-18 |
−0.0406 |
−0.0309 |
−0.0097 |
30 |
| 2025-04-08 |
−0.0947 |
−0.0309 |
−0.0638 |
31 |
| 2025-06-24 |
−0.0472 |
−0.0309 |
−0.0163 |
32 |
| 2025-07-22 |
−0.0451 |
−0.0309 |
−0.0142 |
33 |
| 2025-09-08 |
−0.0404 |
−0.0309 |
−0.0095 |
34 |
| Catatan: Pelanggaran terjadi saat return < VaR (ambang return minimum harian). |
var_breach_table_gt("ABC_k4")
| Pelanggaran VaR ABC_k4 |
| Total observasi: 893 · Jumlah pelanggaran: 37 · Violation rate: 4.14% |
| Tanggal |
Return |
VaR |
Selisih |
Pelanggaran ke- |
| 2022-01-25 |
−0.0340 |
−0.0305 |
−0.0036 |
1 |
| 2022-02-16 |
−0.0338 |
−0.0305 |
−0.0033 |
2 |
| 2022-02-21 |
−0.0353 |
−0.0305 |
−0.0048 |
3 |
| 2022-03-14 |
−0.0388 |
−0.0305 |
−0.0083 |
4 |
| 2022-06-13 |
−0.0351 |
−0.0305 |
−0.0046 |
5 |
| 2022-07-01 |
−0.0584 |
−0.0305 |
−0.0279 |
6 |
| 2022-10-18 |
−0.0325 |
−0.0305 |
−0.0020 |
7 |
| 2022-11-10 |
−0.0329 |
−0.0305 |
−0.0025 |
8 |
| 2023-01-04 |
−0.0333 |
−0.0305 |
−0.0029 |
9 |
| 2023-02-02 |
−0.0378 |
−0.0305 |
−0.0073 |
10 |
| 2023-02-09 |
−0.0342 |
−0.0305 |
−0.0037 |
11 |
| 2023-02-10 |
−0.0337 |
−0.0305 |
−0.0033 |
12 |
| 2023-02-17 |
−0.0355 |
−0.0305 |
−0.0050 |
13 |
| 2023-03-06 |
−0.0350 |
−0.0305 |
−0.0046 |
14 |
| 2023-03-14 |
−0.0483 |
−0.0305 |
−0.0179 |
15 |
| 2023-03-16 |
−0.0354 |
−0.0305 |
−0.0049 |
16 |
| 2023-03-20 |
−0.0410 |
−0.0305 |
−0.0105 |
17 |
| 2023-05-02 |
−0.0489 |
−0.0305 |
−0.0184 |
18 |
| 2023-08-01 |
−0.0377 |
−0.0305 |
−0.0073 |
19 |
| 2023-10-19 |
−0.0308 |
−0.0305 |
−0.0003 |
20 |
| 2023-10-23 |
−0.0330 |
−0.0305 |
−0.0026 |
21 |
| 2023-10-30 |
−0.0575 |
−0.0305 |
−0.0271 |
22 |
| 2023-11-01 |
−0.0367 |
−0.0305 |
−0.0063 |
23 |
| 2023-11-08 |
−0.0401 |
−0.0305 |
−0.0097 |
24 |
| 2023-11-17 |
−0.0309 |
−0.0305 |
−0.0004 |
25 |
| 2024-06-14 |
−0.0357 |
−0.0305 |
−0.0052 |
26 |
| 2024-08-05 |
−0.0734 |
−0.0305 |
−0.0429 |
27 |
| 2024-09-30 |
−0.0340 |
−0.0305 |
−0.0035 |
28 |
| 2024-11-18 |
−0.0361 |
−0.0305 |
−0.0057 |
29 |
| 2025-01-02 |
−0.0618 |
−0.0305 |
−0.0314 |
30 |
| 2025-01-21 |
−0.0309 |
−0.0305 |
−0.0004 |
31 |
| 2025-01-30 |
−0.0388 |
−0.0305 |
−0.0083 |
32 |
| 2025-02-28 |
−0.0338 |
−0.0305 |
−0.0033 |
33 |
| 2025-03-18 |
−0.0405 |
−0.0305 |
−0.0100 |
34 |
| 2025-04-08 |
−0.1083 |
−0.0305 |
−0.0779 |
35 |
| 2025-06-24 |
−0.0407 |
−0.0305 |
−0.0102 |
36 |
| 2025-07-22 |
−0.0382 |
−0.0305 |
−0.0077 |
37 |
| Catatan: Pelanggaran terjadi saat return < VaR (ambang return minimum harian). |
var_breach_table_gt("ABC_k5")
| Pelanggaran VaR ABC_k5 |
| Total observasi: 893 · Jumlah pelanggaran: 36 · Violation rate: 4.03% |
| Tanggal |
Return |
VaR |
Selisih |
Pelanggaran ke- |
| 2022-01-25 |
−0.0317 |
−0.0288 |
−0.0029 |
1 |
| 2022-02-21 |
−0.0329 |
−0.0288 |
−0.0041 |
2 |
| 2022-03-14 |
−0.0314 |
−0.0288 |
−0.0026 |
3 |
| 2022-06-13 |
−0.0328 |
−0.0288 |
−0.0039 |
4 |
| 2022-07-01 |
−0.0511 |
−0.0288 |
−0.0223 |
5 |
| 2022-09-07 |
−0.0296 |
−0.0288 |
−0.0008 |
6 |
| 2022-10-18 |
−0.0293 |
−0.0288 |
−0.0005 |
7 |
| 2022-11-10 |
−0.0297 |
−0.0288 |
−0.0008 |
8 |
| 2023-01-04 |
−0.0333 |
−0.0288 |
−0.0045 |
9 |
| 2023-01-05 |
−0.0327 |
−0.0288 |
−0.0038 |
10 |
| 2023-02-02 |
−0.0322 |
−0.0288 |
−0.0033 |
11 |
| 2023-02-09 |
−0.0310 |
−0.0288 |
−0.0021 |
12 |
| 2023-02-10 |
−0.0310 |
−0.0288 |
−0.0022 |
13 |
| 2023-02-17 |
−0.0318 |
−0.0288 |
−0.0029 |
14 |
| 2023-03-06 |
−0.0305 |
−0.0288 |
−0.0017 |
15 |
| 2023-03-14 |
−0.0481 |
−0.0288 |
−0.0192 |
16 |
| 2023-03-16 |
−0.0352 |
−0.0288 |
−0.0064 |
17 |
| 2023-03-20 |
−0.0406 |
−0.0288 |
−0.0117 |
18 |
| 2023-05-02 |
−0.0460 |
−0.0288 |
−0.0172 |
19 |
| 2023-08-01 |
−0.0324 |
−0.0288 |
−0.0036 |
20 |
| 2023-10-19 |
−0.0295 |
−0.0288 |
−0.0007 |
21 |
| 2023-10-23 |
−0.0307 |
−0.0288 |
−0.0019 |
22 |
| 2023-10-30 |
−0.0553 |
−0.0288 |
−0.0265 |
23 |
| 2023-11-01 |
−0.0340 |
−0.0288 |
−0.0052 |
24 |
| 2023-11-08 |
−0.0419 |
−0.0288 |
−0.0130 |
25 |
| 2024-06-14 |
−0.0343 |
−0.0288 |
−0.0055 |
26 |
| 2024-08-05 |
−0.0704 |
−0.0288 |
−0.0415 |
27 |
| 2024-09-30 |
−0.0334 |
−0.0288 |
−0.0046 |
28 |
| 2024-11-18 |
−0.0335 |
−0.0288 |
−0.0047 |
29 |
| 2025-01-02 |
−0.0566 |
−0.0288 |
−0.0277 |
30 |
| 2025-01-30 |
−0.0360 |
−0.0288 |
−0.0072 |
31 |
| 2025-02-28 |
−0.0320 |
−0.0288 |
−0.0031 |
32 |
| 2025-03-18 |
−0.0381 |
−0.0288 |
−0.0093 |
33 |
| 2025-04-08 |
−0.1046 |
−0.0288 |
−0.0758 |
34 |
| 2025-06-24 |
−0.0382 |
−0.0288 |
−0.0094 |
35 |
| 2025-07-22 |
−0.0371 |
−0.0288 |
−0.0082 |
36 |
| Catatan: Pelanggaran terjadi saat return < VaR (ambang return minimum harian). |