Example of how to use Newton’s method to calculate risk parity weights. The algorithm is described in Efficient Algorithms for Computing Risk Parity Portfolio Weights
Set up F(y), system of equations and its jacobian matrix.
# F(y), system of equations
eval_f <- function(x, sigma, lambda){
# x = weights vector, sigma = covariance matrix, lambda = a constant to be found
# sets up n+1 x 1 matrix of equations, where n is the amount of assets or variables
# and the +1 is the restriction that the weights sum to 1
x <- as.vector(x)
x <- x[1:nrow(sigma)]
f0 <- matrix(nrow = NROW(sigma) + 1, ncol = 1)
f0[1:NROW(sigma), 1] <- (sigma %*% x) - (lambda * 1/x)
f0[(NROW(sigma) + 1), 1] <- sum(x) - 1
return(f0)
}
# Jacobian matrix of F(y)
jacob_f <- function(x, sigma, lambda){
# x = weights vector, sigma = covariance matrix, lambda = a constant to be found
# sets up n+1 x n+1 jacobian matrix of F(y)
x <- as.vector(x)
x <- x[1:NROW(sigma)]
g <- matrix(nrow = NROW(sigma) + 1, ncol = NCOL(sigma)+1)
g[1:NROW(sigma), 1:NCOL(sigma)] <- sigma + as.vector(lambda) * diag(1/x^2)
g[1:NROW(sigma), (NCOL(sigma)+1)] <- 1/x
g[(NROW(sigma)+1), 1:NCOL(sigma)] <- rep(1, ncol(sigma))
g[(NROW(sigma)+1), (NCOL(sigma)+1)] <- 0
return(g)
}Get returns for a hypothetical sector portfolio to test the algorithm
# requires fImport package to get stock returns
if(require(fImport) == FALSE){
install.packages("fImport")
require(fImport)
}
# use fidelity sector ETFs as example
assets <- c("FDIS", "FSTA", "FENY", "FNCL", "FHLC",
"FIDU", "FTEC", "FMAT", "FCOM", "FUTY")
nAssets <- length(assets)
# get historical price data (for past year) from yahoo finance data base
# prices are in every 6th column
prices <- yahooSeries(assets)[, seq(6, nAssets*6, 6)]
colnames(prices) <- c("ConsumerDisc", "ConsumerStap", "Energy",
"Financials", "HealthCare", "Industrials",
"InfoTech", "Materials", "Telecomm", "Utilities")
head(prices)## GMT
## ConsumerDisc ConsumerStap Energy Financials HealthCare
## 2014-01-06 26.38 25.06 24.87 25.59 26.07
## 2014-01-07 26.54 25.19 25.06 25.62 26.37
## 2014-01-08 26.50 25.01 24.87 25.68 26.63
## 2014-01-09 26.44 25.12 24.80 25.74 26.85
## 2014-01-10 26.52 25.14 24.83 25.75 27.00
## 2014-01-13 25.98 24.98 24.38 25.40 26.79
## Industrials InfoTech Materials Telecomm Utilities
## 2014-01-06 26.58 26.38 25.67 24.82 23.65
## 2014-01-07 26.71 26.62 25.64 24.94 23.83
## 2014-01-08 26.67 26.65 25.75 24.91 23.73
## 2014-01-09 26.78 26.53 25.67 24.54 23.86
## 2014-01-10 26.85 26.59 25.75 24.62 24.15
## 2014-01-13 26.54 26.26 25.43 24.39 23.94# convert prices to continuously compounded returns
ret <- timeSeries::returns(prices, method = "continuous")
head(ret)## GMT
## ConsumerDisc ConsumerStap Energy Financials
## 2014-01-07 0.006046882 0.0051741409 0.007610692 0.0011716463
## 2014-01-08 -0.001508296 -0.0071713455 -0.007610692 0.0023391824
## 2014-01-09 -0.002266718 0.0043885967 -0.002818605 0.0023337233
## 2014-01-10 0.003021150 0.0007958616 0.001208946 0.0003884249
## 2014-01-13 -0.020572154 -0.0063846985 -0.018289475 -0.0136854531
## 2014-01-14 0.009576783 0.0047923414 0.012230075 0.0070616221
## HealthCare Industrials InfoTech Materials
## 2014-01-07 0.011441772 0.004878974 0.009056666 -0.001169363
## 2014-01-08 0.009811399 -0.001498689 0.001126338 0.004280995
## 2014-01-09 0.008227421 0.004116002 -0.004512982 -0.003111632
## 2014-01-10 0.005571045 0.002610481 0.002259037 0.003111632
## 2014-01-13 -0.007808182 -0.011612792 -0.012488336 -0.012505048
## 2014-01-14 0.014084740 0.010494849 0.019608471 0.013669418
## Telecomm Utilities
## 2014-01-07 0.004823160 0.0075821762
## 2014-01-08 -0.001203611 -0.0042052207
## 2014-01-09 -0.014964890 0.0054633462
## 2014-01-10 0.003254681 0.0120809634
## 2014-01-13 -0.009385909 -0.0087336800
## 2014-01-14 0.007759892 0.0008350731# covariance matrix of returns
sigma <- cov(ret)Run algorithm (Newton’s method) and test results
maxIterations <- 100
tol <- 1e-6
# time algo
timeStart <- proc.time()
for(i in 1:maxIterations){
if (i == 1){
# 1/N first weight guess
x_i <- rep(1/nAssets, nAssets)
# lambda = portfolio standard deviation / number of assets
lambda <- t(x_i) %*% sigma %*% x_i / nAssets
}
else{
# all other iterations after 1st guess
jInv <- solve(jacob_f(c(x_i, 1), sigma, lambda))
Fx <- eval_f(c(x_i, 1), sigma, lambda)
ans <- c(x_i, 1) - jInv %*% Fx
# if solution is found exit the loop
if(norm(c(x_i,1) - ans) <= tol) break
# otherwise keep trying new weight vector guess
x_i <- ans[1:nAssets]
lambda <- t(x_i) %*% sigma %*% x_i / nAssets
}
}
# time elapsed
proc.time() - timeStart## user system elapsed
## 0.004 0.000 0.004# check solution
ans <- ans[1:nAssets]
num <- t(ans) %*% sigma
denom <- t(ans) %*% sigma %*% ans
cctr <- num / denom[1,1] * ans
cctr <- t(cctr)
colnames(cctr) <- "Percent.Contribution.to.Risk"
cctr## Percent.Contribution.to.Risk
## ConsumerDisc 0.1
## ConsumerStap 0.1
## Energy 0.1
## Financials 0.1
## HealthCare 0.1
## Industrials 0.1
## InfoTech 0.1
## Materials 0.1
## Telecomm 0.1
## Utilities 0.1