Asset Allocation ModelをRで組んでみた。

gmv <- function(r_dat,r_cov){
  library(MASS)
  i <- matrix(1,NCOL(r_dat),1)
  r_weight <- (ginv(r_cov)%*%i)/as.numeric(t(i)%*%ginv(r_cov)%*%i)
  wr_dat <- r_dat*as.numeric(r_weight)
  portfolio <- apply(wr_dat,1,sum)
  pr_dat <- data.frame(wr_dat,portfolio)
  sd <- sd(portfolio)
  result <- list(r_weight,pr_dat,sd)
  names(result) <- c("weight","return","portfolio risk") 
  return(result)
}

nlgmv <- function(r_dat,r_cov){
  qp.out <- solve.QP(Dmat=r_cov,dvec=rep(0,NCOL(r_dat)),Amat=cbind(rep(1,NCOL(r_dat)),diag(NCOL(r_dat))),
           bvec=c(1,rep(0,NCOL(r_dat))),meq=1)
  r_weight <- qp.out$solution
  wr_dat <- r_dat*r_weight
  portfolio <- apply(wr_dat,1,sum)
  pr_dat <- data.frame(wr_dat,portfolio)
  sd <- sd(portfolio)
  result <- list(r_weight,pr_dat,sd)
  names(result) <- c("weight","return","portfolio risk")
  return(result)
}

carroll <- function(r_dat,FLG){
  
  library(rmgarch)
  
  if(FLG == "benchmark"){
    H <- cov(r_dat)
  }else{
    #1. define variables
    N <- NCOL(r_dat) # the number of assets
    
    #2. estimate covariance matrix
    basic_garch = ugarchspec(mean.model = list(armaOrder = c(0, 0),include.mean=TRUE), variance.model = list(garchOrder = c(1,1), model = 'sGARCH'), distribution.model = 'norm')
    multi_garch = multispec(replicate(N, basic_garch))
    dcc_set = dccspec(uspec = multi_garch, dccOrder = c(1, 1), distribution = "mvnorm",model = "DCC")
    fit_dcc_garch = dccfit(dcc_set, data = r_dat, fit.control = list(eval.se = TRUE))
    forecast_dcc_garch <- dccforecast(fit_dcc_garch)
    if (FLG == "CCC"){
      #Constant conditional correlation (CCC) model
      D <- sigma(forecast_dcc_garch)
      R_ccc <- cor(r_dat)
      H <- diag(D[,,1])%*%R_ccc%*%diag(D[,,1])
      colnames(H) <- colnames(r_dat)
      rownames(H) <- colnames(r_dat)
    }
    else{
      #Dynamic Conditional Correlation (DCC) model
      H <- as.matrix(rcov(forecast_dcc_garch)[[1]][,,1])
      if (FLG == "DECO"){
        #Dynamic Equicorrelation (DECO) model
        one <- matrix(1,N,N)
        iota <- rep(1,N)
        Q_dcc <- rcor(forecast_dcc_garch,type="Q")[[1]][,,1]
        rho <- as.vector((N*(N-1))^(-1)*(t(iota)%*%Q_dcc%*%iota-N))
        D <- sigma(forecast_dcc_garch)
        R_deco <- (1-rho)*diag(1,N,N) + rho*one
        H <- diag(D[,,1])%*%R_deco%*%diag(D[,,1])
        colnames(H) <- colnames(r_dat)
        rownames(H) <- colnames(r_dat)
      }
    }
  }
  return(H)
}

library(quantmod)

#**************************
# ★8 ASSETS SIMULATION
# SPY - S&P 500 
# QQQ - Nasdaq 100
# EEM - Emerging Markets
# IWM - Russell 2000
# EFA - EAFE
# TLT - 20 Year Treasury
# IYR - U.S. Real Estate
# GLD - Gold
#**************************

# load historical prices from Yahoo Finance
symbol.names = c("S&P 500","Nasdaq 100","Emerging Markets","Russell 2000","EAFE","20 Year Treasury","U.S. Real Estate","Gold")
symbols = c("SPY","QQQ","EEM","IWM","EFA","TLT","IYR","GLD")
getSymbols(symbols, from = '1980-01-01', auto.assign = TRUE)

#gn dates for all symbols & convert to monthly
hist.prices = merge(SPY,QQQ,EEM,IWM,EFA,TLT,IYR,GLD)
month.ends = endpoints(hist.prices, 'day')
hist.prices = Cl(hist.prices)[month.ends, ]
colnames(hist.prices) = symbols

# remove any missing data
hist.prices = na.omit(hist.prices['1995::'])

# compute simple returns
hist.returns = na.omit( ROC(hist.prices, type = 'discrete') )

# compute historical returns, risk, and correlation
ia = list()
ia$expected.return = apply(hist.returns, 2, mean, na.rm = T)
ia$risk = apply(hist.returns, 2, sd, na.rm = T)
ia$correlation = cor(hist.returns, use = 'complete.obs', method = 'pearson')

ia$symbols = symbols
ia$symbol.names = symbol.names
ia$n = length(symbols)
ia$hist.returns = hist.returns

# convert to annual, year = 12 months
annual.factor = 12
ia$expected.return = annual.factor * ia$expected.return
ia$risk = sqrt(annual.factor) * ia$risk

rm(SPY,QQQ,EEM,IWM,EFA,TLT,IYR,GLD)

PerformanceAnalytics::charts.PerformanceSummary(hist.returns, main = "パフォーマンスサマリー")

# BACK TEST
backtest <- function(r_dat,FLG,start_date,span,learning_term,port){
  #-----------------------------------------
  # BACKTEST
  # r_dat - return data(xts object) 
  # FLG - flag(CCC,DCC,DECO)
  # start_date - start date for backtest
  # span - rebalance frequency
  # learning_term - learning term (days)
  # port - method of portfolio optimization
  #-----------------------------------------
  
  library(stringi)

  initial_dat <- r_dat[stri_c(as.Date(start_date)-learning_term,"::",as.Date(start_date))]
  for (i in NROW(initial_dat):NROW(r_dat)) {
    if (i == NROW(initial_dat)){
      H <- carroll(initial_dat[1:(NROW(initial_dat)-1),],FLG)
      if (port == "nlgmv"){
        result <- nlgmv(initial_dat,H)
      }else if (port == "risk parity"){
        result <- risk_parity(initial_dat,H)
      }
      weight <- t(result$weight)
      colnames(weight) <- colnames(initial_dat)
      p_return <- initial_dat[NROW(initial_dat),]*result$weight
    } else {
      if (i %in% endpoints(r_dat,span)){
        H <- carroll(test_dat[1:(NROW(test_dat)-1),],FLG)
        if (port == "nlgmv"){
          result <- nlgmv(test_dat,H)
        }else if (port == "risk parity"){
          result <- risk_parity(test_dat,H)
        }
        
      }
      weight <- rbind(weight,t(result$weight))
      p_return <- rbind(p_return,test_dat[NROW(test_dat),]*result$weight)
    }
    if (i != NROW(r_dat)){
      term <- stri_c(index(r_dat[i+1,])-learning_term,"::",index(r_dat[i+1,])) 
      test_dat <- r_dat[term]
    }
  }
  p_return$portfolio <- xts(apply(p_return,1,sum),order.by = index(p_return))
  weight.xts <- xts(weight,order.by = index(p_return))

  result <- list(p_return,weight.xts)
  names(result) <- c("return","weight")
  return(result)
}

CCC <- backtest(hist.returns,"CCC","2007-01-04","months",365,"risk parity")
DCC <- backtest(hist.returns,"DCC","2007-01-04","months",365,"risk parity")
DECO <- backtest(hist.returns,"DECO","2007-01-04","months",365,"risk parity")
benchmark <- backtest(hist.returns,"benchmark","2007-01-04","months",365,"risk parity")

result <- merge(CCC$return$portfolio,DCC$return$portfolio,DECO$return$portfolio,benchmark$return$portfolio)
colnames(result) <- c("CCC","DCC","DECO","benchmark")

PerformanceAnalytics::charts.PerformanceSummary(result,main = "BACKTEST")

# plot allocation weighting
d_allocation <- function(ggweight,title){
  #install.packages("tidyverse")
  library(tidyverse)
  ggweight <- gather(ggweight,key=ASSET,value=weight,-Date,-method)
  ggplot(ggweight, aes(x=Date, y=weight,fill=ASSET)) +
    geom_area(colour="black",size=.1) +
    scale_y_continuous(limits = c(0,1)) +
    labs(title=title) + facet_grid(method~.)
}

gmv_weight <- rbind(data.frame(CCC$weight,method="CCC",Date=index(CCC$weight)),data.frame(DCC$weight,method="DCC",Date=index(DCC$weight)),data.frame(DECO$weight,method="DECO",Date=index(DECO$weight)))

# plot allocation weighting
d_allocation(gmv_weight,"GMV Asset Allocation")

risk_parity <- function(r_dat,r_cov){
  fn <- function(weight, r_cov) {
    N <- NROW(r_cov)
    risks <-  weight * (r_cov %*% weight)
    g <- rep(risks, times = N) - rep(risks, each = N)
  return(sum(g^2))
  }
  dfn <- function(weight,r_cov){
    out <- weight
    for (i in 0:length(weight)) {
      up <- dn <- weight
      up[i] <- up[i]+.0001
      dn[i] <- dn[i]-.0001
      out[i] = (fn(up,r_cov) - fn(dn,r_cov))/.0002
    }
    return(out)
  }
  std <- sqrt(diag(r_cov)) 
  x0 <- 1/std/sum(1/std)
  res <- nloptr::nloptr(x0=x0,
                 eval_f=fn,
                 eval_grad_f=dfn,
                 eval_g_eq=function(weight,r_cov) { sum(weight) - 1 },
                 eval_jac_g_eq=function(weight,r_cov) { rep(1,length(std)) },
                 lb=rep(0,length(std)),ub=rep(1,length(std)),
                 opts = list("algorithm"="NLOPT_LD_SLSQP","print_level" = 0,"xtol_rel"=1.0e-8,"maxeval" = 1000),
                 r_cov = r_cov)
  r_weight <- res$solution
  names(r_weight) <- colnames(r_cov)
  wr_dat <- r_dat*r_weight
  portfolio <- apply(wr_dat,1,sum)
  pr_dat <- data.frame(wr_dat,portfolio)
  sd <- sd(portfolio)
  result <- list(r_weight,pr_dat,sd)
  names(result) <- c("weight","return","portfolio risk")
  return(result)
  }

CCC <- backtest(hist.returns,"CCC","2007-01-04","months",365,"risk parity")
DCC <- backtest(hist.returns,"DCC","2007-01-04","months",365,"risk parity")
DECO <- backtest(hist.returns,"DECO","2007-01-04","months",365,"risk parity")
benchmark <- backtest(hist.returns,"benchmark","2007-01-04","months",365,"risk parity")

result <- merge(CCC$return$portfolio,DCC$return$portfolio,DECO$return$portfolio,benchmark$return$portfolio)
colnames(result) <- c("CCC","DCC","DECO","benchmark")

PerformanceAnalytics::charts.PerformanceSummary(result, main = "BACKTEST COMPARISON")

library(plotly)

# plot allocation weighting
riskparity_weight <- rbind(data.frame(CCC$weight,method="CCC",Date=index(CCC$weight)),data.frame(DCC$weight,method="DCC",Date=index(DCC$weight)),data.frame(DECO$weight,method="DECO",Date=index(DECO$weight)),data.frame(benchmark$weight,method="benchmark",Date=index(benchmark$weight)))

PerformanceAnalytics::charts.PerformanceSummary(result, main = "BACKTEST COMPARISON")

riskparity_weight <- rbind(data.frame(CCC$weight,method="CCC",Date=index(CCC$weight)),data.frame(DCC$weight,method="DCC",Date=index(DCC$weight)),data.frame(DECO$weight,method="DECO",Date=index(DECO$weight)),data.frame(benchmark$weight,method="benchmark",Date=index(benchmark$weight)))

# plot allocation weighting
d_allocation(riskparity_weight, "Risk Parity Asset Allocation")