library(xts)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
library(portfolioBacktest)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(riskParityPortfolio)

# download price data
faang_data <- stockDataDownload(c("BA", "FDX", "GOOG", "NFLX", "AAPL", "AMZN", "NKE"),
                                from = "2016-04-03", to = "2021-04-03")
## Loading stock data from local file /Users/h/Documents/hu/515/project/stockdata_from_2016-04-03_to_2021-04-03_(7dc03631429d9c651e28337561272a5a).RData
risk_parity <- function(dataset) {
  prices <- dataset$adjusted
  log_returns <- diff(log(prices))[-1]
  return(riskParityPortfolio(cov(log_returns))$w)
}

# tangency portfolio (maximum sharpe ratio)
library(quadprog)
max_sharpe_ratio <- function(dataset) {
    prices <- dataset$adjusted
    log_returns <- diff(log(prices))[-1]
    N <- ncol(prices)
    Sigma <- cov(log_returns)
    mu <- colMeans(log_returns)
    if (all(mu <= 1e-8))
        return(rep(0, N))
    Dmat <- 2 * Sigma
    Amat <- diag(N)
    Amat <- cbind(mu, Amat)
    bvec <- c(1, rep(0, N))
    dvec <- rep(0, N)
    res <- solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec, meq = 1)
    w <- res$solution
    return(w/sum(w))
}

# call portfolioBacktest and benchmark against the uniform (1/N) portfolio
bt <- portfolioBacktest(list("risk parity portfolio" = risk_parity,
                             "tangency portfolio"    = max_sharpe_ratio),
                        list(faang_data),
                        T_rolling_window = 12*20, 
                        optimize_every = 3*20, rebalance_every = 3*20)
## 
##  [Backtesting 2 portfolios over 1 datasets (periodicity = daily data)]
# dates of the designed portfolios
index(bt$tangency$data1$w_designed)
##  [1] "2017-03-15" "2017-06-09" "2017-09-05" "2017-11-29" "2018-02-27"
##  [6] "2018-05-23" "2018-08-17" "2018-11-12" "2019-02-11" "2019-05-08"
## [11] "2019-08-02" "2019-10-28" "2020-01-24" "2020-04-21" "2020-07-16"
## [16] "2020-10-09" "2021-01-06"
# check performance summary
backtestSummary(bt)$performance
##                   risk parity portfolio tangency portfolio
## Sharpe ratio               1.115135e+00       1.023492e+00
## max drawdown               3.436311e-01       3.142726e-01
## annual return              2.762041e-01       2.936504e-01
## annual volatility          2.476867e-01       2.869103e-01
## Sterling ratio             8.037806e-01       9.343813e-01
## Omega ratio                1.236890e+00       1.227325e+00
## ROT (bps)                  8.213921e+03       9.556370e+02
## VaR (0.95)                 2.371748e-02       2.757762e-02
## CVaR (0.95)                3.958884e-02       4.500072e-02
## cpu time                   8.000000e-03       3.823529e-03
## failure rate               0.000000e+00       0.000000e+00
backtestChartCumReturns(bt)

backtestChartDrawdown(bt)

backtestChartStackedBar(bt, portfolio = "risk parity portfolio", legend = TRUE)

backtestChartStackedBar(bt, portfolio = "tangency portfolio" , legend = TRUE)

Comparison with other Packages

library(FinCovRegularization)
library(cccp)
library(RiskPortfolios)

# generate synthetic data
set.seed(42)
N <- 10
V <- matrix(rnorm(N*(N+N/5)), N+N/5, N)
Sigma <- cov(V)

# uniform initial guess for the portfolio weights
w0 <- rep(1/N, N)

# compute risk parity portfolios using different methods
rpp <- riskParityPortfolio(Sigma, w0 = w0, formulation = "rc-double-index")
## Warning in riskParityPortfolio(Sigma, w0 = w0, formulation = "rc-double-index"):
## The problem is a vanilla risk parity portofolio, but a nonconvex formulation has
## been chosen. Consider not specifying the formulation argument in order to use
## the convex formulation and get a guaranteed global solution.
#> Warning in riskParityPortfolio(Sigma, w0 = w0, formulation = "rc-double-index"): The problem is a vanilla risk parity
#> portofolio, but a nonconvex formulation has been chosen. Consider not specifying the formulation argument in order to
#> use the convex formulation and get a guaranteed global solution.
fincov_w <- FinCovRegularization::RiskParity(Sigma)
riskport_w <- optimalPortfolio(Sigma = Sigma,
                               control = list(type = "erc",
                                              constraint = "lo"))
cccp_w <- c(getx(cccp::rp(w0, Sigma, mrc = w0, optctrl = ctrl(trace = FALSE))))

# plot
w_all <- cbind("riskParityPortfolio"  = rpp$w,
               "FinCovRegularization" = fincov_w,
               "cccp"                 = cccp_w,
               "RiskPortfolios"       = riskport_w)
barplotPortfolioRisk(w_all, Sigma)

computational time

library(microbenchmark)
library(ggplot2)
library(cccp)
library(RiskPortfolios)
library(riskParityPortfolio)


b <- rep(1/N, N)

op <- microbenchmark(
          cccp = rp(b, Sigma, b, optctrl = ctrl(trace = FALSE)),
          riskParityPortfolio = 
            riskParityPortfolio:::risk_parity_portfolio_ccd_spinu(Sigma, b, 1e-6, 50),
          RiskPortfolios = 
            optimalPortfolio(Sigma = Sigma,
                             control = list(type = 'erc', constraint = 'lo')),
          times = 10L)
print(op)
## Unit: microseconds
##                 expr      min       lq      mean    median       uq       max
##                 cccp 1251.313 1376.285 1700.7725 1518.3305 1585.527  3816.827
##  riskParityPortfolio   17.058   26.919   43.5547   40.5285   42.708   115.056
##       RiskPortfolios 6657.256 6675.168 7520.6573 6836.0285 7583.377 11153.130
##  neval
##     10
##     10
##     10
autoplot(op)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

https://cran.r-project.org/web/packages/riskParityPortfolio/vignettes/RiskParityPortfolio.html