#' @name covEstimation
#' @aliases covEstimation
#' @title Covariance matrix estimation
#' @description Function which performs various estimations of covariance matrices.
#' @details The argument \code{control} is a list that can supply any of the following
#' components:
#' \itemize{
#' \item \code{type} method used to compute the
#' covariance matrix, among \code{'naive'}, \code{'ewma'}, \code{'lw'},
#' \code{'factor'},\code{'const'}, \code{'cor'}, \code{'oneparm'},
#' \code{'diag'} and \code{'large'} where:
#'
#' \code{'naive'} is used to compute
#' the naive (standard) covariance matrix.
#'
#' \code{'ewma'} is used to compute the exponential weighting moving average covariance matrix. The following formula is used
#' to compute the ewma covariance matrix:
#' \deqn{\Sigma_t := \lambda \Sigma_{t-1} + (1-\lambda)r_{t-1}r_{t-1}}{Sigma[t]
#' := lambda * Sigma[t-1] + (1-lambda) r[t-1]'r[t-1]}
#' where \eqn{r_t} is the \eqn{(N \times 1)}{(N x 1)} vector of returns at time
#' \eqn{t}. Note that the data must be sorted from the oldest to the latest. See RiskMetrics (1996)
#'
#' \code{'factor'} is used to compute the covariance matrix estimation using a
#' K-factor approach. See Harman (1976).
#'
#' \code{'lw'} is a weighted average of the sample covariance matrix and a
#' 'prior' or 'shrinkage target'. The prior is given by a one-factor model and
#' the factor is equal to the cross-sectional average of all the random
#' variables. See Ledoit and Wolf (2003).
#'
#' \code{'const'} is a weighted average of the sample covariance matrix and a
#' 'prior' or 'shrinkage target'. The prior is given by constant correlation
#' matrix. See Ledoit and Wolf (2002).
#'
#' \code{'cor'} is a weighted average of the sample covariance matrix and a
#' 'prior' or 'shrinkage target'. The prior is given by the constant
#' correlation covariance matrix given by Ledoit and Wolf (2003).
#'
#' \code{'oneparm'} is a weighted average of the sample covariance matrix and a
#' 'prior' or 'shrinkage target'. The prior is given by the one-parameter
#' matrix. All variances are the same and all covariances are zero.
#' See Ledoit and Wolf (2004).
#'
#' \code{'diag'} is a weighted average of the sample covariance matrix and a
#' 'prior' or 'shrinkage target'. The prior is given by a diagonal matrix.
#' See Ledoit and Wolf (2002).
#'
#' \code{'large'} This estimator is a weighted average of the sample covariance
#' matrix and a 'prior' or 'shrinkage target'. Here, the prior is given by a
#' one-factor model. The factor is equal to the cross-sectional average of all
#' the random variables. The weight, or 'shrinkage intensity' is chosen to
#' minimize quadratic loss measured by the Frobenius norm. The estimator is
#' valid as the number of variables and/or the number of observations go to
#' infinity, but Monte-Carlo simulations show that it works well for values as
#' low as 10. The main advantage is that this estimator is guaranteed to be
#' invertible and well-conditioned even if variables outnumber observations.
#' See Ledoit and Wolf (2004).
#'
#' \code{'bs'} is the Bayes-Stein estimator for the covariance matrix given by
#' Jorion (1986).
#'
#' Default: \code{type = 'naive'}.
#'
#' \item \code{lambda} decay parameter. Default: \code{lambda = 0.94}.
#'
#' \item \code{K} number of factors to use when the K-factor approach is
#' chosen to estimate the covariance matrix. Default: \code{K = 1}.}
#'
#' @param rets a matrix \eqn{(T \times N)}{(T x N)} of returns.
#' @param control control parameters (see *Details*).
#' @return A \eqn{(N \times N)}{(N x N)} covariance matrix.
#' @note Part of the code is adapted from the Matlab code by Ledoit and Wolf (2014).
#' @author David Ardia, Kris Boudt and Jean-Philippe Gagnon Fleury.
#' @references
#' Jorion, P. (1986).
#' Bayes-Stein estimation for portfolio analysis.
#' \emph{Journal of Financial and Quantitative Analysis} \bold{21}(3), pp.279-292.
#'
#' Harman, H.H. (1976)
#' \emph{Modern Factor Analysis}.
#' 3rd Ed. Chicago: University of Chicago Press.
#'
#' Ledoit, O., Wolf, M. (2002).
#' Improved estimation of the covariance matrix of stock returns with an application to portfolio selection.
#' \emph{Journal of Empirical Finance} \bold{10}(5), pp.603-621.
#'
#' Ledoit, O., Wolf, M. (2003).
#' Honey, I Shrunk the Sample Covariance Matrix.
#' \emph{Journal of Portfolio Management} \bold{30}(4), pp.110-119.
#'
#' Ledoit, O., Wolf, M. (2004).
#' A well-conditioned estimator for large-dimensional covariance matrices.
#' \emph{Journal of Multivariate Analysis} \bold{88}(2), pp.365-411.
#'
#' RiskMetrics (1996)
#' \emph{RiskMetrics Technical Document}.
#' J. P. Morgan/Reuters.
#' @keywords htest
#' @examples
#' # Load returns of assets or portfolios
#' data("Industry_10")
#' rets = Industry_10
#'
#' # Naive covariance estimation
#' covEstimation(rets)
#'
#' # Ewma estimation of the covariance with default lambda = 0.94
#' covEstimation(rets, control = list(type = 'ewma'))
#'
#' # Ewma estimation of the covariance with default lambda = 0.90
#' covEstimation(rets, control = list(type = 'ewma', lambda = 0.9))
#'
#' # Factor estimation of the covariance with dafault K = 1
#' covEstimation(rets, control = list(type = 'factor'))
#'
#' # Factor estimation of the covariance with K = 3
#' covEstimation(rets, control = list(type = 'factor', K = 3))
#'
#' # Ledot-Wolf's estimation of the covariance
#' covEstimation(rets, control = list(type = 'lw'))
#'
#' # Shrinkage of the covariance matrix using constant correlation matrix
#' covEstimation(rets, control = list(type = 'const'))
#'
#' # Shrinkage of the covariance matrix towards constant correlation matrix by
#' # Ledoit-Wolf.
#' covEstimation(rets, control = list(type = 'cor'))
#'
#' # Shrinkage of the covariance matrix towards one-parameter matrix
#' covEstimation(rets, control = list(type = 'oneparm'))
#'
#' # Shrinkage of the covaraince matrix towards diagonal matrix
#' covEstimation(rets, control = list(type = 'diag'))
#'
#' # Shrinkage of the covariance matrix for large data set
#' covEstimation(rets, control = list(type = 'large'))
#' @export
#' @importFrom stats cor cov factanal lm median quantile sd
covEstimation <- function(rets, control = list()) {
if (missing(rets)) {
stop("rets is missing")
}
if (!is.matrix(rets)) {
stop("rets must be a (T x N) matrix")
}
ctr <- .ctrCov(control)
if (ctr$type[1] == "naive") {
Sigma <- .naiveCov(rets = rets)
} else if (ctr$type[1] == "ewma") {
Sigma <- .ewmaCov(rets = rets, lambda = ctr$lambda)
} else if (ctr$type[1] == "lw") {
Sigma <- .lwCov(rets = rets)
} else if (ctr$type[1] == "factor") {
Sigma <- .factorCov(rets = rets, K = ctr$K)
} else if (ctr$type[1] == "const") {
Sigma <- .constCov(rets = rets)
} else if (ctr$type[1] == "cor") {
Sigma <- .corCov(rets = rets)
} else if (ctr$type[1] == "oneparm") {
Sigma <- .oneparmCov(rets = rets)
} else if (ctr$type[1] == "diag") {
Sigma <- .diagCov(rets = rets)
} else if (ctr$type[1] == "large") {
Sigma <- .largeCov(rets = rets)
} else if (ctr$type[1] == "bs") {
Sigma <- .bsCov(rets = rets)
} else {
stop("control$type is not well defined")
}
return(Sigma)
}
.ctrCov <- function(control = list()) {
#Function used to control the list input INPUTS control : a control
#list The argument control is a list that can supply any of the
#following components type : default = 'naive' lambda : default = 0.94
#K : default = 1 OUTPUTs control : a list
if (!is.list(control)) {
stop("control must be a list")
}
if (length(control) == 0) {
control <- list(type = "naive", lambda = 0.94, K = 1)
}
nam <- names(control)
if (!("type" %in% nam) || is.null(control$type)) {
control$type <- c("naive", "ewma", "lw", "factor", "rtm", "const",
"cor", "oneparm", "diag", "large", "bs")
}
if (!("lambda" %in% nam) || is.null(control$lambda)) {
control$lambda <- 0.94
}
if (!("K" %in% nam) || is.null(control$K)) {
control$K <- 1
}
return(control)
}
.naiveCov <- function(rets) {
## Compute the naive covariance matrix INPUTs rets : matrix (T x N)
## returns OUTPUTs Sigma : matrix (N x N) covariance matrix DA here we
## could check that the dimension doesn't lead to bullshit outputs
Sigma <- cov(rets)
return(Sigma)
}
.ewmaCov <- function(rets, lambda) {
## Compute the exponential weighted moving average covariance matrix
## INPUTs rets : matrix (T x N) returns lambda : OUTPUTs Sigma : matrix
## (N x N) covariance matrix DA check speed here (use just-in-time
## compilation or compiler package)
# # Infinite sample t = nrow(rets) n = ncol(rets) Sigma = matrix(0,
# ncol = n, nrow = n) Sigma = cov(rets)
# for (i in 1 : t){ Sigma = lambda * Sigma + (1 - lambda) * outer(
# as.double(rets[i, ]), as.double(rets[i, ]) ) } Finite sample ewma
t <- nrow(rets)
Sigma <- cov(rets)
# Sigma = matrix(0, ncol = dim(rets)[2], dim(rets)[2])
mu <- colMeans(rets)
shiftRets <- sweep(rets, 2, mu, "-")
for (i in 1:t) {
r <- as.double(shiftRets[i, ])
r2 <- outer(r, r)
Sigma <- (1 - lambda)/(1 - lambda^t) * r2 + lambda * Sigma
}
return(Sigma)
}
.constCov <- function(rets) {
## shrinkage of the covariance matrix using constant correlation matrix
## INPUTs rets : matrix (T x N) returns OUTPUTs Sigma : matrix (N x N)
## covariance matrix
n <- dim(rets)[2]
tmpMat <- matrix(rep(1, n^2), ncol = n)
rho <- mean(cor(rets)[lower.tri(tmpMat, diag = FALSE)])
R <- rho * tmpMat
diag(R) <- 1
std <- apply(rets, 2, sd)
diagStd <- diag(std)
Sigma <- diagStd %*% R %*% diagStd
return(Sigma)
}
.factorCov <- function(rets, K) {
## Compute the covariance matrix using K-factor approach INPUTs rets :
## matrix (T x N) returns K : [scalar] number of factors OUTPUTs Sigma :
## matrix (N x N) covariance matrix NOTE Matlab function is factoran
## (statistical toolbox) !!! TOFIX !!! valider avec Matlab
std <- apply(rets, 2, sd)
sigma <- cov(rets)
loading <- factanal(rets, K)$loadings
uniquenesses <- factanal(rets, K)$uniquenesses
R <- tcrossprod(loading) + diag(uniquenesses)
diagStd <- diag(std)
Sigma <- diagStd %*% R %*% diagStd
return(Sigma)
}
.lwCovElement <- function(rets, type) {
## Computes the common elements of the functions for Ledoit-Wolf INPUTs
## rets : matrix (T x N) returns type : 'large', 'lw' or else OUTPUTs
## list : useful outputs
t <- dim(rets)[1]
n <- dim(rets)[2]
mu <- colMeans(rets)
shiftRets <- sweep(rets, 2, mu, "-")
y <- shiftRets^2
if (type == "large" || type == "lw") {
# mkt = rowMeans(rets)
mkt <- rowMeans(shiftRets)
smple <- cov(cbind(rets, mkt)) * (t - 1)/t
covmkt <- smple[1:n, n + 1]
covmkt_ <- matrix(rep(covmkt, n), ncol = n, byrow = FALSE)
varmkt <- as.numeric(smple[n + 1, n + 1])
smple <- smple[-(n + 1), -(n + 1)]
prior <- outer(covmkt, covmkt)/varmkt
diag(prior) <- diag(smple)
lwCovElement <- list(t = t, n = n, mu = mu, shiftRets = shiftRets,
y = y, smple = smple, mkt = mkt, covmkt = covmkt, covmkt_ = covmkt_,
varmkt = varmkt, prior = prior)
} else {
smple <- (1/t) * crossprod(shiftRets)
lwCovElement <- list(t = t, n = n, mu = mu, shiftRets = shiftRets,
y = y, smple = smple, mkt = NULL, covmkt = NULL, covmkt_ = NULL,
varmkt = NULL, prior = NULL)
}
return(lwCovElement)
}
.lwCov <- function(rets) {
## Shrinkage of the covariance matrix towards market INPUTs rets :
## matrix (T x N) returns OUTPUTs Sigma : matrix (N x N) covariance
## matrix
## Adapted from covMarket.m by Olivier Ledoit and Michael Wolf (2014)
lwCovElement <- .lwCovElement(rets, type = "lw")
t <- lwCovElement$t
n <- lwCovElement$n
mu <- lwCovElement$mu
shiftRets <- lwCovElement$shiftRets
mkt <- lwCovElement$mkt
covmkt <- lwCovElement$covmkt
varmkt <- lwCovElement$varmkt
smple <- lwCovElement$smple
prior <- lwCovElement$prior
y <- lwCovElement$y
z <- sweep(shiftRets, 1, mkt, "*")
# Phi hat
phiMat <- crossprod(y)/t - 2 * crossprod(shiftRets) * smple/t + smple^2
phi <- sum(apply(phiMat, 2, sum))
# Rho hat
rhoMat1 <- 1/t * sweep(crossprod(y, z), 2, covmkt, "*")/varmkt
rhoMat3 <- 1/t * crossprod(z) * outer(covmkt, covmkt)/varmkt^2
rhoMat <- 2 * rhoMat1 - rhoMat3 - prior * smple
diag(rhoMat) <- diag(phiMat)
rho <- sum(apply(rhoMat, 2, sum))
# Gamma hat
gamma <- norm(smple - prior, "F")^2
# Kappa hat
kappa <- (phi - rho)/gamma
# shrinkage value
shrinkage <- pmax(0, pmin(1, kappa/t))
# Sigma hat
Sigma <- shrinkage * prior + (1 - shrinkage) * smple
return(Sigma)
}
.largeCov <- function(rets) {
## INPUTs rets : matrix (T x N) returns OUTPUTs Sigma : matrix (N x N)
## covariance matrix !!! TOFIX !!! largeCov ou puis-je trouver les
## sources pour formules ?
lwCovElement <- .lwCovElement(rets, type = "large")
t <- lwCovElement$t
n <- lwCovElement$n
mu <- lwCovElement$mu
shiftRets <- lwCovElement$shiftRets
y <- lwCovElement$y
mkt <- lwCovElement$mkt
covmkt <- lwCovElement$covmkt
covmkt_ <- lwCovElement$covmkt_
varmkt <- lwCovElement$varmkt
smple <- lwCovElement$smple
prior <- lwCovElement$prior
z <- sweep(shiftRets, 1, mkt, "*")
d <- 1/n * norm(smple - prior, "F")^2
r2 <- 1/n/t^2 * sum(apply(crossprod(y, y), 2, sum)) - 1/n/t * sum(apply(smple^2,
2, sum))
phidiag <- 1/n/t^2 * sum(apply(y^2, 2, sum)) - 1/n/t * sum(diag(smple)^2)
v1 <- 1/t^2 * crossprod(y, z) - 1/t * covmkt_ * smple
phioff1 <- 1/n * sum(apply(v1 * t(covmkt_), 2, sum))/varmkt - 1/n *
sum(diag(v1) * covmkt)/varmkt
v3 <- 1/t^2 * crossprod(z, z) - 1/t * varmkt * smple
phioff3 <- 1/n * sum(apply(v3 * tcrossprod(covmkt, covmkt), 2, sum))/varmkt^2 -
1/n * sum(diag(v3) * covmkt^2)/varmkt^2
phioff <- 2 * phioff1 - phioff3
phi <- phidiag + phioff
Sigma <- (r2 - phi)/d * prior + (1 - (r2 - phi)/d) * smple
return(Sigma)
}
.corCov <- function(rets) {
## Shrinkage of the covariance matrix towards constant correlation
## matrix by Ledoit-Wolf INPUTs rets : matrix (T x N) returns OUTPUTs
## Sigma : matrix (N x N) covariance matrix
## Adapted from covCor.m by Olivier Ledoit and Michael Wolf (2014)
lwCovElement <- .lwCovElement(rets, type = "cor")
t <- lwCovElement$t
n <- lwCovElement$n
mu <- lwCovElement$mu
shiftRets <- lwCovElement$shiftRets
y <- lwCovElement$y
smple <- lwCovElement$smple
var <- diag(smple)
sqrtvar <- sqrt(var)
outerSqrtVar <- outer(sqrtvar, sqrtvar)
rBar <- (sum(sum(smple/outerSqrtVar)) - n)/(n * (n - 1))
prior <- rBar * outerSqrtVar
diag(prior) <- var
# phi hat
phiMat <- crossprod(y)/t - 2 * crossprod(shiftRets) * smple/t + smple^2
phi <- sum(apply(phiMat, 2, sum))
# rho hat
term1 <- crossprod(shiftRets^3, shiftRets)/t
term2 <- sweep(smple, 1, diag(smple), "*")
# term3 = sweep(smple, 1, var, '*') # lorsqu'on developpe on pourrait
# annuler les term3 et term4 ... R term4 = sweep(smple, 1, var, '*')
rhoMat <- term1 - term2 #- term3 + term4
diag(rhoMat) <- 0
rho <- sum(diag(phiMat)) + rBar * sum(sum(outer(1/sqrtvar, sqrtvar) *
rhoMat))
# gamma hat
gamma <- norm(smple - prior, type = "F")^2
# kappa hat
kappa <- (phi - rho)/gamma
# shrinkage value
shrinkage <- pmax(0, pmin(1, kappa/t))
# Sigma hat
Sigma <- shrinkage * prior + (1 - shrinkage) * smple
return(Sigma)
}
.diagCov <- function(rets) {
## Shrinks towards diagonal matrix INPUTs rets : matrix (T x N) returns
## OUTPUTs Sigma : matrix (N x N) covariance matrix
lwCovElement <- .lwCovElement(rets, type = "diag")
t <- lwCovElement$t
n <- lwCovElement$n
mu <- lwCovElement$mu
shiftRets <- lwCovElement$shiftRets
y <- lwCovElement$y
smple <- lwCovElement$smple
prior <- diag(diag(smple))
# phi hat
phiMat <- crossprod(y)/t - 2 * (crossprod(shiftRets)) * smple/t + smple^2
phi <- sum(apply(phiMat, 2, sum))
# rho hat
rho <- sum(diag(phiMat))
# gamma hat
gamma <- norm(smple - prior, "F")^2
# kappa hat
kappa <- (phi - rho)/gamma
# shrinkage value
shrinkage <- pmax(0, pmin(1, kappa/t))
# Sigma hat
Sigma <- shrinkage * prior + (1 - shrinkage) * smple
return(Sigma)
}
.oneparmCov <- function(rets) {
## Shrinks towards one-parameter matrix INPUTs rets : matrix (T x N)
## returns OUTPUTs Sigma : matrix (N x N) covariance matrix
## Adapted from cov1para.m by Olivier Ledoit and Michael Wolf (2014)
lwCovElement <- .lwCovElement(rets, type = "oneparm")
t <- lwCovElement$t
n <- lwCovElement$n
mu <- lwCovElement$mu
shiftRets <- lwCovElement$shiftRets
smple <- lwCovElement$smple
y <- lwCovElement$y
meanvar <- mean(diag(smple))
prior <- meanvar * diag(n)
# phi hat
phiMat <- crossprod(y)/t - 2 * (crossprod(shiftRets)) * smple/t + smple^2
phi <- sum(apply(phiMat, 2, sum))
# gamma hat
gamma <- norm(smple - prior, type = "F")^2
# kappa hat
kappa <- phi/gamma
# shrinkage value
shrinkage <- pmax(0, pmin(1, kappa/t))
# Sigma hat
Sigma <- shrinkage * prior + (1 - shrinkage) * smple
return(Sigma)
}
.bsCov <- function(rets) {
## Compute the Bayes-Stein covariance matrix INPUTs rets : matrix (T x
## returns OUTPUTs Sigma : Matrix (N x N) covariance David : see
## Kolusheva, Daniela. (July 2008) Out-of-sample Performance of Asset
## Allocation Strategies
lwCovElement <- .lwCovElement(rets, type = "lw")
t <- lwCovElement$t
n <- lwCovElement$n
mu <- colMeans(rets)
Sigma <- cov(rets)
invSigma <- solve(Sigma)
i <- rep(1, n)
invSigmai <- crossprod(invSigma, i)
w_min <- (invSigmai)/as.numeric(crossprod(i, invSigmai))
mu_min <- crossprod(mu, w_min)
invSigmaMu <- crossprod(invSigma, mu - mu_min)
phi <- (n + 2)/((n + 2) + t * crossprod(mu - mu_min, invSigmaMu))
phi <- max(min(phi, 1), 0)
tau <- t * phi/(1 - phi)
Sigma <- Sigma * (1 + 1/(t + tau)) +
tau/(t * (t + 1 + tau)) * outer(i, i)/as.numeric(crossprod(i, invSigmai))
return(Sigma)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.