lwShrink <- function(x, shrink=NULL){
# port of matlab code from http://www.econ.uzh.ch/faculty/wolf/publications.html#9
# Ledoit, O. and Wolf, M. (2004).
# Honey, I shrunk the sample covariance matrix.
# Journal of Portfolio Management 30, Volume 4, 110-119.
# De-mean returns
n <- nrow(x)
p <- ncol(x)
meanx <- colMeans(x)
x <- x - matrix(rep(meanx, n), ncol=p, byrow=TRUE)
# Compute sample covariance matrix using the de-meaned returns
sample <- (1 / n) * (t(x) %*% x)
# Compute prior
var <- matrix(diag(sample), ncol=1)
sqrtvar <- sqrt(var)
tmpMat <- matrix(rep(sqrtvar, p), nrow=p)
rBar <- (sum(sum(sample / (tmpMat * t(tmpMat)))) - p) / (p * (p - 1))
prior <- rBar * tmpMat * t(tmpMat)
diag(prior) <- var
if(is.null(shrink)){
# What is called pi-hat
y <- x^2
phiMat <- t(y) %*% y / n - 2 * (t(x) %*% x) * sample / n + sample^2
phi <- sum(phiMat)
# What is called rho-hat
term1 <- (t(x^3) %*% x) / n
help <- t(x) %*% x / n
helpDiag <- matrix(diag(help), ncol=1)
term2 <- matrix(rep(helpDiag, p), ncol=p, byrow=FALSE) * sample
term3 <- help * matrix(rep(var, p), ncol=p, byrow=FALSE)
term4 <- matrix(rep(var, p), ncol=p, byrow=FALSE) * sample
thetaMat <- term1 - term2 - term3 + term4
diag(thetaMat) <- 0
rho <- sum(diag(phiMat)) + rBar * sum(sum(((1 / sqrtvar) %*% t(sqrtvar)) * thetaMat))
# What is called gamma-hat
gamma <- norm(sample - prior, "F")^2
# Compute shrinkage constant
kappa <- (phi - rho) / gamma
shrinkage <- max(0, min(1, kappa / n))
} else {
shrinkage <- shrink
}
# Compute the estimator
sigma <- shrinkage * prior + (1 - shrinkage) * sample
out <- list(cov=sigma, prior=prior, shrinkage=shrinkage)
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.