Efficient Frontier (Mean Variance Optimization) using ggplot2 in R

library(quantmod)
library(PerformanceAnalytics)
library(magrittr)
library(ggplot2)

symbols = c('005930.KS', '068270.KS', '005380.KS', '055550.KS', '017670.KS')
names = c('삼성전자', '셀트리온', '현대차', '신한지주', 'SK텔레콤')
getSymbols(symbols, src='yahoo')

prices = do.call(cbind, lapply(symbols, function(x) Cl(get(x)))) %>% setNames(symbols)
rets = Return.calculate(prices) %>% na.omit()
names(rets) = names
covs = cov(rets) * 252

먼저 필요한 패키지들을 불러온 후, 삼성전자, 셀트리온, 현대차, 신한지주, SK텔레콤에 해당하는 종목 가격을 불러옵니다. 그 후 수익률과 공분산을 계산하도록 합니다. 공분산은 연율화를 해주기 위해 252를 곱해줍니다.

library(quantmod)
charts.PerformanceSummary(rets, main = 'rets')

수익률 그래프를 그려보면 셀트리온의 누적수익률이 압도적으로 높습니다.

yearly_mean = Return.annualized(rets, geometric = FALSE) %>% c()
yearly_vol = StdDev.annualized(rets) %>% c()
# diag(covs) ^ (1/2)

yearly_stat = cbind(yearly_mean, yearly_vol) %>%
  data.frame()

연율화 수익률과 연율화 변동성을 계산하도록 합니다. 연율화 변동성의 경우 위에서 구한 연율화 분산-공분산의 대각 부분의 제곱근을 구해도 같습니다. 이를 cbind()를 이용해 열로 묶습니다.

yearly_stat %>%
  ggplot(aes(x = yearly_vol, y = yearly_mean, label = names, color = names)) +
  geom_point(shape = 18) +
  scale_x_continuous(expand = c(0.02, 0.02)) + 
  scale_y_continuous(expand = c(0.02, 0.02)) +
  geom_text(size = 3, vjust = -0.7) +
  theme(legend.position = 'none') +
  xlab('Volatility') +
  ylab('Return')

ggplot을 이용해 x축에는 연율화 변동성, y축에는 연율화 수익률을 나타내 줍니다. 셀트리온의 수익률과 변동성이 압도적으로 높습니다.

# Simulation
numAssets = length(symbols)
numPortfolio = 10000
sim = list()

for (i in 1 : numPortfolio) {
  
  wt = runif(numAssets, min=0, max=1)
  wt = wt / sum(wt)
  expected_return = wt %*% yearly_mean
  expected_vol = (t(wt) %*% covs %*% wt) %>% sqrt()
  expected_sharpe = (expected_return-0.05) / expected_vol
  
  sim[[i]] = c(wt, expected_return, expected_vol, expected_sharpe)
  
}

sim = do.call(rbind, sim) %>% data.frame()
colnames(sim) = c(names, 'Return', 'Volatility', 'Sharpe')

본격적으로 시뮬레이션을 실행합니다. 시뮬레이션 횟수를 총 10,000번이며 많을 수록 좋습니다. 먼저 runif() 함수를 이용해 난수를 생성하며, 공매도 방지를 위해 0에서 1 사이의 값에서 난수를 생성하도록 합니다. 그 후 합이 1이 되도록 정규화를 해줍니다. (개별 종목별 비중제한을 두고 싶으면 while 구문을 이용해 정규화를 계속해주면 됩니다.)

기대수익률은 비중과 수익률의 단순 가중평균합이 됩니다.
포트폴리오의 분산은 w'Ωw로 나타낼 수 있으며, 해당값에 sqrt()를 치해주면 포트폴리오의 표준편차가 됩니다.
마지막으로 샤프지수는 포트폴리오 수익률에서 무위험 수익률을 차감한 값을, 포트폴리오 변동성으로 나눠줍니다. 무위험 수익률은 5%로 가정합니다.

해당 값들을 c()를 이용해 벡터로 묶어준 후 sim 리스트에 저장합니다.
시뮬레이션이 끝나면 do.call() 함수를 이용해 행으로 묶습니다.

target_minvol = sim[which.min(sim$Volatility), ]
target_maxsharpe = sim[which.max(sim$Sharpe), ]

> target_minvol
      삼성전자   셀트리온    현대차  신한지주  SK텔레콤    Return Volatility    Sharpe
4657 0.2447684 0.09979261 0.0992029 0.1240248 0.4322112 0.1306394  0.1896045 0.4253035
> target_maxsharpe
      삼성전자  셀트리온    현대차    신한지주   SK텔레콤    Return Volatility    Sharpe
3457 0.3303538 0.5068116 0.1004232 0.009327201 0.05308417 0.3044934  0.3103827 0.8199342

효율적 경계선에서 중요한 지점은 크게 최소분산 지점과 최대샤프 지점입니다. 최소분산 지점은 which.min() 함수를 통해 표준편차가 가장 작은 지점, 최대샤프 지점은 which.max() 함수를 통해 샤프지수가 최대인 지점을 찾으면 됩니다.

# plot simulation
p = sim %>%
  ggplot(aes(x = `Volatility`, y = `Return`, color = `Sharpe`)) +
  geom_point() +
  scale_color_gradient(low = 'blue', high = 'red') +
  theme_classic()

p

ggplot을 이용해 시뮬레이션 결과를 그려주도록 합니다. aes() 내부에 color 인자를 Sharpe로 지정하여 샤프지수에 따라 색이 다르며, 샤프지수가 높을수록 붉게 표시됩니다.

# plot min vol and max sharpe
p2 = p +
  geom_point(data = target_minvol,
             aes(x = Volatility, y = Return),
             col = 'black', shape = 18, size = 3) +
  geom_point(data = target_maxsharpe,
             aes(x = Volatility, y = Return),
             col = 'black', shape = 18, size = 3) +
  annotate('text',
           x = target_minvol$Volatility + 0.03,
           y = target_minvol$Return,
           label='Min Vol',
           fontface=2) +
  annotate('text',
           x = target_maxsharpe$Volatility + 0.03,
           y = target_maxsharpe$Return,
          label='Max Sharpe',
          fontface=2) +
  geom_point(data = yearly_stat,
             aes(x = yearly_vol, y = yearly_mean),
             color = 'black', shape = 4, size = 3)
p2

앞선 그림에 추가적으로 최소분산 및 최대 샤프지수 지점을 표시해 주었습니다. 또한 검은색 X는 개별 주식의 변동성과 수익률에 해당하는 부분으로써, 포트폴리오 구성을 통해 훨씬 효율적인 투자가 가능해짐이 보입니다.

# plot CML

p2 + 
  geom_abline(intercept = 0.05,
              slope = target_maxsharpe$Sharpe,
              size = 1) +
  expand_limits(x=0) +
  expand_limits(y=0) 

마지막으로 무위험 수익률과 최대샤프지수 지점를 이은 CML을 그려줍니다. 인터셉트는 무위험 수익률, 기울기는 샤프지수가 됩니다. Max Sharpe 지점에 투자하는 것이 가장 이상적인 포트폴리오이며, 투자자의 risk averse에 따라 직선의 왼쪽이나 오른쪽으로 이동하며 위험자산과 무위험자산의 투자 비중을 조절하게 됩니다.