# basic finance
library(xts) # to manipulate time series of stock data
library(pob) # book package with financial data
# install with: devtools::install_github("dppalomar/pob")
# data
library(mvtnorm) # to generate heavy-tailed data
library(fitHeavyTail) # to fit heavy-tailed distributions
library(ICSNP) # to calculate spatial mean
library(covFactorModel) # factor model estimation
# devtools::install_github("dppalomar/covFactorModel")
# plotting
library(ggplot2) # for nice plots
library(patchwork) # for combining plots
library(reshape2) # to reshape data
library(dplyr); conflictRules('dplyr', exclude = 'lag'); options(xts.warn_dplyr_breaks_lag = FALSE)
library(scales)
library(latex2exp) # for latex symbols with TeX()Portfolio Optimization
Chapter 3 - Financial Data: IID Modeling
R code
R code examples for Chapter 3 of the book:
Daniel P. Palomar (2024). Portfolio Optimization: Theory and Application. Cambridge University Press.
Loading packages
The following packages are used in the examples:
IID model
Example of a synthetic Gaussian i.i.d. time series:
library(mvtnorm)
# get realistic mu and variance
sp500_prices <- SP500_2015to2020$index
sp500_returns <- diff(log(sp500_prices))[-1]
mu <- mean(sp500_returns)
var <- var(sp500_returns)
# generate synthetic data
set.seed(42)
T <- 500
x <- rmvnorm(n = T, mean = mu, sigma = var)
data.frame(t = 1:T, x = x) |>
ggplot(aes(x = t, y = x)) +
geom_line(linewidth = 0.8, color = "blue", show.legend = FALSE) +
labs(title = "i.i.d. time series", x = NULL, y = NULL)
Sample estimators
Estimation error of sample estimators versus number of observations (for Gaussian data with \(N=100\)):
library(mvtnorm)
# get realistic mu and Sigma
X <- diff(log(SP500_2015to2020$stocks))[-1, 1:100]
mu_true <- colMeans(X)
Sigma_true <- cov(X)
# generate synthetic data
set.seed(42)
T_max <- 2000
X <- rmvnorm(n = 10*T_max, mean = mu_true, sigma = Sigma_true)
# main loop
df <- data.frame()
T_sweep <- seq(from = 50, by = 50, to = T_max)
for(T in T_sweep) {
for (i in 1:100) {
X_ <- X[sample(nrow(X), T), ]
mu <- colMeans(X_)
Sigma <- cov(X_)
df <- rbind(df,
data.frame("T" = T,
"method" = "sample mean",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
"error Sigma" = norm(Sigma - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
}
}
df <- df |>
group_by(method, T) |>
summarize("error mu" = 100*mean(`error mu`),
"error Sigma" = 100*mean(`error Sigma`)) |>
ungroup()
# plot
p_error_mu <- df |>
ggplot(aes(x = T, y = `error mu`, color = method)) +
geom_line(linewidth = 1.2, show.legend = FALSE) +
labs(title = "Estimation error of sample mean", x = "T", y = "normalized error (%)")
p_error_Sigma <- df |>
ggplot(aes(x = T, y = `error Sigma`, color = method)) +
geom_line(linewidth = 1.2, show.legend = FALSE) +
labs(title = "Estimation error of sample covariance matrix", x = "T", y = "normalized error (%)")
p_error_mu / p_error_Sigma
Location estimators
Estimation error of location estimators versus number of observations (with \(N=100\)):
library(mvtnorm)
# get realistic mu and Sigma
X <- diff(log(SP500_2015to2020$stocks))[-1, 1:100]
mu_true <- colMeans(X)
Sigma_true <- cov(X)
nu <- 4
scatter_true <- (nu-2)/nu * Sigma_true # cov=nu/(nu-2)*scatter
# generate synthetic data
set.seed(42)
T_max <- 1000
X_Gaussian <- rmvnorm(n = 10*T_max, mean = mu_true, sigma = Sigma_true)
X_heavy_tail <- rmvt(n = 10*T_max, delta = mu_true, sigma = scatter_true, df = nu)
# main loop
df <- data.frame()
T_sweep <- seq(from = 50, by = 50, to = T_max)
for(T in T_sweep) {
for (distribution in c("Gaussian", "heavy-tailed")) {
for (i in 1:100) {
if (distribution == "Gaussian")
X_ <- X_Gaussian[sample(nrow(X_Gaussian), T), ]
else if (distribution == "heavy-tailed")
X_ <- X_heavy_tail[sample(nrow(X_heavy_tail), T), ]
else stop("Unknown distribution.")
# sample mean
mu <- colMeans(X_)
df <- rbind(df, data.frame("T" = T,
"distribution" = distribution,
"method" = "sample mean",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
check.names = FALSE))
# median
mu <- apply(X_, 2, median)
df <- rbind(df, data.frame("T" = T,
"distribution" = distribution,
"method" = "median",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
check.names = FALSE))
# spatial median
mu <- ICSNP::spatial.median(X_)
df <- rbind(df, data.frame("T" = T,
"distribution" = distribution,
"method" = "spatial median",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
check.names = FALSE))
}
}
}
df$method <- factor(df$method, levels = unique(df$method))
df <- df |>
group_by(method, T, distribution) |>
summarize("error mu" = 100*mean(`error mu`)) |>
ungroup()
# plot0
p_Gaussian <- df[df$distribution == "Gaussian", ] |>
ggplot(aes(x = T, y = `error mu`, color = method,)) +
geom_line(linewidth = 1.2) +
labs(color = "estimator") +
coord_cartesian(xlim = c(1, T_max)) +
labs(title = "Error of location estimators under Gaussian data", x = "T", y = "normalized error (%)")
p_heavy_tail <- df[df$distribution == "heavy-tailed", ] |>
ggplot(aes(x = T, y = `error mu`, color = method)) +
geom_line(linewidth = 1.2) +
labs(color = "estimator") +
coord_cartesian(xlim = c(1, T_max)) +
labs(title = "Error of location estimators under heavy-tailed data", x = "T", y = "normalized error (%)")
p_Gaussian / p_heavy_tail + plot_layout(guides = "collect")
Estimation error of location estimators versus degrees of freedom in a \(t\) distribution (with \(T=200\) and \(N=100\)):
library(mvtnorm)
# get realistic mu and Sigma
X <- diff(log(SP500_2015to2020$stocks))[-1, 1:100]
mu_true <- colMeans(X)
Sigma_true <- cov(X)
T <- 200
# main loop
set.seed(42)
df <- data.frame()
nu_sweep <- seq(from = 3, by = 2, to = 30)
for(nu in nu_sweep) {
for (i in 1:200) {
# generate synthetic data
X_ <- rmvt(n = T, delta = mu_true, sigma = (nu-2)/nu * Sigma_true, df = nu)
# sample mean
mu <- colMeans(X_)
df <- rbind(df, data.frame("T" = T,
"nu" = nu,
"method" = "sample mean",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
check.names = FALSE))
# median
mu <- apply(X_, 2, median)
df <- rbind(df, data.frame("T" = T,
"nu" = nu,
"method" = "median",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
check.names = FALSE))
# spatial median
mu <- ICSNP::spatial.median(X_)
df <- rbind(df, data.frame("T" = T,
"nu" = nu,
"method" = "spatial median",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
check.names = FALSE))
}
}
df$method <- factor(df$method, levels = unique(df$method))
df <- df |>
group_by(method, T, nu) |>
summarize("error mu" = 100*mean(`error mu`)) |>
ungroup()
# plot
df |>
ggplot(aes(x = nu, y = `error mu`, color = method)) +
geom_line(linewidth = 1.2) +
labs(color = "estimator") +
labs(title = "Error of location estimators", x = TeX("$\\nu$"), y = "normalized error (%)")
Gaussian ML estimators
Estimation error of Gaussian ML estimators versus number of observations (with \(N=100\)):
library(mvtnorm)
# get realistic mu and Sigma
X <- diff(log(SP500_2015to2020$stocks))[-1, 1:100]
mu_true <- colMeans(X)
Sigma_true <- cov(X)
nu <- 4
scatter_true <- (nu-2)/nu * Sigma_true # cov=nu/(nu-2)*scatter
# generate synthetic data
set.seed(42)
T_max <- 1000
X_Gaussian <- rmvnorm(n = 10*T_max, mean = mu_true, sigma = Sigma_true)
X_heavy_tail <- rmvt(n = 10*T_max, delta = mu_true, sigma = scatter_true, df = nu)
# main loop
df <- data.frame()
T_sweep <- seq(from = 50, by = 50, to = T_max)
for(T in T_sweep) {
for (distribution in c("Gaussian", "heavy-tailed")) {
for (i in 1:100) {
if (distribution == "Gaussian")
X_ <- X_Gaussian[sample(nrow(X_Gaussian), T), ]
else if (distribution == "heavy-tailed")
X_ <- X_heavy_tail[sample(nrow(X_heavy_tail), T), ]
else stop("Unknown distribution.")
# sample mean
mu <- colMeans(X_)
Sigma <- (T-1)/T * cov(X_)
df <- rbind(df, data.frame("T" = T,
"distribution" = distribution,
"method" = "Gaussian ML",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
"error Sigma" = norm(Sigma - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
}
}
}
df$method <- factor(df$method, levels = unique(df$method))
df <- df |>
group_by(method, T, distribution) |>
summarize("error mu" = 100*mean(`error mu`),
"error Sigma" = 100*mean(`error Sigma`)) |>
ungroup()
# plot
p_error_mu <- df |>
ggplot(aes(x = T, y = `error mu`, color = distribution)) +
geom_line(linewidth = 1.2) +
labs(color = "data distribution") +
labs(title = "Error of Gaussian ML mean estimator", x = "T", y = "normalized error (%)")
p_error_Sigma <- df |>
ggplot(aes(x = T, y = `error Sigma`, color = distribution)) +
geom_line(linewidth = 1.2) +
labs(color = "data distribution") +
labs(title = "Error of Gaussian ML covariance estimator", x = "T", y = "normalized error (%)")
p_error_mu / p_error_Sigma + plot_layout(guides = "collect")
Estimation error of Gaussian ML estimators versus degrees of freedom in a \(t\) distribution (with \(T=200\) and \(N=100\)):
library(mvtnorm)
# get realistic mu and Sigma
X <- diff(log(SP500_2015to2020$stocks))[-1, 1:100]
mu_true <- colMeans(X)
Sigma_true <- cov(X)
T <- 200
# main loop
set.seed(42)
df <- data.frame()
nu_sweep <- seq(from = 3, by = 2, to = 30)
for(nu in nu_sweep) {
for (i in 1:200) {
# generate synthetic data
X <- rmvt(n = T, delta = mu_true, sigma = (nu-2)/nu * Sigma_true, df = nu)
# sample mean
mu <- colMeans(X)
Sigma <- (T-1)/T * cov(X)
df <- rbind(df, data.frame("T" = T,
"nu" = nu,
"method" = "Gaussian ML",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
"error Sigma" = norm(Sigma - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
}
}
df$method <- factor(df$method, levels = unique(df$method))
df <- df |>
group_by(method, T, nu) |>
summarize("error mu" = 100*mean(`error mu`),
"error Sigma" = 100*mean(`error Sigma`)) |>
ungroup()
# plot
p_error_mu <- df |>
ggplot(aes(x = nu, y = `error mu`, color = method)) +
geom_line(linewidth = 1.2, show.legend = FALSE) +
coord_cartesian(ylim = c(200, 250)) +
labs(title = "Estimation error of Gaussian ML mean estimator", x = TeX("$\\nu$"), y = "normalized error (%)")
p_error_Sigma <- df |>
ggplot(aes(x = nu, y = `error Sigma`, color = method)) +
geom_line(linewidth = 1.2, show.legend = FALSE) +
labs(title = "Estimation error of Gaussian ML covariance estimator", x = TeX("$\\nu$"), y = "normalized error (%)")
p_error_mu / p_error_Sigma
Heavy-tailed ML estimators
Estimation error of different ML estimators versus number of observations (for \(t\)-distributed heavy-tailed data with \(\nu=4\) and \(N=100\)):
library(mvtnorm)
library(fitHeavyTail)
library(ICSNP)
# get realistic mu and Sigma
X <- diff(log(SP500_2015to2020$stocks))[-1, 1:100]
mu_true <- colMeans(X)
Sigma_true <- cov(X)
nu <- 4
# generate synthetic data
set.seed(42)
T_max <- 1000
X <- rmvt(n = 10*T_max, delta = mu_true, sigma = (nu-2)/nu * Sigma_true, df = nu)
# main loop
df <- data.frame()
T_sweep <- seq(from = 150, by = 50, to = T_max)
for(T in T_sweep) {
for (i in 1:100) {
X_ <- X[sample(nrow(X), T), ]
# Gaussian MLE
mu <- colMeans(X_)
Sigma <- (T-1)/T * cov(X_)
df <- rbind(df, data.frame("T" = T,
"method" = "Gaussian MLE",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
"MAE mu" = sum(abs(mu - mu_true))/sum(abs(mu_true)),
"error Sigma" = norm(Sigma - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
# heavy-tailed MLE
fitted <- fit_mvt(X_, nu = "iterative", nu_iterative_method = "POP")
df <- rbind(df, data.frame("T" = T,
"method" = "heavy-tailed MLE",
"error mu" = norm(fitted$mu - mu_true, "2")/norm(mu_true, "2"),
"MAE mu" = sum(abs(fitted$mu - mu_true))/sum(abs(mu_true)),
"error Sigma" = norm(fitted$cov - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
# Tyler (with spatial median)
fitted <- fit_Tyler(X_, initial = list(mu = ICSNP::spatial.median(X_)))
df <- rbind(df, data.frame("T" = T,
"method" = "Tyler estimator",
"error mu" = norm(fitted$mu - mu_true, "2")/norm(mu_true, "2"),
"MAE mu" = sum(abs(fitted$mu - mu_true))/sum(abs(mu_true)),
"error Sigma" = norm(fitted$cov - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
}
}
df$method <- factor(df$method, levels = unique(df$method))
df <- df |>
group_by(method, T) |>
summarize("error mu" = 100*mean(`error mu`),
"RMSE mu" = 100*sqrt(mean(`error mu`)),
"MAE mu" = 100*mean(`MAE mu`),
"error Sigma" = 100*mean(`error Sigma`)) |>
ungroup()
# plot
p_error_mu <- df |>
ggplot(aes(x = T, y = `error mu`, color = method)) +
geom_line(linewidth = 1.2) +
scale_color_discrete(name = "estimator", labels = c("Gaussian MLE", "heavy-tailed MLE", "spatial median")) +
labs(title = "Error of mean estimators", x = "T", y = "normalized error (%)")
p_error_Sigma <- df |>
ggplot(aes(x = T, y = `error Sigma`, color = method)) +
geom_line(linewidth = 1.2) +
labs(color = "estimator") +
labs(title = "Error of covariance estimators", x = "T", y = "normalized error (%)")
p_error_mu / p_error_Sigma
Estimation error of different ML estimators versus degrees of freedom in a \(t\) distribution (with \(T=200\) and \(N=100\)):
library(mvtnorm)
library(fitHeavyTail)
library(ICSNP)
# get realistic mu and Sigma
X <- diff(log(SP500_2015to2020$stocks))[-1, 1:100]
mu_true <- colMeans(X)
Sigma_true <- cov(X)
T <- 200
# main loop
set.seed(42)
df <- data.frame()
nu_sweep <- seq(from = 3, by = 2, to = 30)
for(nu in nu_sweep) {
for (i in 1:200) {
# generate synthetic data
X <- rmvt(n = T, delta = mu_true, sigma = (nu-2)/nu * Sigma_true, df = nu)
# Gaussian MLE
mu <- colMeans(X)
Sigma <- (T-1)/T * cov(X)
df <- rbind(df, data.frame("T" = T,
"nu" = nu,
"method" = "Gaussian MLE",
"error mu" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
"error Sigma" = norm(Sigma - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
# heavy-tailed MLE
fitted <- fit_mvt(X, nu = "iterative", nu_iterative_method = "POP")
df <- rbind(df, data.frame("T" = T,
"nu" = nu,
"method" = "heavy-tailed MLE",
"error mu" = norm(fitted$mu - mu_true, "2")/norm(mu_true, "2"),
"error Sigma" = norm(fitted$cov - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
# Tyler (with spatial median)
fitted <- fit_Tyler(X, initial = list(mu = ICSNP::spatial.median(X)))
df <- rbind(df, data.frame("T" = T,
"nu" = nu,
"method" = "Tyler estimator",
"error mu" = norm(fitted$mu - mu_true, "2")/norm(mu_true, "2"),
"error Sigma" = norm(fitted$cov - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
}
}
df$method <- factor(df$method, levels = unique(df$method))
df <- df |>
group_by(method, T, nu) |>
summarize("error mu" = 100*mean(`error mu`),
"error Sigma" = 100*mean(`error Sigma`)) |>
ungroup()
# plot
p_error_mu <- df |>
ggplot(aes(x = nu, y = `error mu`, color = method)) +
geom_line(linewidth = 1.2) +
scale_color_discrete(name = "estimator", labels = c("Gaussian MLE", "heavy-tailed MLE", "spatial median")) +
labs(title = "Error of mean estimators", x = TeX("$\\nu$"), y = "normalized error (%)")
p_error_Sigma <- df |>
ggplot(aes(x = nu, y = `error Sigma`, color = method)) +
geom_line(linewidth = 1.2) +
labs(color = "estimator") +
labs(title = "Error of covariance estimators", x = TeX("$\\nu$"), y = "normalized error (%)")
p_error_mu / p_error_Sigma
Convergence of iterative robust heavy-tailed ML estimators:
library(mvtnorm)
library(fitHeavyTail)
# get realistic mu and Sigma
X <- diff(log(SP500_2015to2020$stocks))[-1, 1:100]
mu_true <- colMeans(X)
Sigma_true <- cov(X)
T <- 200
nu <- 4
# generate synthetic data
set.seed(42)
X <- rmvt(n = T, delta = mu_true, sigma = (nu-2)/nu * Sigma_true, df = nu)
# heavy-tailed MLE
fitted<- fit_mvt(X, nu = "iterative", nu_iterative_method = "ECM", ftol = 1e-3, return_iterates = TRUE)
# plot
data.frame("k" = 0:fitted$num_iterations,
"negative log-likelihood" = -fitted$log_likelihood_vs_iterations,
check.names = FALSE) |>
ggplot(aes(x = k, y = `negative log-likelihood`)) +
geom_line(linewidth = 1.2, color = "blue") +
ggtitle("Converge of objective function") + xlab("iterations")
Prior information
Shrinkage
First, we define a simple and convenient function to shrink the mean vector:
shrink_mu <- function(X, mu = NULL, cov = NULL, mu_sh_use_cor = TRUE,
mu_target = c("grand mean", "volatility weighted", "covmat weighted")) {
T <- nrow(X)
N <- ncol(X)
if (is.null(mu))
mu <- colMeans(X)
if (is.null(cov))
cov <- cov(X)
# target
if (!is.numeric(mu_target))
mu_target <- switch(match.arg(mu_target),
"grand mean" = mean(mu),
"volatility weighted" = {
vol <- sqrt(diag(cov))
mean(mu/vol) * vol
},
"covmat weighted" = {
invSigma_ones <- if (mu_sh_use_cor) solve(cov2cor(cov), rep(1, N))
else solve(cov, rep(1, N))
invSigma_ones <- pmax(0, invSigma_ones)
as.numeric(mu %*% invSigma_ones)/sum(invSigma_ones) # return achieved by GMVP
},
stop("Shrinkage estimation method for mu not recognized!"))
# shrinkage
rho <- (N + 2) / ((N + 2) + T * as.numeric((mu - mu_target) %*% solve(cov, mu - mu_target))) # improved James-Stein estimator
#rho <- (N - 2)/(T*as.numeric((mu - mu_target) %*% solve(cov, mu - mu_target))) # original James-Stein estimator
rho <- max(0, min(1, rho))
mu_shrinked <- (1 - rho) * mu + rho * mu_target
return(mu_shrinked)
}Estimation error of different shrinkage estimators versus number of observations (for Gaussian data with \(N=100\)):
library(mvtnorm)
# get realistic mu and Sigma
X <- diff(log(SP500_2015to2020$stocks))[-1, 1:100]
N <- ncol(X)
mu_true <- colMeans(X)
Sigma_true <- cov(X)
nu <- 8
# generate synthetic data
set.seed(42)
T_max <- 550
X <- rmvt(n = 10*T_max, delta = mu_true, sigma = (nu-2)/nu * Sigma_true, df = nu)
# main loop
df <- data.frame()
T_sweep <- seq(from = 50, by = 50, to = T_max)
for(T in T_sweep) {
for (i in 1:200) {
X_ <- X[sample(nrow(X), T), ]
############# #############
############# Mean estimation #############
############# #############
Sigma <- cov(X_)
Sigma_reg <- Sigma + mean(diag(Sigma)) * diag(N) # regularize Sigma for stability
#
# plain sample mean
#
mu <- colMeans(X_)
df <- rbind(df, data.frame("T" = T,
"parameter" = "mu",
"method" = "sample mean",
"error" = norm(mu - mu_true, "2")/norm(mu_true, "2"),
check.names = FALSE))
#
# sample mean + shrinkage to zero
#
mu_JS <- shrink_mu(X_, mu = colMeans(X_), mu_target = 0, cov = Sigma_reg)
df <- rbind(df, data.frame("T" = T,
"parameter" = "mu",
"method" = "shrinkage to zero",
"error" = norm(mu_JS - mu_true, "2")/norm(mu_true, "2"),
check.names = FALSE))
#
# sample mean + shrinkage to grand mean
#
mu_JS <- shrink_mu(X_, mu = colMeans(X_), mu_target = "grand mean", cov = Sigma_reg)
df <- rbind(df, data.frame("T" = T,
"parameter" = "mu",
"method" = "shrinkage to grand mean",
"error" = norm(mu_JS - mu_true, "2")/norm(mu_true, "2"),
check.names = FALSE))
#
# sample mean + shrinkage to volatility weighted mean
#
mu_JS <- shrink_mu(X_, mu = colMeans(X_), mu_target = "covmat weighted", cov = Sigma_reg, mu_sh_use_cor = FALSE)
df <- rbind(df, data.frame("T" = T,
"parameter" = "mu",
"method" = "shrinkage to volatility-weighted mean",
"error" = norm(mu_JS - mu_true, "2")/norm(mu_true, "2"),
check.names = FALSE))
############# #############
############# Covariance matrix estimation #############
############# #############
#
# plain SCM
#
Sigma <- cov(X_)
df <- rbind(df, data.frame("T" = T,
"parameter" = "Sigma",
"method" = "sample covariance matrix",
"error" = norm(Sigma - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
#
# SCM + Ledoit-Wolf shrinked to scaled identity
#
S <- cov(X_)
Sigma_target <- mean(diag(S)) * diag(N)
Xc <- X_ - matrix(colMeans(X_), T, N, byrow = TRUE)
rho <- min(1, ((1/T^2) * sum(rowSums(Xc^2)^2) - 1/T * sum(S^2)) / sum((S - Sigma_target)^2))
Sigma_sh <- (1 - rho)*S + rho*Sigma_target
df <- rbind(df, data.frame("T" = T,
"parameter" = "Sigma",
"method" = "shrinkage to scaled identity",
"error" = norm(Sigma_sh - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
#
# SCM + Ledoit-Wolf shrinked to diagonal matrix
#
S <- cov(X_)
Sigma_target <- diag(diag(S))
Xc <- X_ - matrix(colMeans(X_), T, N, byrow = TRUE)
rho <- min(1, ((1/T^2) * sum(rowSums(Xc^2)^2) - 1/T * sum(S^2)) / sum((S - Sigma_target)^2))
Sigma_sh <- (1 - rho)*S + rho*Sigma_target
df <- rbind(df, data.frame("T" = T,
"parameter" = "Sigma",
"method" = "shrinkage to diagonal matrix",
"error" = norm(Sigma_sh - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
#
# SCM + Ledoit-Wolf shrinked to equal-correlation matrix
#
S <- cov(X_)
sigma <- sqrt(diag(S))
C <- cov2cor(S)
ave_corr <- mean(C[lower.tri(C)])
C[lower.tri(C)] <- ave_corr
C[upper.tri(C)] <- ave_corr
Sigma_target <- t(C * sigma) * sigma
Xc <- X_ - matrix(colMeans(X_), T, N, byrow = TRUE)
rho <- min(1, ((1/T^2) * sum(rowSums(Xc^2)^2) - 1/T * sum(S^2)) / sum((S - Sigma_target)^2))
Sigma_sh <- (1 - rho)*S + rho*Sigma_target
df <- rbind(df, data.frame("T" = T,
"parameter" = "Sigma",
"method" = "shrinkage to equal-correlation matrix",
"error" = norm(Sigma_sh - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
}
}
df$method <- factor(df$method, levels = unique(df$method))
df <- df |>
group_by(parameter, method, T) |>
summarize("error" = 100*mean(`error`)) |>
ungroup()
# plot
p_error_mu <- df[df$parameter == "mu", ] |>
ggplot(aes(x = T, y = `error`, color = method)) +
geom_line(linewidth = 1.2) +
labs(color = "estimator") +
scale_x_continuous(breaks = seq(from = T_sweep[1], to = last(T_sweep), by = 100)) +
labs(title = "Error of mean estimators", x= "T", y = "normalized error (%)")
p_error_Sigma <- df[df$parameter == "Sigma", ] |>
ggplot(aes(x = T, y = `error`, color = method)) +
geom_line(linewidth = 1.2) +
labs(color = "estimator") +
scale_x_continuous(breaks = seq(from = T_sweep[1], to = last(T_sweep), by = 100)) +
labs(title = "Error of covariance estimators", x= "T", y = "normalized error (%)")
p_error_mu / p_error_Sigma
Factor models
Estimation error of covariance matrix under factor modeling versus number of observations (with \(N=100\)):
library(mvtnorm)
library(covFactorModel)
# get realistic mu and Sigma
X <- diff(log(SP500_2015to2020$stocks))[-1, 1:100]
N <- ncol(X)
mu_true <- rep(0, N) #colMeans(X)
# Forcing low-rank covmat - method #1
Sigma_true <- covFactorModel(X, type = "Stat-PCA", K = 3)
# generate synthetic data
set.seed(42)
T_max <- 500
X <- rmvt(n = 10*T_max, delta = mu_true, sigma = Sigma_true, df = Inf)
X <- as.xts(X, order.by = as.Date("1995-03-15") + 1:nrow(X))
# main loop
df <- data.frame()
T_sweep <- seq(from = 100, by = 50, to = T_max)
eigenvalues_true <- eigen(Sigma_true)$values
for(T in T_sweep) {
for (i in 1:200) {
X_ <- X[sample(nrow(X), T), ]
# SCM
Sigma <- cov(X_)
df <- rbind(df, data.frame("T" = T,
"method" = "sample covariance matrix",
"error eigenvalues" = norm(eigen(Sigma)$values - eigenvalues_true, "2")/norm(eigenvalues_true, "2"),
"error Sigma" = norm(Sigma - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
# 1-factor model
Sigma <- covFactorModel(X_, type = "Stat-PCA", K = 1)
df <- rbind(df, data.frame("T" = T,
"method" = "factor model (K=1)",
"error eigenvalues" = norm(eigen(Sigma)$values - eigenvalues_true, "2")/norm(eigenvalues_true, "2"),
"error Sigma" = norm(Sigma - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
# 3-factor model
Sigma <- covFactorModel(X_, type = "Stat-PCA", K = 3)
df <- rbind(df, data.frame("T" = T,
"method" = "factor model (K=3)",
"error eigenvalues" = norm(eigen(Sigma)$values - eigenvalues_true, "2")/norm(eigenvalues_true, "2"),
"error Sigma" = norm(Sigma - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
# 5-factor model
Sigma <- covFactorModel(X_, type = "Stat-PCA", K = 5)
df <- rbind(df, data.frame("T" = T,
"method" = "factor model (K=5)",
"error eigenvalues" = norm(eigen(Sigma)$values - eigenvalues_true, "2")/norm(eigenvalues_true, "2"),
"error Sigma" = norm(Sigma - Sigma_true, "F")/norm(Sigma_true, "F"),
check.names = FALSE))
}
}
df$method <- factor(df$method, levels = unique(df$method))
df <- df |>
group_by(method, T) |>
summarize("error Sigma" = 100*mean(`error Sigma`),
"error eigenvalues" = 100*mean(`error eigenvalues`)) |>
ungroup()`summarise()` has grouped output by 'method'. You can override using the
`.groups` argument.# plot
df |>
ggplot(aes(x = T, y = `error Sigma`, color = method)) +
geom_line(linewidth = 1.2) +
labs(color = "estimator") +
labs(title = "Error of covariance estimators", x= "T", y = "normalized error (%)")