factors_r <- c("SP500", "DTWEXAFEGS") # "SP500" does not contain dividends; note: "DTWEXM" discontinued as of Jan 2020
factors_d <- c("DGS10", "BAMLH0A0HYM2")tickers <- c("BAICX") # fund inception date is "2011-11-28"intercept <- TRUERegression analysis
A factor model expresses fund returns as exposures to a set of market factors, which can be estimated with regression analysis. To measure the risk of a fund, estimate factor exposures for an allocation fund, decompose volatility into factor contributions, and then stress the factors under scenarios.
The example uses equity, currency, interest rate, and credit factors with a trailing one-year window of overlapping returns and equal weights (note: rates and spreads are negated yield changes, i.e. duration-style returns, and volatility is annualized using the square root of time).
Ordinary least squares
Coefficients
Solve for the weighted least squares coefficients in matrix form:
\[ \begin{aligned} \hat{\beta}=(X^\mathrm{T}WX)^{-1}X^\mathrm{T}Wy \end{aligned} \]
lm_coef <- function(x, y, weights, intercept) {
# cbind.xts() changes names, e.g., "(Intercept)" => "X.Intercept."
if (intercept) x <- cbind("(Intercept)" = 1, as.matrix(x))
result <- solve(crossprod(x, diag(weights)) %*% x) %*% crossprod(x, diag(weights) %*% y)
return(result)
}t(lm_coef(overlap_x_xts, overlap_y_xts, weights, intercept)) (Intercept) SP500 DTWEXAFEGS DGS10 BAMLH0A0HYM2
BAICX 0.0002087485 0.2236805 -0.1242992 3.283418 1.247491if (intercept) {
form <- reformulate(termlabels = factors, response = tickers)
} else {
form <- reformulate(termlabels = factors, response = tickers, intercept = FALSE)
}
fit <- lm(form, data = overlap_xts, weights = weights)
coef(fit) (Intercept) SP500 DTWEXAFEGS DGS10 BAMLH0A0HYM2
0.0002087485 0.2236804831 -0.1242991704 3.2834179825 1.2474909705 R-squared
Then compute the weighted R-squared, i.e. the proportion of variance explained by the factors:
\[ \begin{aligned} R^{2}=\frac{\hat{\beta}^\mathrm{T}(X^\mathrm{T}WX)\hat{\beta}}{y^\mathrm{T}Wy} \end{aligned} \]
lm_rsq <- function(x, y, weights, intercept) {
coef <- lm_coef(x, y, weights, intercept)
if (intercept) {
x <- cbind(1, x)
x <- sweep(x, 2, apply(x, 2, weighted.mean, w = weights), "-")
y <- sweep(y, 2, apply(y, 2, weighted.mean, w = weights), "-")
}
result <- (t(coef) %*% (crossprod(x, diag(weights)) %*% x) %*% coef) / (crossprod(y, diag(weights)) %*% y)
colnames(result) <- "R-squared"
return(result)
}lm_rsq(overlap_x_xts, overlap_y_xts, weights, intercept) R-squared
BAICX 0.7123818summary(fit)$r.squared[1] 0.7123818Standard errors
Also compute the standard errors of the intercept and coefficients using the residual variance:
\[ \begin{aligned} \sigma_{\hat{\beta}}^{2}&=\sigma_{\varepsilon}^{2}(X^\mathrm{T}WX)^{-1}\\ &=\frac{(1-R^{2})(y^\mathrm{T}Wy)}{n-p}(X^\mathrm{T}WX)^{-1}\\ &=\frac{SSE}{df_{E}}(X^\mathrm{T}WX)^{-1}\\ \sigma_{\hat{\alpha}}^{2}&=\sigma_{\varepsilon}^{2}\left(\frac{1}{\sum_{i=1}^{n}w_{i}}+\mu^\mathrm{T}(X^\mathrm{T}WX)^{-1}\mu\right) \end{aligned} \]
lm_se <- function(x, y, weights, intercept) {
n_rows <- nrow(x)
n_cols <- ncol(x)
rsq <- lm_rsq(x, y, weights, intercept)
if (intercept) {
# cbind.xts() changes names, e.g., "(Intercept)" => "X.Intercept."
x <- cbind("(Intercept)" = 1, as.matrix(x))
y <- sweep(y, 2, apply(y, 2, weighted.mean, w = weights), "-")
df_resid <- n_rows - n_cols - 1
} else {
df_resid <- n_rows - n_cols
}
var_y <- crossprod(y, diag(weights)) %*% y
var_resid <- as.numeric((1 - rsq) * var_y / df_resid)
result <- sqrt(var_resid * diag(solve(crossprod(x, diag(weights)) %*% x)))
return(result)
}lm_se(overlap_x_xts, overlap_y_xts, weights, intercept) (Intercept) SP500 DTWEXAFEGS DGS10 BAMLH0A0HYM2
5.719259e-05 2.156491e-02 4.231182e-02 3.444225e-01 3.359372e-01 coef(summary(fit))[ , "Std. Error"] (Intercept) SP500 DTWEXAFEGS DGS10 BAMLH0A0HYM2
5.719259e-05 2.156491e-02 4.231182e-02 3.444225e-01 3.359372e-01 Shapley values
Since the factors are correlated, allocate the R-squared to each factor using the Shapley value, i.e. the average marginal contribution across all orderings of factors:
\[ R^{2}_{i}=\sum_{S\subseteq N\setminus\{i\}}{\frac{|S|!\;(n-|S|-1)!}{n!}}(R^{2}(S\cup\{i\})-R^{2}(S)) \]
lm_shap <- function(x, y, weights, intercept) {
n_rows <- nrow(x)
n_cols <- ncol(x)
n_combn <- 2 ^ n_cols
n_vec <- array(0, n_combn)
ix_mat <- matrix(0, nrow = n_cols, ncol = n_combn)
rsq <- array(0, n_combn)
result <- array(0, n_cols)
# number of binary combinations
for (k in 1:n_combn) {
n <- 0
n_size <- k - 1
# find the binary combination
for (j in 1:n_cols) {
if (n_size %% 2 == 0) {
n <- n + 1
ix_mat[j, k] <- j
}
n_size <- n_size %/% 2
}
n_vec[k] <- n
if (n > 0) {
ix_subset <- which(ix_mat[ , k] != 0)
x_subset <- x[ , ix_subset]
rsq[k] <- lm_rsq(x_subset, y, weights, intercept)
}
}
# calculate the exact Shapley value for r-squared
for (j in 1:n_cols) {
ix_pos <- which(ix_mat[j, ] != 0)
ix_neg <- which(ix_mat[j, ] == 0)
ix_n <- n_vec[ix_neg]
rsq_diff <- rsq[ix_pos] - rsq[ix_neg]
for (k in 1:(n_combn / 2)) {
s <- ix_n[k]
weight <- factorial(s) * factorial(n_cols - s - 1) / factorial(n_cols)
result[j] <- result[j] + weight * rsq_diff[k]
}
}
names(result) <- colnames(x)
return(result)
}lm_shap(overlap_x_xts, overlap_y_xts, weights, intercept) SP500 DTWEXAFEGS DGS10 BAMLH0A0HYM2
0.32815100 0.09371267 0.15141857 0.13909952 Principal component regression
Alternatively, to address correlated factors, regress on the principal components of the factors and then map the coefficients to the original variables:
library(pls) # "Error in mvrValstats(object = fit, estimate = 'train'): could not find function 'mvrValstats'"comps <- 1Coefficients
\[ \begin{aligned} W_{k}&=\mathbf{X}V_{k}=[\mathbf{X}\mathbf{v}_{1},\ldots,\mathbf{X}\mathbf{v}_{k}]\\ {\widehat{\gamma}}_{k}&=\left(W_{k}^\mathrm{T}W_{k}\right)^{-1}W_{k}^\mathrm{T}\mathbf{Y}\\ {\widehat{\boldsymbol{\beta}}}_{k}&=V_{k}{\widehat{\gamma}}_{k} \end{aligned} \]
pcr_coef <- function(x, y, comps) {
x <- sweep(x, 2, colMeans(x), "-")
LV <- eigen(cov(x))
V <- LV[["vectors"]]
W <- x %*% V
gamma <- solve(crossprod(W)) %*% (crossprod(W, y))
result <- V[ , 1:comps] %*% as.matrix(gamma[1:comps])
return(result)
}scale_x_xts <- scale(overlap_x_xts)t(pcr_coef(scale_x_xts, overlap_y_xts, comps)) [,1] [,2] [,3] [,4]
[1,] 0.0005653212 -0.0004504335 0.0001678333 0.0005012898t(pcr_coef(overlap_x_xts, overlap_y_xts, comps)) [,1] [,2] [,3] [,4]
[1,] 0.3147428 -0.06562128 0.002804879 0.01264036fit <- pls::pcr(reformulate(termlabels = ".", response = tickers),
data = merge(scale_x_xts, overlap_y_xts), ncomp = comps)
coef(fit)[ , , 1] SP500 DTWEXAFEGS DGS10 BAMLH0A0HYM2
0.0005653212 -0.0004504335 0.0001678333 0.0005012898 R-squared
pcr_rsq <- function(x, y, comps) {
coef <- pcr_coef(x, y, comps)
x <- sweep(x, 2, colMeans(x), "-")
y <- sweep(y, 2, colMeans(y), "-")
result <- (t(coef) %*% crossprod(x) %*% coef) / crossprod(y)
colnames(result) <- "R-squared"
return(result)
}pcr_rsq(scale_x_xts, overlap_y_xts, comps) R-squared
BAICX 0.6190123pcr_rsq(overlap_x_xts, overlap_y_xts, comps) R-squared
BAICX 0.5958848pls::R2(fit)$val[comps + 1][1] 0.6190123Standard errors
\[ \begin{aligned} \text{Var}({\widehat{\boldsymbol{\beta}}}_{k})&=\sigma^{2}V_{k}(W_{k}^\mathrm{T}W_{k})^{-1}V_{k}^\mathrm{T}\\ &=\sigma^{2}V_{k}\text{diag}\left(\lambda_{1}^{-1},\ldots,\lambda_{k}^{-1}\right)V_{k}^\mathrm{T}\\ &=\sigma^{2}\sum_{j=1}^{k}{\frac{\mathbf{v}_{j}\mathbf{v}_{j}^\mathrm{T}}{\lambda_{j}}} \end{aligned} \]
Standard errors are equivalent to the standard errors from a regression on the principal component scores, mapped to the original coefficients, i.e. use the residual degrees of freedom for the \(k\) retained components:
pcr_se <- function(x, y, comps) {
n_rows <- nrow(x)
rsq <- pcr_rsq(x, y, comps)
y <- sweep(y, 2, colMeans(y), "-")
df_resid <- n_rows - comps - 1
var_y <- crossprod(y)
var_resid <- as.numeric((1 - rsq) * var_y / df_resid)
LV <- eigen(cov(x))
L <- LV$values[1:comps] * (n_rows - 1)
V <- LV$vectors[ , 1:comps]
result <- sqrt(var_resid * diag(V %*% sweep(t(V), 1, 1 / L, "*")))
return(result)
}pcr_se(scale_x_xts, overlap_y_xts, comps)[1] 2.804989e-05 2.234943e-05 8.327487e-06 2.487280e-05pcr_se(overlap_x_xts, overlap_y_xts, comps)[1] 0.0163929530 0.0034177958 0.0001460883 0.0006583562LV <- eigen(cov(scale_x_xts))
V <- as.matrix(LV[["vectors"]][ , 1:comps])
W <- scale_x_xts %*% V
fit <- lm(zoo::coredata(overlap_y_xts) ~ W)
sqrt(diag(V %*% as.matrix(vcov(fit)[-1, -1]) %*% t(V)))[1] 2.804989e-05 2.234943e-05 8.327487e-06 2.487280e-05Partial least squares
Risk decomposition
Next, decompose the volatility of the fund using the factor model:
Standalone risk
Standalone risk is the volatility of each factor exposure in isolation, i.e. ignoring correlations, so the components do not sum to the total:
\[ \begin{aligned} \text{SAR}_{k}&=\sqrt{w_{k}^{2}\sigma_{k}^{2}}\\ \text{SAR}_{\varepsilon}&=\sqrt{(1-R^{2})\sigma_{y}^{2}} \end{aligned} \]
lm_sar <- function(x, y, weights, intercept) {
coef <- lm_coef(x, y, weights, intercept)
rsq <- lm_rsq(x, y, weights, intercept)
if (intercept) x <- cbind(1, x)
sigma <- cov.wt(cbind(x, y), wt = weights, center = intercept)$cov
sar <- coef ^ 2 * diag(sigma[-ncol(sigma), -ncol(sigma)])
sar_eps <- (1 - rsq) * sigma[ncol(sigma), ncol(sigma)]
result <- sqrt(c(sigma[ncol(sigma), ncol(sigma)],
sar,
sar_eps))
names(result) <- c("total", rownames(coef), "eps")
return(result)
}lm_sar(overlap_x_xts, overlap_y_xts, weights, intercept) * sqrt(scale[["periods"]] * scale[["overlap"]]) total (Intercept) SP500 DTWEXAFEGS DGS10 BAMLH0A0HYM2
0.057469905 0.000000000 0.030238254 0.006591822 0.019855872 0.010113892
eps
0.030821139 Risk contribution
Marginal contribution to risk allocates the total volatility to each factor exposure, i.e. including correlations, so the components sum to the total:
\[ \begin{aligned} \text{MCR}_{k}&=w_{k}\frac{\partial\sigma_{y}}{\partial w_{k}}\\ &=w_{k}\frac{(\Sigma w)_{k}}{\sigma_{y}}\\ \text{MCR}_{\varepsilon}&=\sigma_{y}-\sum_{k=1}^{n}\text{MCR}_{k} \end{aligned} \]
lm_mcr <- function(x, y, weights, intercept) {
coef <- lm_coef(x, y, weights, intercept)
rsq <- lm_rsq(x, y, weights, intercept)
if (intercept) x <- cbind(1, x)
sigma <- cov.wt(cbind(x, y), wt = weights, center = intercept)$cov
mcr <- coef * sigma[-ncol(sigma), -ncol(sigma)] %*% coef / sqrt(sigma[ncol(sigma), ncol(sigma)])
mcr_eps <- sqrt(sigma[ncol(sigma), ncol(sigma)]) - sum(mcr)
result <- c(sqrt(sigma[ncol(sigma), ncol(sigma)]),
mcr,
mcr_eps)
names(result) <- c("total", rownames(coef), "eps")
return(result)
}lm_mcr(overlap_x_xts, overlap_y_xts, weights, intercept) * sqrt(scale[["periods"]] * scale[["overlap"]]) total (Intercept) SP500 DTWEXAFEGS DGS10 BAMLH0A0HYM2
0.057469905 0.000000000 0.023082945 0.003223079 0.009119969 0.005514519
eps
0.016529393 Scenario analysis
Finally, estimate the impact of scenarios on the fund using the factor exposures:
Implied shocks
Implied shocks translate a scenario for a subset of factors into consistent shocks for all factors using weighted least squares:
\[ \begin{aligned} \hat{\beta}&=(Z^\mathrm{T}WZ)^{-1}Z^\mathrm{T}WX \end{aligned} \]
implied_shocks <- function(shocks, x, z, weights) {
beta <- solve(crossprod(z, diag(weights) %*% z)) %*% crossprod(z, diag(weights) %*% x)
result <- shocks %*% beta
return(result)
}shocks <- c(-0.1, 0.1)
overlap_z_xts <- overlap_x_xts[ , 1:2]implied_shocks(shocks, overlap_x_xts, overlap_z_xts, weights) SP500 DTWEXAFEGS DGS10 BAMLH0A0HYM2
[1,] -0.1 0.1 -0.002248735 -0.004155653Stress P&L
Stress P&L applies the implied shocks to the factor exposures. Note: if an intercept is included, then alpha is stressed along with the factors, i.e. the constant is passed through the implied shocks regression:
pnl_stress <- function(shocks, x, y, z, weights, intercept) {
coef <- lm_coef(x, y, weights, intercept)
# cbind.xts() changes names, e.g., "(Intercept)" => "X.Intercept."
if (intercept) x <- cbind("(Intercept)" = 1, as.matrix(x))
result <- t(coef) * implied_shocks(shocks, x, z, weights)
return(result)
}pnl_stress(shocks, overlap_x_xts, overlap_y_xts, overlap_z_xts, weights, intercept) (Intercept) SP500 DTWEXAFEGS DGS10 BAMLH0A0HYM2
BAICX 0.0004116162 -0.02236805 -0.01242992 -0.007383537 -0.00518414