몬테카를로 시뮬레이션은 불확실한 사건의 가능한 결과를 추정하는데 사용되는 수학적 기법이다. 이 방법은 제2차 세계대전 중 존 폰 노이만과 스타니스와프 울람이 불확실한 상황에서 의사결정을 개선하기 위해 개발한 것으로 알려져 있다. 무작위 난수 추출(random sampling)을 통해 복잡하고 불확실한 미래의 특정 확률 분포를 만들 수 있다. 예를 들어, 특정 주식의 100일 후 가격 분포를 시뮬레이션하는 데 활용하는 식이다. 물론 모델링에 있어서 많은 한계점 존재하지만 여전히 많은 분야에서 광범위하게 사용된다는 점에서 본 기법에 대한 학습은 충분한 의의를 가진다.
library(purrr)
package <- c("tidyverse", "quantmod", "broom", "pander", "PerformanceAnalytics", "PortfolioAnalytics", "scales")
pack_map <- map(package, ~ library(.x, character.only = TRUE))
options(scipen = 999)
PF_ticker <- tibble(ticker = c("005930.KS", "005380.KS", "032830.KS", "034020.KS", "139480.KS"))
PF_name <- c("삼성전자", "현대차", "삼성생명", "두산에너빌리티", "이마트")
PF_df <- PF_ticker %>%
mutate(data = map(ticker, ~getSymbols(Symbols = .x,
auto.assign = FALSE,
warnings = FALSE) %>%
Ad() %>%
na.omit()),
daily_ret = map(data, ~dailyReturn(.x, type = "log")))
PF_daily_ret_matrix <- purrr::reduce(PF_df$daily_ret, cbind) %>%
na.omit() %>%
`colnames<-`(PF_name)
mu <- colMeans(PF_daily_ret_matrix)
sigma <- cov(PF_daily_ret_matrix)
n_assets <- length(PF_name)
n_sim <- 10000
n_days <- 250
initial_price <- rep(100, n_assets) # 모든 자산 시작가 100 고정
chol함수는 상삼각행렬 반환L <- chol(sigma)
Z_pos <- array(rnorm(n_assets * (n_sim / 2) * n_days), dim = c(n_assets, (n_sim / 2), n_days))
Z_neg <- - Z_pos
Z_full <- array(0, dim = c(n_assets, n_sim, n_days))
Z_full[, 1:(n_sim/2), ] <- Z_pos
Z_full[, (n_sim/2 + 1):n_sim, ] <- Z_neg
sim_returns <- array(0, dim = c(n_assets, n_sim, n_days))
for (t in 1:n_days) {
sim_returns[,,t] <- matrix(rep(mu - 0.5 * diag(sigma), n_sim), nrow = n_assets) + L %*% Z_full[,,t]
}
cum_rets <- apply(sim_returns, c(1,2), sum)
final_prices <- initial_price * exp(cum_rets)
# 동일 비중 포트폴리오
w <- rep(1/n_assets, n_assets)
PF_values <- colSums(cum_rets * w)
hist(PF_values, breaks = 50, probability = TRUE, col = "lightblue", xlab = "Portfoili Return", main = "250일 후 포트폴리오 누적 수익률 분포")
var_95 <- quantile(PF_values, probs = 0.05)
abline(v = var_95, col = "red", lwd = 2)
text(x = var_95 - 0.15, y = 1.6, labels = paste0("95% VaR", "\n", percent(var_95, accuracy = 0.01)))