Risk

analysis
finance
r
Author

Jason Foster

Published

April 6, 2025

factors_r <- c("SP500", "DTWEXAFEGS") # "SP500" does not contain dividends; note: "DTWEXM" discontinued as of Jan 2020
factors_d <- c("DGS10", "BAMLH0A0HYM2")
tickers <- "BAICX" # fund inception date is "2011-11-28" 
intercept <- TRUE

Regression analysis

Ordinary least squares

Coefficients

\[ \begin{aligned} \hat{\beta}=(X^\mathrm{T}WX)^{-1}X^\mathrm{T}Wy \end{aligned} \]

lm_coef <- function(x, y, weights, intercept) {
    
    if (intercept) x <- model.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)    xSP500 xDTWEXAFEGS   xDGS10 xBAMLH0A0HYM2
BAICX 7.30012e-05 0.1913491   -0.145976 3.230495      1.220434
if (intercept) {
    fit <- lm(overlap_y_xts ~ overlap_x_xts, weights = weights)
} else {
    fit <- lm(overlap_y_xts ~ overlap_x_xts - 1, weights = weights)
}
    
coef(fit)
              (Intercept)        overlap_x_xtsSP500   overlap_x_xtsDTWEXAFEGS 
             0.0000730012              0.1913490751             -0.1459759507 
       overlap_x_xtsDGS10 overlap_x_xtsBAMLH0A0HYM2 
             3.2304951748              1.2204339943

R-squared

\[ \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 <- model.matrix(~ 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)
    
    return(result)
    
}
lm_rsq(overlap_x_xts, overlap_y_xts, weights, intercept)
          BAICX
BAICX 0.7576384
summary(fit)$r.squared
[1] 0.7576384

Standard errors

\[ \begin{aligned} \sigma_{\hat{\beta}}^{2}&=\sigma_{\varepsilon}^{2}(X^\mathrm{T}WX)^{-1}\\ &=\frac{(1-R^{2})}{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}{n}+\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) {
        
        x <- model.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)        xSP500   xDTWEXAFEGS        xDGS10 xBAMLH0A0HYM2 
 5.121335e-05  1.940601e-02  4.495683e-02  2.662767e-01  2.582900e-01 
coef(summary(fit))[ , "Std. Error"]
              (Intercept)        overlap_x_xtsSP500   overlap_x_xtsDTWEXAFEGS 
             5.121335e-05              1.940601e-02              4.495683e-02 
       overlap_x_xtsDGS10 overlap_x_xtsBAMLH0A0HYM2 
             2.662767e-01              2.582900e-01

Shapley values

\[ 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 - 1 + 1
        
      }
      
      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]

    }
    
  }

  return(result)
  
}
lm_shap(overlap_x_xts, overlap_y_xts, weights, intercept)
[1] 0.26603748 0.13080190 0.26329100 0.09750806

Principal component regression

library(pls)
comps <- 1

Coefficients

\[ \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.0006077055 -0.0004782704 0.0003097628 0.0004990806
t(pcr_coef(overlap_x_xts, overlap_y_xts, comps))
          [,1]        [,2]        [,3]       [,4]
[1,] 0.2976711 -0.02689162 0.003195201 0.01504472
fit <- pcr(reformulate(termlabels = ".", response = tickers), 
           data = merge(scale_x_xts, overlap_y_xts), ncomp = comps)
coef(fit)[ , , 1]
        SP500    DTWEXAFEGS         DGS10  BAMLH0A0HYM2 
 0.0006077055 -0.0004782704  0.0003097628  0.0004990806

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)
    
    return(result)
    
}
pcr_rsq(scale_x_xts, overlap_y_xts, comps)
          BAICX
BAICX 0.6555066
pcr_rsq(overlap_x_xts, overlap_y_xts, comps)
         BAICX
BAICX 0.456853
R2(fit)$val[comps + 1]
[1] 0.6555066

Standard 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} \]

# unable to verify the result
pcr_se <- function(x, y, comps) {
    
    n_rows <- nrow(x)
    n_cols <- ncol(x)
    
    rsq <- pcr_rsq(x, y, comps)
    
    y <- sweep(y, 2, colMeans(y), "-")
    
    df_resid <- n_rows - n_cols - 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.803152e-05 2.206109e-05 1.428837e-05 2.302100e-05
pcr_se(overlap_x_xts, overlap_y_xts, comps)
[1] 0.0206518316 0.0018656873 0.0002216767 0.0010437727

Partial least squares

Risk decomposition

Standalone risk

\[ \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 <- model.matrix(~ 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))
    
    return(result)
    
}
lm_sar(overlap_x_xts, overlap_y_xts, weights, intercept) * sqrt(scale[["periods"]] * scale[["overlap"]])
[1] 0.056078921 0.000000000 0.024119624 0.007299644 0.029654129 0.012305724
[7] 0.027607781

Risk contribution

\[ \begin{aligned} \text{MCR}&=w^\mathrm{T}\frac{\partial\sigma_{y}}{\partial w}\\ &=w^\mathrm{T}\frac{\Sigma w}{\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 <- model.matrix(~ 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)
    
    return(result)
    
}
lm_mcr(overlap_x_xts, overlap_y_xts, weights, intercept) * sqrt(scale[["periods"]] * scale[["overlap"]])
[1] 0.056078921 0.000000000 0.015956899 0.003914750 0.017939151 0.004676747
[7] 0.013591374

Scenario analysis

Implied shocks

\[ \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.0103585 -0.005015706

Stress P&L

pnl_stress <- function(shocks, x, y, z, weights, intercept) {
    
    coef <- lm_coef(x, y, weights, intercept)
    
    if (intercept) x <- model.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)      xSP500 xDTWEXAFEGS      xDGS10 xBAMLH0A0HYM2
BAICX -0.0001049635 -0.01913491  -0.0145976 -0.03346309  -0.006121339