Optimasi Portofolio
Fedhilla Fieldhathama
2025-06-12
Menampilkan packages untuk mendukung komputasi
#packages
library(quantmod)## Warning: package 'quantmod' was built under R version 4.3.3## Loading required package: xts## Warning: package 'xts' was built under R version 4.3.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.3.1##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: TTR## Warning: package 'TTR' was built under R version 4.3.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoolibrary(purrr)## Warning: package 'purrr' was built under R version 4.3.1library(quadprog)## Warning: package 'quadprog' was built under R version 4.3.1library(PerformanceAnalytics)## Warning: package 'PerformanceAnalytics' was built under R version 4.3.3##
## Attaching package: 'PerformanceAnalytics'## The following object is masked from 'package:graphics':
##
## legendlibrary(dplyr)## Warning: package 'dplyr' was built under R version 4.3.1##
## ######################### Warning from 'xts' package ##########################
## # #
## # The dplyr lag() function breaks how base R's lag() function is supposed to #
## # work, which breaks lag(my_xts). Calls to lag(my_xts) that you type or #
## # source() into this session won't work correctly. #
## # #
## # Use stats::lag() to make sure you're not using dplyr::lag(), or you can add #
## # conflictRules('dplyr', exclude = 'lag') to your .Rprofile to stop #
## # dplyr from breaking base R's lag() function. #
## # #
## # Code in packages is not affected. It's protected by R's namespace mechanism #
## # Set `options(xts.warn_dplyr_breaks_lag = FALSE)` to suppress this warning. #
## # #
## #################################################################################
## Attaching package: 'dplyr'## The following objects are masked from 'package:xts':
##
## first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(NMOF)## Warning: package 'NMOF' was built under R version 4.3.3library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.3library(tseries)## Warning: package 'tseries' was built under R version 4.3.3library(psych)## Warning: package 'psych' was built under R version 4.3.3##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(mco)## Warning: package 'mco' was built under R version 4.3.3PERSIAPAN DATA YANG AKAN DIGUNAKAN
Load Data saham saham yang ada di Economic30
list saham yang termasuk pada Economic30:
1. PT Ace Hardware Indonesia Tbk (ACES)
2. PT Amman Mineral Internasional Tbk (AMMN)
3. PT Aneka Tambang Tbk (ANTM)
4. PT Bank Jago Tbk (ARTO)
5. PT Bank Central Asia Tbk (BBCA)
6. PT Bank Negara Indonesia Tbk (BBNI)
7. PT Bank Rakyat Indonesia Tbk (BBRI)
8. PT Bank Tabungan Negara Tbk (BBTN)
9. PT Bank Mandiri Tbk (BMRI)
10. PT Bank Syariah Indonesia Tbk (BRIS)
11. PT Bumi Resources Minerals Tbk (BRMS)
12. PT Barito Pacific Tbk (BRPT)
13. PT Ciputra Development Tbk (CTRA)
14. PT ESSA Industries Indonesia Tbk (ESSA)
15. PT MD Pictures Tbk (FILM)
16. PT Vale Indonesia Tbk (INCO)
17. PT Indah Kiat Pulp and Paper Tbk (INKP)
18. PT Indocement Tunggal Prakarsa Tbk (INTP)
19. PT Mitra Adiperkasa Tbk (MAPA)
20. PT MAP Aktif Adiperkasa Tbk (MAPA)
21. PT Merdeka Copper Gold Tbk (MDKA)
22. PT Media Nusantara Citra Tbk (MNCN)
23. PT Panin Financial Tbk (PNLF)
24. PT Surya Citra Media Tbk (SCMA)
25. PT Semen Indonesia Tbk (SMGR)
26. PT Saratoga Investama Sedaya Tbk (SRTG)
27. PT Chandra Asri Pacific Tbk (TPIA)
28. PT Petrindo Jaya Kreasi Tbk (CUAN)
29. PT MBMA Tbk (MBMA)
30. PT Pakuwon Jati Tbk (PWON)
saham <- c(
"ACES.JK", "AMMN.JK", "ANTM.JK", "ARTO.JK", "BBCA.JK", "BBNI.JK", "BBRI.JK", "BBTN.JK", "BMRI.JK", "BRIS.JK", "BRMS.JK", "BRPT.JK", "CTRA.JK", "ESSA.JK", "FILM.JK", "INCO.JK", "INKP.JK", "INTP.JK", "MAPA.JK","MAPI.JK", "MBMA.JK", "MDKA.JK", "MNCN.JK", "PANI.JK", "PNLF.JK", "PWON.JK", "SMGR.JK", "SRTG.JK", "TPIA.JK", "SMRA.JK"
)
getSymbols(saham, src = "yahoo", from = "2020-01-01", to = "2025-01-01")## [1] "ACES.JK" "AMMN.JK" "ANTM.JK" "ARTO.JK" "BBCA.JK" "BBNI.JK" "BBRI.JK"
## [8] "BBTN.JK" "BMRI.JK" "BRIS.JK" "BRMS.JK" "BRPT.JK" "CTRA.JK" "ESSA.JK"
## [15] "FILM.JK" "INCO.JK" "INKP.JK" "INTP.JK" "MAPA.JK" "MAPI.JK" "MBMA.JK"
## [22] "MDKA.JK" "MNCN.JK" "PANI.JK" "PNLF.JK" "PWON.JK" "SMGR.JK" "SRTG.JK"
## [29] "TPIA.JK" "SMRA.JK"getSymbols("^JKSE", from = "2020-01-01", to = Sys.Date(), src = "yahoo")## [1] "JKSE"Setelah data saham sudah berhasil untuk diimport, lalu selanjutnya menghitung nilai return tiap data. Jika data memiliki nilai missing value maka drop kolom/hapus saham yang memiliki nilai missing value. karena AMMN.JK memiliki nilai missing value sebanyak 853 dan MBMA.JK memiliki nilai missing value sebanyak 807 baris maka AMMN.JK dan MBMA.JK dihapus dalam daftar analisis kita.
Langkah selanjutnya yakni, menghapus semua saham yang memiliki rata-rata return negatif. Lalu selanjutnya menghitung semua matriks korelasi pada tiap saham yang tersisa.
# Hitung return harian untuk tiap saham
returns_data <- do.call(merge, lapply(saham, function(x) dailyReturn(Cl(get(x)))))
colnames(returns_data) <- saham
# Bersihkan kolom dengan NA
cat("berikut merupakan tampilan jumlah data yang kosong pada masing-masing saham: \n")## berikut merupakan tampilan jumlah data yang kosong pada masing-masing saham:colSums(is.na(returns_data))## ACES.JK AMMN.JK ANTM.JK ARTO.JK BBCA.JK BBNI.JK BBRI.JK BBTN.JK BMRI.JK BRIS.JK
## 0 853 0 0 0 0 0 0 0 0
## BRMS.JK BRPT.JK CTRA.JK ESSA.JK FILM.JK INCO.JK INKP.JK INTP.JK MAPA.JK MAPI.JK
## 0 0 0 0 0 0 0 0 0 0
## MBMA.JK MDKA.JK MNCN.JK PANI.JK PNLF.JK PWON.JK SMGR.JK SRTG.JK TPIA.JK SMRA.JK
## 807 0 0 0 0 0 0 0 0 0returns_data <- subset(returns_data, select = -c(`AMMN.JK`, `MBMA.JK`))
# Hitung rata-rata return harian tiap saham
average_returns <- sort(colMeans(returns_data, na.rm = TRUE), decreasing = TRUE)
cat("Berikut merupakan hasil dari rata-rata return yang telah diurutkan:\n")## Berikut merupakan hasil dari rata-rata return yang telah diurutkan:for (i in seq_along(average_returns)) {
cat(names(average_returns)[i], ": ", round(average_returns[i], 8), "\n", sep = "")
}## PANI.JK: 0.00780177
## FILM.JK: 0.00370862
## ARTO.JK: 0.00255227
## BRIS.JK: 0.00247745
## BRMS.JK: 0.00238954
## ESSA.JK: 0.00174274
## TPIA.JK: 0.00133382
## SRTG.JK: 0.00132649
## MAPA.JK: 0.00112875
## ANTM.JK: 0.00097623
## MDKA.JK: 0.00080665
## PNLF.JK: 0.00076175
## MAPI.JK: 0.00066184
## BMRI.JK: 0.00055212
## BBCA.JK: 0.0004349
## INCO.JK: 0.00040636
## CTRA.JK: 0.00033993
## INKP.JK: 0.00033453
## BBNI.JK: 0.00033437
## BRPT.JK: 0.00030867
## BBRI.JK: 0.00023763
## PWON.JK: 1.3e-06
## BBTN.JK: -0.00010261
## SMRA.JK: -0.00013221
## ACES.JK: -0.00018224
## INTP.JK: -0.00045253
## SMGR.JK: -0.0007555
## MNCN.JK: -0.00113481#menghapus saham yang memiliki rata-rata return negatif
positive_cols <- names(average_returns[average_returns >= 0])
returns_data_clean <- returns_data[, positive_cols]
#buat cek aja
average_returns <- sort(colMeans(returns_data_clean, na.rm = TRUE), decreasing = TRUE)
cat("\n Berikut merupakan hasil dari rata-rata return positif:\n")##
## Berikut merupakan hasil dari rata-rata return positif:for (i in seq_along(average_returns)) {
cat(names(average_returns)[i], ": ", round(average_returns[i], 8), "\n", sep = "")
}## PANI.JK: 0.00780177
## FILM.JK: 0.00370862
## ARTO.JK: 0.00255227
## BRIS.JK: 0.00247745
## BRMS.JK: 0.00238954
## ESSA.JK: 0.00174274
## TPIA.JK: 0.00133382
## SRTG.JK: 0.00132649
## MAPA.JK: 0.00112875
## ANTM.JK: 0.00097623
## MDKA.JK: 0.00080665
## PNLF.JK: 0.00076175
## MAPI.JK: 0.00066184
## BMRI.JK: 0.00055212
## BBCA.JK: 0.0004349
## INCO.JK: 0.00040636
## CTRA.JK: 0.00033993
## INKP.JK: 0.00033453
## BBNI.JK: 0.00033437
## BRPT.JK: 0.00030867
## BBRI.JK: 0.00023763
## PWON.JK: 1.3e-06melakukan seleksi saham menggunakan uji korelasi pearson.
# Hitung korelasi Pearson antar saham
cor_matrix <- cor(returns_data_clean, method = "pearson")
# Filter pasangan saham dengan korelasi < 0.002 dan return keduanya positif
saham_pairs <- combn(colnames(cor_matrix), 2, simplify = FALSE)
low_cor_positive_return <- list()
for (pair in saham_pairs) {
saham1 <- pair[1]
saham2 <- pair[2]
cor_val <- cor_matrix[saham1, saham2]
if (cor_val < 0.002 && average_returns[saham1] > 0 && average_returns[saham2] > 0) {
low_cor_positive_return[[length(low_cor_positive_return) + 1]] <- list(
Saham1 = saham1,
Saham2 = saham2,
Korelasi = round(cor_val, 3),
Return1 = round(average_returns[saham1], 6),
Return2 = round(average_returns[saham2], 6)
)
}
}
# Konversi ke data.frame
low_cor_df <- do.call(rbind, lapply(low_cor_positive_return, as.data.frame))
low_cor_df## Saham1 Saham2 Korelasi Return1 Return2
## PANI.JK PANI.JK TPIA.JK -0.014 0.007802 0.001334
## FILM.JK FILM.JK ARTO.JK -0.001 0.003709 0.002552
## BRMS.JK BRMS.JK MAPA.JK -0.016 0.002390 0.001129# Tampilkan jumlah pasangan dan tabel hasil
cat("📊 Jumlah pasangan saham dengan korelasi < 0.002 dan return positif:", nrow(low_cor_df), "\n\n")## 📊 Jumlah pasangan saham dengan korelasi < 0.002 dan return positif: 3print(low_cor_df, row.names = FALSE)## Saham1 Saham2 Korelasi Return1 Return2
## PANI.JK TPIA.JK -0.014 0.007802 0.001334
## FILM.JK ARTO.JK -0.001 0.003709 0.002552
## BRMS.JK MAPA.JK -0.016 0.002390 0.001129filter data sesuai yang ada pada output di atas
# Daftar nama saham yang ingin diambil
stocks <- c("PANI.JK", "TPIA.JK", "FILM.JK", "BRMS.JK", "ARTO.JK", "MAPA.JK")
# Saring hanya kolom yang sesuai dari returns_data_clean
df<- returns_data_clean[, stocks]
head(df,10)## PANI.JK TPIA.JK FILM.JK BRMS.JK ARTO.JK
## 2020-01-02 0.00000000 0.000000000 0.000000000 0.00000000 0.000000000
## 2020-01-03 0.00000000 0.000000000 0.005586592 0.03921571 0.000000000
## 2020-01-06 -0.12844036 -0.051094891 0.011111111 -0.01886794 -0.080645165
## 2020-01-07 0.03157884 -0.012820513 -0.010989011 0.00000000 0.000000000
## 2020-01-08 0.13265310 -0.033766234 0.200000000 -0.01923078 0.000000000
## 2020-01-09 0.00000000 0.018817204 -0.046296296 0.00000000 0.007017603
## 2020-01-10 0.00000000 -0.005277045 0.048543689 0.01960785 -0.006968699
## 2020-01-13 -0.06306303 0.000000000 0.250000000 0.01923078 0.007017603
## 2020-01-14 0.00000000 0.002652520 0.066666667 0.00000000 0.003484307
## 2020-01-15 0.00000000 0.000000000 -0.083333333 -0.01886794 0.003472209
## MAPA.JK
## 2020-01-02 0.000000000
## 2020-01-03 -0.029702970
## 2020-01-06 0.000000000
## 2020-01-07 0.000000000
## 2020-01-08 -0.002040816
## 2020-01-09 0.000000000
## 2020-01-10 0.000000000
## 2020-01-13 0.000000000
## 2020-01-14 0.000000000
## 2020-01-15 -0.018404908STATISTIKA DESKRIPTIF
# Matriks kovarian return
cov_matrix <- var(df, na.rm = TRUE)
cov_matrix## PANI.JK TPIA.JK FILM.JK BRMS.JK ARTO.JK
## PANI.JK 3.161754e-03 -2.458334e-05 3.006478e-04 4.513138e-05 1.179506e-04
## TPIA.JK -2.458334e-05 9.100323e-04 5.388248e-05 5.281716e-05 4.239436e-05
## FILM.JK 3.006478e-04 5.388248e-05 2.462521e-03 1.524592e-04 -1.958028e-06
## BRMS.JK 4.513138e-05 5.281716e-05 1.524592e-04 1.459029e-03 4.102494e-05
## ARTO.JK 1.179506e-04 4.239436e-05 -1.958028e-06 4.102494e-05 2.187548e-03
## MAPA.JK 6.353722e-05 2.993349e-05 9.978211e-05 -1.942866e-05 6.491745e-05
## MAPA.JK
## PANI.JK 6.353722e-05
## TPIA.JK 2.993349e-05
## FILM.JK 9.978211e-05
## BRMS.JK -1.942866e-05
## ARTO.JK 6.491745e-05
## MAPA.JK 1.039668e-03# Hitung statistik deskriptif
stats_summary <- data.frame(
Mean = colMeans(df, na.rm = TRUE),
Stdev = apply(df, 2, sd, na.rm = TRUE),
Variance = apply(df, 2, var, na.rm = TRUE),
Skewness = apply(df, 2, skewness, na.rm = TRUE),
Kurtosis = apply(df, 2, kurtosis, na.rm = TRUE)
)
# Tampilkan statistik
print(round(stats_summary, 6))## Mean Stdev Variance Skewness Kurtosis
## PANI.JK 0.007802 0.056229 0.003162 1.833731 6.154939
## TPIA.JK 0.001334 0.030167 0.000910 0.961615 14.778444
## FILM.JK 0.003709 0.049624 0.002463 1.658249 5.709237
## BRMS.JK 0.002390 0.038197 0.001459 2.113047 11.018433
## ARTO.JK 0.002552 0.046771 0.002188 1.993495 8.359635
## MAPA.JK 0.001129 0.032244 0.001040 1.108037 5.990681#visualisasi data plot time series
FILM <- ts(df$FILM.JK, freq=12)
TPIA <- ts(df$TPIA.JK, freq= 12)
PANI <- ts(df$PANI.JK, freq= 12)
MAPA <- ts(df$MAPA.JK, freq= 12)
ARTO <- ts(df$ARTO.JK, freq= 12)
BRMS <- ts(df$BRMS.JK, freq= 12)
plot(cbind(FILM, TPIA, PANI, MAPA, ARTO, BRMS))
# Normalitas
library(nortest)
# Pastikan df jadi data.frame
df_df <- as.data.frame(df)
for (col in colnames(df_df)) {
if (is.numeric(df_df[[col]])) {
pval <- ad.test(df_df[[col]])$p.value
if (pval < 0.05) {
cat(col, ": Tidak memenuhi asumsi normalitas (p-value =", pval, ")\n")
} else {
cat(col, ": Memenuhi asumsi normalitas (p-value =", pval, ")\n")
}
}
}## PANI.JK : Tidak memenuhi asumsi normalitas (p-value = 3.7e-24 )
## TPIA.JK : Tidak memenuhi asumsi normalitas (p-value = 3.7e-24 )
## FILM.JK : Tidak memenuhi asumsi normalitas (p-value = 3.7e-24 )
## BRMS.JK : Tidak memenuhi asumsi normalitas (p-value = 3.7e-24 )
## ARTO.JK : Tidak memenuhi asumsi normalitas (p-value = 3.7e-24 )
## MAPA.JK : Tidak memenuhi asumsi normalitas (p-value = 3.7e-24 )library(MVN)## Warning: package 'MVN' was built under R version 4.3.3mvn_result1 <- mvn(data = df_df, mvnTest = "mardia")
mvn_result2 <- mvn(data = df_df, mvnTest = "hz")
mvn_result3 <- mvn(data = df_df, mvnTest = "royston")
mvn_result1$multivariateNormality## Test Statistic p value Result
## 1 Mardia Skewness 3652.78759872382 0 NO
## 2 Mardia Kurtosis 93.6399270581825 0 NO
## 3 MVN <NA> <NA> NOmvn_result2$multivariateNormality## Test HZ p value MVN
## 1 Henze-Zirkler 9.52231 0 NOmvn_result3$multivariateNormality## Test H p value MVN
## 1 Royston 802.3277 4.837073e-170 NOMVEP
# Hitung matriks kovarians return harian
cov_matrix <- cov(df)
# Hitung banyak aset
n_assets <- ncol(df)
# Buat vektor satu
one_vec <- rep(1, n_assets)
# Hitung bobot MVEP (Minimum Variance)
weights_mvep <- solve(cov_matrix) %*% one_vec / as.numeric(t(one_vec) %*% solve(cov_matrix) %*% one_vec)
# Hitung return harian rata-rata dari tiap aset
mean_returns <- colMeans(df)
# Hitung return portofolio harian
port_return_mvep <- as.numeric(t(weights_mvep) %*% mean_returns)
# Hitung variansi portofolio harian
port_var_mvep <- as.numeric(t(weights_mvep) %*% cov_matrix %*% weights_mvep)
# Hitung standar deviasi (volatilitas harian)
port_sd_mvep <- sqrt(port_var_mvep)
# Asumsikan risk-free rate harian, misalnya 4% tahunan = 0.04 / 252
rf <- 0.05 / 252
# Hitung Sharpe Ratio
sharpe_mvep <- (port_return_mvep - rf) / port_sd_mvep
daily_port_returns <- df %*% weights_mvep
skew_p <- skewness(daily_port_returns)
kurt_p <- kurtosis(daily_port_returns)
z <- qnorm(0.95) # z-score normal
# Cornish-Fisher adjustment
z_cf <- z + (1/6)*(z^2 - 1)*skew_p + (1/24)*(z^3 - 3*z)*kurt_p - (1/36)*(2*z^3 - 5*z)*skew_p^2
# VaR Cornish-Fisher
VaR_cf_mvep <- -(port_return_mvep + z_cf * port_sd_mvep)
cat("Return Portofolio MVEP:", round(port_return_mvep, 6), "\n",
"Risiko Portofolio MVEP:", round(port_sd_mvep, 6), "\n",
"Sharpe Ratio MVEP:", round(sharpe_mvep, 4), "\n",
"Cornish-Fisher VaR (95%) MVEP:", round(VaR_cf_mvep, 6), "\n")## Return Portofolio MVEP: 0.002288
## Risiko Portofolio MVEP: 0.0171
## Sharpe Ratio MVEP: 0.1222
## Cornish-Fisher VaR (95%) MVEP: -0.031334Modelling SIM (Single Index Model)
pada analisis ini menggunakan stock market dari IHSG karena dari indeks ECONOMIC30 tidak dapat diupload oleh device saya.
# Return harian IHSG
market_return <- dailyReturn(Cl(JKSE))
market_return <- market_return[1:1211, ]
head(market_return,5)## daily.returns
## 2020-01-02 0.000000000
## 2020-01-03 0.006347458
## 2020-01-06 -0.010447275
## 2020-01-07 0.003506784
## 2020-01-08 -0.008545501# Gabungkan semua return
df_sim <- data.frame(df, market_return)
colnames(df_sim)[colnames(df_sim) == "daily.returns"] <- "IHSG"
head(df_sim,5)## PANI.JK TPIA.JK FILM.JK BRMS.JK ARTO.JK
## 2020-01-02 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
## 2020-01-03 0.00000000 0.00000000 0.005586592 0.03921571 0.00000000
## 2020-01-06 -0.12844036 -0.05109489 0.011111111 -0.01886794 -0.08064517
## 2020-01-07 0.03157884 -0.01282051 -0.010989011 0.00000000 0.00000000
## 2020-01-08 0.13265310 -0.03376623 0.200000000 -0.01923078 0.00000000
## MAPA.JK IHSG
## 2020-01-02 0.000000000 0.000000000
## 2020-01-03 -0.029702970 0.006347458
## 2020-01-06 0.000000000 -0.010447275
## 2020-01-07 0.000000000 0.003506784
## 2020-01-08 -0.002040816 -0.008545501# Estimasi Alpha, Beta, dan Variansi Residual via regresi linier
alphas <- c()
betas <- c()
var_res <- c()
for (stock in stocks) {
model <- lm(df_sim[[stock]] ~ df_sim$IHSG)
alphas <- c(alphas, coef(model)[1])
betas <- c(betas, coef(model)[2])
var_res <- c(var_res, var(residuals(model)))
}
parameters_df <- data.frame(Code = stocks, Beta = betas, Alpha = alphas, Residual_Var = var_res)
parameters_df## Code Beta Alpha Residual_Var
## 1 PANI.JK 0.5008211 0.007725795 0.0031350618
## 2 TPIA.JK 0.8495241 0.001204952 0.0008332312
## 3 FILM.JK 0.9157358 0.003569706 0.0023732817
## 4 BRMS.JK 0.6023577 0.002298161 0.0014204167
## 5 ARTO.JK 0.8395056 0.002424924 0.0021125479
## 6 MAPA.JK 0.7877972 0.001009247 0.0009736219# Hitung ERB (Expected Return Beta-adjusted)
Rf <- 0.05 / 252 # Risk-free rate harian
parameters_df <- parameters_df %>%
mutate(
ERB = (Alpha + Beta * mean(df_sim$IHSG, na.rm = TRUE) - Rf) / Beta
) %>%
arrange(desc(ERB))
parameters_df## Code Beta Alpha Residual_Var ERB
## 1 PANI.JK 0.5008211 0.007725795 0.0031350618 0.015181776
## 2 FILM.JK 0.9157358 0.003569706 0.0023732817 0.003833207
## 3 BRMS.JK 0.6023577 0.002298161 0.0014204167 0.003637576
## 4 ARTO.JK 0.8395056 0.002424924 0.0021125479 0.002803864
## 5 TPIA.JK 0.8495241 0.001204952 0.0008332312 0.001336522
## 6 MAPA.JK 0.7877972 0.001009247 0.0009736219 0.001180936# Hitung Ai, Bi, Ci
sigma_m2 <- var(df_sim$IHSG, na.rm = TRUE)
parameters_df <- parameters_df %>%
mutate(
A_i = (Alpha + Beta * mean(df_sim$IHSG, na.rm = TRUE) - ERB) / Residual_Var,
B_i = Beta^2 / Residual_Var,
C_i = map_dbl(row_number(), ~ (sigma_m2 * sum(A_i[1:.x])) / (1 + sigma_m2 * sum(B_i[1:.x])))
)
parameters_df## Code Beta Alpha Residual_Var ERB A_i B_i
## 1 PANI.JK 0.5008211 0.007725795 0.0031350618 0.015181776 -2.35402347 80.00537
## 2 FILM.JK 0.9157358 0.003569706 0.0023732817 0.003833207 -0.05249663 353.33864
## 3 BRMS.JK 0.6023577 0.002298161 0.0014204167 0.003637576 -0.87864462 255.44254
## 4 ARTO.JK 0.8395056 0.002424924 0.0021125479 0.002803864 -0.11909396 333.61121
## 5 TPIA.JK 0.8495241 0.001204952 0.0008332312 0.001336522 -0.00324238 866.13548
## 6 MAPA.JK 0.7877972 0.001009247 0.0009736219 0.001180936 -0.05359911 637.43885
## C_i
## 1 -0.0002483959
## 2 -0.0002448079
## 3 -0.0003257256
## 4 -0.0003267264
## 5 -0.0003019381
## 6 -0.0002902909# Tentukan C* (C_star)
C_star <- tail(parameters_df$C_i[parameters_df$ERB > parameters_df$C_i], 1)
# Seleksi saham berdasarkan ERB >= C_star
selected_df <- parameters_df %>%
filter(ERB >= C_star) %>%
mutate(
Z_i = (Beta / Residual_Var) * (ERB - C_star),
w_i = Z_i / sum(Z_i)
)
selected_df## Code Beta Alpha Residual_Var ERB A_i B_i
## 1 PANI.JK 0.5008211 0.007725795 0.0031350618 0.015181776 -2.35402347 80.00537
## 2 FILM.JK 0.9157358 0.003569706 0.0023732817 0.003833207 -0.05249663 353.33864
## 3 BRMS.JK 0.6023577 0.002298161 0.0014204167 0.003637576 -0.87864462 255.44254
## 4 ARTO.JK 0.8395056 0.002424924 0.0021125479 0.002803864 -0.11909396 333.61121
## 5 TPIA.JK 0.8495241 0.001204952 0.0008332312 0.001336522 -0.00324238 866.13548
## 6 MAPA.JK 0.7877972 0.001009247 0.0009736219 0.001180936 -0.05359911 637.43885
## C_i Z_i w_i
## 1 -0.0002483959 2.471638 0.2520271
## 2 -0.0002448079 1.591061 0.1622367
## 3 -0.0003257256 1.665695 0.1698470
## 4 -0.0003267264 1.229586 0.1253780
## 5 -0.0003019381 1.658623 0.1691259
## 6 -0.0002902909 1.190430 0.1213853# Hitung alpha_p, beta_p, return, dan varian portofolio
alpha_p <- sum(selected_df$w_i * selected_df$Alpha)
beta_p <- sum(selected_df$w_i * selected_df$Beta)
portfolio_return <- alpha_p + beta_p * mean(df_sim$IHSG, na.rm = TRUE)
portfolio_variance <- sum(selected_df$w_i^2 * selected_df$Residual_Var) + (beta_p^2 * sigma_m2)
portfolio_summary <- data.frame(
Alpha = alpha_p,
Beta = beta_p,
Return = portfolio_return,
Variance = portfolio_variance
)
print(portfolio_summary)## Alpha Beta Return Variance
## 1 0.003546911 0.7216542 0.003656382 0.0004293831# Optimisasi kuadrat (Quadratic Optimization)
Risk_matrix <- selected_df$Beta %*% t(selected_df$Beta) * sigma_m2 + diag(selected_df$Residual_Var)
N <- nrow(Risk_matrix)
Dmat <- Risk_matrix
dvec <- rep(0, N)
Amat <- cbind(rep(1, N), diag(N))
bvec <- c(1, rep(0, N))
optimal_weights <- solve.QP(Dmat, dvec, Amat, bvec, meq = 1)$solution
alpha_p2 <- sum(optimal_weights * selected_df$Alpha)
beta_p2 <- sum(optimal_weights * selected_df$Beta)
portfolio_return2 <- alpha_p2 + beta_p2 * mean(df_sim$IHSG, na.rm = TRUE)
portfolio_variance2 <- sum(optimal_weights^2 * selected_df$Residual_Var) + (beta_p2^2 * sigma_m2)
opt_summary <- data.frame(
Alpha = alpha_p2,
Beta = beta_p2,
Return = portfolio_return2,
Variance = portfolio_variance2
)
print(opt_summary)## Alpha Beta Return Variance
## 1 0.002265342 0.7659425 0.002381531 0.0003041017# Untuk seluruh saham (tanpa seleksi ERB)
parameters_df <- parameters_df %>%
mutate(
Z_i = (Beta / Residual_Var) * (ERB - C_star),
w_i = Z_i / sum(Z_i)
)
Risk_matrix_all <- parameters_df$Beta %*% t(parameters_df$Beta) * sigma_m2 + diag(parameters_df$Residual_Var)
N_all <- nrow(Risk_matrix_all)
Dmat_all <- Risk_matrix_all
dvec_all <- rep(0, N_all)
Amat_all <- cbind(rep(1, N_all), diag(N_all))
bvec_all <- c(1, rep(0, N_all))
optimal_weights_all <- solve.QP(Dmat_all, dvec_all, Amat_all, bvec_all, meq = 1)$solution
alpha_p3 <- sum(optimal_weights_all * parameters_df$Alpha)
beta_p3 <- sum(optimal_weights_all * parameters_df$Beta)
portfolio_return3 <- alpha_p3 + beta_p3 * mean(df_sim$IHSG, na.rm = TRUE)
portfolio_variance3 <- sum(optimal_weights_all^2 * parameters_df$Residual_Var) + (beta_p3^2 * sigma_m2)
opt2_summary <- data.frame(
Alpha = alpha_p3,
Beta = beta_p3,
Return = portfolio_return3,
Variance = portfolio_variance3
)
print(opt2_summary)## Alpha Beta Return Variance
## 1 0.002265342 0.7659425 0.002381531 0.0003041017# Output
# Copy data frame
parameters_df_rounded <- parameters_df
# Hanya kolom numerik yang dibulatkan
num_cols <- sapply(parameters_df, is.numeric)
parameters_df_rounded[num_cols] <- round(parameters_df[num_cols], 5)
print(parameters_df_rounded)## Code Beta Alpha Residual_Var ERB A_i B_i C_i
## 1 PANI.JK 0.50082 0.00773 0.00314 0.01518 -2.35402 80.00537 -0.00025
## 2 FILM.JK 0.91574 0.00357 0.00237 0.00383 -0.05250 353.33864 -0.00024
## 3 BRMS.JK 0.60236 0.00230 0.00142 0.00364 -0.87864 255.44254 -0.00033
## 4 ARTO.JK 0.83951 0.00242 0.00211 0.00280 -0.11909 333.61121 -0.00033
## 5 TPIA.JK 0.84952 0.00120 0.00083 0.00134 -0.00324 866.13548 -0.00030
## 6 MAPA.JK 0.78780 0.00101 0.00097 0.00118 -0.05360 637.43885 -0.00029
## Z_i w_i
## 1 2.47164 0.25203
## 2 1.59106 0.16224
## 3 1.66570 0.16985
## 4 1.22959 0.12538
## 5 1.65862 0.16913
## 6 1.19043 0.12139cat("\nC* (C_star):", round(C_star, 6), "\n\n")##
## C* (C_star): -0.00029cat("Portofolio (seleksi ERB >= C*):\n")## Portofolio (seleksi ERB >= C*):cat("Return :", round(portfolio_return, 6), "\n")## Return : 0.003656cat("Variance:", round(portfolio_variance, 6), "\n\n")## Variance: 0.000429cat("Portofolio Optimasi (tanpa seleksi):\n")## Portofolio Optimasi (tanpa seleksi):cat("Return :", round(portfolio_return2, 6), "\n")## Return : 0.002382cat("Variance:", round(portfolio_variance2, 6), "\n\n")## Variance: 0.000304cat("Portofolio Optimasi (semua saham):\n")## Portofolio Optimasi (semua saham):cat("Return :", round(portfolio_return3, 6), "\n")## Return : 0.002382cat("Variance:", round(portfolio_variance3, 6), "\n")## Variance: 0.000304returns <- do.call(merge, lapply(stocks, function(x) dailyReturn(Cl(get(x)))))
# Bobot saham sesuai urutan kolom pada 'returns'
#weights <- c(0.08397910, 0.09668465, 0.17952719, 0.11155887, 0.28186275, 0.24638744)
# Pastikan jumlah bobot = 1
#sum(weights) # seharusnya hasilnya 1
# Hitung return portofolio harian
portfolio_returns <- Return.portfolio(df, weights = optimal_weights_all, geometric = FALSE)
portfolio_risk <- sqrt(portfolio_variance)
skew_p <- skewness(portfolio_returns)
kurt_p <- kurtosis(portfolio_returns)
z <- qnorm(0.95) # z-score normal
# Cornish-Fisher adjustment
z_cf <- z + (1/6)*(z^2 - 1)*skew_p + (1/24)*(z^3 - 3*z)*kurt_p - (1/36)*(2*z^3 - 5*z)*skew_p^2
# VaR Cornish-Fisher
VaR_sim <- -(portfolio_return + z_cf * portfolio_risk)
#sharpe ratio
rf <- 0.05 / 252 # 3% per tahun dikonversi ke harian
sharpe_ratio_sim <- (portfolio_return - rf) / portfolio_risk
cat("Return Portofolio SIM:", round(portfolio_return, 6), "\n",
"Risiko Portofolio SIM:", round(portfolio_risk, 6), "\n",
"Sharpe Ratio SIM:", round(sharpe_ratio_sim, 4), "\n",
"Cornish-Fisher VaR (95%) SIM:", round(VaR_sim, 6), "\n")## Return Portofolio SIM: 0.003656
## Risiko Portofolio SIM: 0.020722
## Sharpe Ratio SIM: 0.1669
## Cornish-Fisher VaR (95%) SIM: -0.042687CAPM
# Gabungkan semua return
df_capm <- data.frame(df, market_return)
colnames(df_capm)[colnames(df_capm) == "daily.returns"] <- "IHSG"
head(df_capm,5)## PANI.JK TPIA.JK FILM.JK BRMS.JK ARTO.JK
## 2020-01-02 0.00000000 0.00000000 0.000000000 0.00000000 0.00000000
## 2020-01-03 0.00000000 0.00000000 0.005586592 0.03921571 0.00000000
## 2020-01-06 -0.12844036 -0.05109489 0.011111111 -0.01886794 -0.08064517
## 2020-01-07 0.03157884 -0.01282051 -0.010989011 0.00000000 0.00000000
## 2020-01-08 0.13265310 -0.03376623 0.200000000 -0.01923078 0.00000000
## MAPA.JK IHSG
## 2020-01-02 0.000000000 0.000000000
## 2020-01-03 -0.029702970 0.006347458
## 2020-01-06 0.000000000 -0.010447275
## 2020-01-07 0.000000000 0.003506784
## 2020-01-08 -0.002040816 -0.008545501# Hitung beta dan expected return via CAPM
rf <- 0.05/252
betas <- c()
er <- c()
for (stock in stocks) {
model <- lm(df_capm[[stock]] ~ df_capm$IHSG)
beta <- coef(model)[2]
betas <- c(betas, beta)
er_i <- rf + beta * (mean(df_capm$IHSG, na.rm = TRUE) - rf)
er <- c(er, er_i)
}# Buat dataframe parameter
parameters_df <- data.frame(Code = stocks, Beta = betas, Expected_Return = er)
print(parameters_df)## Code Beta Expected_Return
## 1 PANI.JK 0.5008211 0.0001750151
## 2 TPIA.JK 0.8495241 0.0001587242
## 3 FILM.JK 0.9157358 0.0001556309
## 4 BRMS.JK 0.6023577 0.0001702714
## 5 ARTO.JK 0.8395056 0.0001591923
## 6 MAPA.JK 0.7877972 0.0001616080# Matriks kovarian return
cov_matrix <- var(df, na.rm = TRUE)
cov_matrix## PANI.JK TPIA.JK FILM.JK BRMS.JK ARTO.JK
## PANI.JK 3.161754e-03 -2.458334e-05 3.006478e-04 4.513138e-05 1.179506e-04
## TPIA.JK -2.458334e-05 9.100323e-04 5.388248e-05 5.281716e-05 4.239436e-05
## FILM.JK 3.006478e-04 5.388248e-05 2.462521e-03 1.524592e-04 -1.958028e-06
## BRMS.JK 4.513138e-05 5.281716e-05 1.524592e-04 1.459029e-03 4.102494e-05
## ARTO.JK 1.179506e-04 4.239436e-05 -1.958028e-06 4.102494e-05 2.187548e-03
## MAPA.JK 6.353722e-05 2.993349e-05 9.978211e-05 -1.942866e-05 6.491745e-05
## MAPA.JK
## PANI.JK 6.353722e-05
## TPIA.JK 2.993349e-05
## FILM.JK 9.978211e-05
## BRMS.JK -1.942866e-05
## ARTO.JK 6.491745e-05
## MAPA.JK 1.039668e-03# Return berlebih
excess_returns <- parameters_df$Expected_Return - rf
excess_returns## [1] -2.339761e-05 -3.968849e-05 -4.278181e-05 -2.814125e-05 -3.922044e-05
## [6] -3.680470e-05# Hitung bobot portofolio optimal manual Sharpe
weights <- solve(cov_matrix) %*% excess_returns /
as.numeric(t(rep(1, length(stocks))) %*% solve(cov_matrix) %*% excess_returns)
weights <- as.numeric(weights) # ubah ke vektor biasa
cat("\nBobot portofolio optimal (manual Sharpe):\n")##
## Bobot portofolio optimal (manual Sharpe):print(round(weights, 5))## [1] 0.04047 0.32798 0.11090 0.13246 0.12755 0.26064# Hitung bobot optimal dengan optimasi kuadrat
N <- nrow(cov_matrix)
Dmat <- cov_matrix
dvec <- rep(0, N)
Amat <- cbind(excess_returns)
bvec <- c(1)
sol <- solve.QP(Dmat, dvec, Amat, bvec, meq = 1)$solution
optimal_weights <- sol / sum(sol)
cat("\nBobot portofolio optimal (QP):\n")##
## Bobot portofolio optimal (QP):print(round(optimal_weights, 5))## [1] 0.04047 0.32798 0.11090 0.13246 0.12755 0.26064# Hitung bobot optimal menggunakan NMOF (maxSharpe)
nmof_weights <- maxSharpe(excess_returns, cov_matrix)
cat("\nBobot portofolio optimal (NMOF maxSharpe):\n")##
## Bobot portofolio optimal (NMOF maxSharpe):print(round(nmof_weights, 5))## [1] 0.04047 0.32798 0.11090 0.13246 0.12755 0.26064# Hitung expected return dan risiko portofolio (manual Sharpe)
port_return <- sum(nmof_weights * parameters_df$Expected_Return)
port_risk <- sqrt(as.numeric(t(nmof_weights) %*% cov_matrix %*% weights))
cat("\nKinerja Portofolio:\n")##
## Kinerja Portofolio:cat("Expected Return:", round(port_return, 6), "\n")## Expected Return: 0.000161cat("Risk (Std Dev):", round(port_risk, 6), "\n")## Risk (Std Dev): 0.017378# Hitung return portofolio harian dengan bobot manual Sharpe
portfolio_returns <- Return.portfolio(returns, weights = nmof_weights, geometric = FALSE)
skew_p <- skewness(portfolio_returns)
kurt_p <- kurtosis(portfolio_returns)
z <- qnorm(0.95) # z-score normal
# Cornish-Fisher adjustment
z_cf <- z + (1/6)*(z^2 - 1)*skew_p + (1/24)*(z^3 - 3*z)*kurt_p - (1/36)*(2*z^3 - 5*z)*skew_p^2
# VaR Cornish-Fisher
VaR_capm <- -(port_return + z_cf * port_risk)
#sharpe ratio
rf <- 0.05 / 252 # 3% per tahun dikonversi ke harian
sharpe_ratio_capm <- (port_return - rf) / port_risk
cat("Return Portofolio CAPM:", round(port_return, 6), "\n",
"Risiko Portofolio CAPM:", round(port_risk, 6), "\n",
"Sharpe Ratio CAPM:", round(sharpe_ratio_capm, 6), "\n",
"Cornish-Fisher VaR (95%) CAPM:", round(VaR_capm, 6), "\n")## Return Portofolio CAPM: 0.000161
## Risiko Portofolio CAPM: 0.017378
## Sharpe Ratio CAPM: -0.002131
## Cornish-Fisher VaR (95%) CAPM: -0.032171NSGA-II
# Fungsi objektif untuk NSGA-II
obj_func <- function(weights) {
weights <- weights / sum(weights) # Normalisasi agar total bobot = 1
port_return <- mean(Return.portfolio(returns, weights = weights))
port_risk <- sd(Return.portfolio(returns, weights = weights))
return(c(-port_return, port_risk)) # NSGA-II meminimalkan semuanya
}# NSGA-II parameters
set.seed(1234)
result <- nsga2(
fn = obj_func,
idim = length(stocks),
odim = 2,
lower.bounds = rep(0, length(stocks)),
upper.bounds = rep(1, length(stocks)),
popsize = 16,
generations = 16
)
result## $par
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.913096993 0.17046005 0.0006203759 0.3026934 0.14793105 0.004660528
## [2,] 0.001990187 0.89723928 0.0336138233 0.2286067 0.02133430 0.639518634
## [3,] 0.018114141 0.84994784 0.3127359855 0.6405013 0.10842770 0.058287617
## [4,] 0.162173826 0.04620836 0.9984649836 0.5890102 0.10687156 0.058733070
## [5,] 0.014627256 0.86396626 0.0636536690 0.4150563 0.05914759 0.864288051
## [6,] 0.045050401 0.35244512 0.3203001125 0.4525959 0.12774637 0.674377886
## [7,] 0.186722790 0.36578899 0.3166124548 0.3109750 0.13552313 0.221765019
## [8,] 0.826196069 0.06636104 0.1838698045 0.2401247 0.24578677 0.715342068
## [9,] 0.913096993 0.17046005 0.0006203759 0.3632809 0.22494186 0.004660528
## [10,] 0.001990187 0.89723928 0.1934516035 0.2281267 0.01585616 0.639518634
## [11,] 0.085569783 0.36798781 0.3329479013 0.2078471 0.10969413 0.625279768
## [12,] 0.121955291 0.17046005 0.0010405664 0.2932181 0.14793105 0.004660528
## [13,] 0.002205558 0.46062913 0.0571305258 0.2118783 0.09744628 0.639518634
## [14,] 0.062115167 0.89723928 0.1101810128 0.2286480 0.02162527 0.639518634
## [15,] 0.909099675 0.06570924 0.4249515935 0.4863109 0.24578677 0.715342068
## [16,] 0.909099675 0.81814864 0.4249515935 0.5824875 0.24382658 0.715342068
##
## $value
## [,1] [,2]
## [1,] -0.006751989 0.04240386
## [2,] -0.001633343 0.01839800
## [3,] -0.002987696 0.02388816
## [4,] -0.004816659 0.03198870
## [5,] -0.002606565 0.02237437
## [6,] -0.003654701 0.02680895
## [7,] -0.005143506 0.03353498
## [8,] -0.006249882 0.03988557
## [9,] -0.006645481 0.04147161
## [10,] -0.001852632 0.02139810
## [11,] -0.004326273 0.03014160
## [12,] -0.005472334 0.03617954
## [13,] -0.001883278 0.02194279
## [14,] -0.003937992 0.02958950
## [15,] -0.006103846 0.03872961
## [16,] -0.005853941 0.03776930
##
## $pareto.optimal
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE
##
## attr(,"class")
## [1] "nsga2" "mco"# Ambil solusi dengan return tertinggi
best_idx <- which.max(-result$value[,1])
best_weights <- result$par[best_idx, ] / sum(result$par[best_idx, ])
names(best_weights) <- stocks
print(round(best_weights, 4))## PANI.JK TPIA.JK FILM.JK BRMS.JK ARTO.JK MAPA.JK
## 0.5931 0.1107 0.0004 0.1966 0.0961 0.0030# Plot return vs risk (front pareto)
plot(-result$value[,1], result$value[,2],
xlab = "Expected Return", ylab = "Risk (Std Dev)",
main = "Efficient Frontier (NSGA-II)",
pch = 19, col = "blue")
# Hitung return portofolio
port_returns <- Return.portfolio(returns, weights = best_weights, geometric = FALSE)
# Statistik return
mu <- mean(port_returns)
sigma <- sd(port_returns)
skew <- skewness(port_returns)
kurt <- kurtosis(port_returns)
# Cornish-Fisher Z-score adjustment (95% confidence)
z <- qnorm(0.05)
z_cf <- z + (1/6)*(z^2 - 1)*skew + (1/24)*(z^3 - 3*z)*kurt - (1/36)*(2*z^3 - 5*z)*skew^2
# Cornish-Fisher VaR (1-day)
VaR_cf_nsga2 <- (mu + z_cf * sigma)
print(paste("Cornish-Fisher VaR (95%):", round(VaR_cf_nsga2, 5)))## [1] "Cornish-Fisher VaR (95%): -0.03962"# Sharpe ratio
sharpe_nsga2 <- mean(port_returns) / sd(port_returns)
print(paste("Sharpe Ratio:", round(sharpe_nsga2, 4)))## [1] "Sharpe Ratio: 0.1592"