Risk

analysis
finance
python
Author

Jason Foster

Published

December 13, 2025

factors_r = ["SP500", "DTWEXAFEGS"] # "SP500" does not contain dividends; note: "DTWEXM" discontinued as of Jan 2020
factors_d = ["DGS10", "BAMLH0A0HYM2"]
tickers = ["BAICX"] # fund inception date is "2011-11-28"
intercept = True

Regression analysis

import math
import statsmodels.api as sm

Ordinary least squares

Coefficients

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

def lm_coef(x, y, weights, intercept):
  
  if (intercept): x = sm.add_constant(x)
      
  result = np.dot(np.linalg.inv(np.dot(x.T, np.multiply(weights, x))),
                  np.dot(x.T, np.multiply(weights, y)))
  
  return np.ravel(result)
lm_coef(overlap_x_df, overlap_y_df, weights, intercept)
array([ 1.78898582e-04,  1.94426816e-01, -1.73838822e-01,  2.41633767e+00,
        5.40986160e-01])
if (intercept): overlap_x_df = sm.add_constant(overlap_x_df)
    
fit = sm.WLS(overlap_y_df, overlap_x_df, weights = weights).fit()

if (intercept): overlap_x_df = overlap_x_df.iloc[:, 1:]

np.array(fit.params)
array([ 1.78898582e-04,  1.94426816e-01, -1.73838822e-01,  2.41633767e+00,
        5.40986160e-01])

R-squared

\[ \begin{aligned} R^{2}=\frac{\hat{\beta}^\mathrm{T}(X^\mathrm{T}WX)\hat{\beta}}{y^\mathrm{T}Wy} \end{aligned} \]

def lm_rsq(x, y, weights, intercept):
          
  coef = np.matrix(lm_coef(x, y, weights, intercept))
  
  if (intercept):
    
    x = sm.add_constant(x)
    x = x - np.average(x, axis = 0, weights = weights.reshape(-1))
    y = y - np.average(y, axis = 0, weights = weights.reshape(-1))
      
  result = np.dot(coef, np.dot(np.dot(x.T, np.multiply(weights, x)), coef.T)) / \
    np.dot(y.T, np.multiply(weights, y))
  
  return np.float64(result.item())
lm_rsq(overlap_x_df, overlap_y_df, weights, intercept)
np.float64(0.5875188575304348)
fit.rsquared
np.float64(0.5875188575304344)

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

def lm_se(x, y, weights, intercept):
  
  n_rows = x.shape[0]
  n_cols = x.shape[1]
  
  rsq = lm_rsq(x, y, weights, intercept)
  
  if (intercept):
    
    x = sm.add_constant(x)
    y = y - np.average(y, axis = 0, weights = weights.reshape(-1))
    
    df_resid = n_rows - n_cols - 1 
    
  else:
    df_resid = n_rows - n_cols        
  
  var_y = np.dot(y.T, np.multiply(weights, y))
  var_resid = (1 - rsq) * var_y / df_resid
  
  result = np.sqrt(var_resid * np.linalg.inv(np.dot(x.T, np.multiply(weights, x))).diagonal())
  
  return np.ravel(result)
lm_se(overlap_x_df, overlap_y_df, weights, intercept)
array([6.57047023e-05, 2.12783402e-02, 3.94319790e-02, 3.51282714e-01,
       2.69524602e-01])
np.array(fit.bse)
array([6.57047023e-05, 2.12783402e-02, 3.94319790e-02, 3.51282714e-01,
       2.69524602e-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)) \]

def lm_shap(x, y, weights, intercept):

  n_rows = x.shape[0]
  n_cols = x.shape[1]
  n_combn = 2 ** n_cols
  n_vec = np.zeros(n_combn)
  ix_mat = np.zeros((n_cols, n_combn))
  rsq = np.zeros(n_combn)
  result = np.zeros(n_cols)
  
  # number of binary combinations
  for k in range(n_combn):
    
    n = 0
    n_size = k
    
    # find the binary combination
    for j in range(n_cols):
      
      if n_size % 2 == 0:
        
        n += 1
        
        ix_mat[j, k] = j + 1
          
      n_size //= 2
    
    n_vec[k] = n
    
    if n > 0:
      
      ix_subset = np.where(ix_mat[:, k] != 0)[0]
      x_subset = x.iloc[:, ix_subset]
      
      rsq[k] = lm_rsq(x_subset, y, weights, intercept)

  # calculate the exact Shapley value for r-squared
  for j in range(n_cols):
    
    ix_pos = np.where(ix_mat[j, :] != 0)[0]
    ix_neg = np.where(ix_mat[j, :] == 0)[0]
    ix_n = n_vec[ix_neg]
    rsq_diff = rsq[ix_pos] - rsq[ix_neg]

    for k in range(int(n_combn / 2)):
      
      s = int(ix_n[k])
      weight = math.factorial(s) * math.factorial(n_cols - s - 1) \
        / math.factorial(n_cols)
      result[j] += weight * rsq_diff[k]

  return result
lm_shap(overlap_x_df, overlap_y_df, weights, intercept)
array([0.28785857, 0.06534543, 0.12088291, 0.11343195])

Principal component regression

from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression
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} \]

def pcr_coef(x, y, comps):
  
  x = x - np.average(x, axis = 0)
  L, V = np.linalg.eig(np.cov(x.T, ddof = 1))
  idx = L.argsort()[::-1]
  V = V[:, idx]
  
  W = np.dot(x, V)
  gamma = np.dot(np.dot(np.linalg.inv(np.dot(W.T, W)), W.T), y)
  
  result = np.dot(V[:, :comps], gamma[:comps])
  
  return np.ravel(result)
scale_x_df = (overlap_x_df - np.average(overlap_x_df, axis = 0)) \
  / np.std(overlap_x_df, axis = 0, ddof = 1)
pcr_coef(scale_x_df, overlap_y_df, comps)
array([ 5.54019179e-04, -1.97715761e-04,  8.15907881e-05,  5.40199199e-04])
pcr_coef(overlap_x_df, overlap_y_df, comps)
array([ 0.24520833, -0.01359597,  0.00134507,  0.01331329])
pca = PCA(n_components = len(factors))
pca_x_df = pca.fit_transform(scale_x_df)

fit = LinearRegression(fit_intercept = False).fit(pca_x_df, overlap_y_df)

gamma = fit.coef_
np.dot(pca.components_.T[:, :comps], gamma.T[:comps]).ravel()
array([ 5.54019179e-04, -1.97715761e-04,  8.15907881e-05,  5.40199199e-04])

R-squared

def pcr_rsq(x, y, comps):
  
  coef = np.matrix(pcr_coef(x, y, comps))
  
  x = x - np.average(x, axis = 0)
  y = y - np.average(y, axis = 0)
  
  result = np.dot(np.dot(coef, np.dot(x.T, x)), coef.T) / np.dot(y.T, y)
  
  return np.float64(result.item())
pcr_rsq(scale_x_df, overlap_y_df, comps)
np.float64(0.46508317129844273)
pcr_rsq(overlap_x_df, overlap_y_df, comps)
np.float64(0.44814853224120876)

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
def pcr_se(x, y, comps):
  
  n_rows = x.shape[0]
  n_cols = x.shape[1]
  
  rsq = pcr_rsq(x, y, comps)
  
  y = y - np.average(y, axis = 0)
  
  df_resid = n_rows - n_cols - 1
  
  var_y = np.dot(y.T, y)   
  var_resid = (1 - rsq) * var_y / df_resid
  
  # uses statsmodels for illustrative purposes
  pca = sm.multivariate.PCA(x, standardize = False, demean = True)
  L = pca.eigenvals[:comps]
  V = pca.eigenvecs.iloc[:, :comps]
  
  result = np.sqrt(var_resid * np.dot(V, np.dot(np.diag(1 / L), V.T)).diagonal())
  
  return np.ravel(result)
pcr_se(scale_x_df, overlap_y_df, comps)
array([3.78054264e-05, 1.34918229e-05, 5.56763132e-06, 3.68623720e-05])
pcr_se(overlap_x_df, overlap_y_df, comps)
array([1.73135794e-02, 9.59979334e-04, 9.49721846e-05, 9.40019726e-04])

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

def cov_wt(x, weights, center):
  
  sum_w = sum(weights)
  sumsq_w = sum(np.power(weights, 2))
  
  if (center):
  
    x = x - np.average(x, axis = 0, weights = weights.reshape(-1))
  
  result = np.dot(x.T, np.multiply(weights, x)) / (sum_w - sumsq_w / sum_w)
  
  return result
def lm_sar(x, y, weights, intercept):
    
  coef = lm_coef(x, y, weights, intercept)
  rsq = lm_rsq(x, y, weights, intercept)
  
  if (intercept): x = sm.add_constant(x)
  
  # sigma = np.cov(np.concatenate((x, y), axis = 1).T,
  #                aweights = weights.reshape(-1))
  sigma = cov_wt(np.concatenate((x, y), axis = 1), weights, intercept)
  sar = np.multiply(np.power(coef, 2).T, sigma[:-1, :-1].diagonal())
  sar_eps = (1 - rsq) * sigma[-1, -1]
  
  result = np.sqrt(np.concatenate((np.matrix(sigma[-1, -1]),
                                   np.matrix(sar),
                                   np.matrix(sar_eps)), axis = 1))
  
  return np.ravel(result)
lm_sar(overlap_x_df, overlap_y_df, weights, intercept) * np.sqrt(scale["periods"] * scale["overlap"])
array([0.05487256, 0.        , 0.02895832, 0.01029764, 0.01654892,
       0.00638248, 0.03524174])

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

def lm_mcr(x, y, weights, intercept):
    
  coef = np.matrix(lm_coef(x, y, weights, intercept)).T
  rsq = lm_rsq(x, y, weights, intercept)
      
  if (intercept): x = sm.add_constant(x)
  
#     sigma = np.cov(np.concatenate((x, y), axis = 1).T,
#                    aweights = weights.reshape(-1))
  sigma = cov_wt(np.concatenate((x, y), axis = 1), weights, intercept)
  mcr = np.multiply(coef, np.dot(sigma[:-1, :-1], coef)) / np.sqrt(sigma[-1, -1])
  mcr_eps = np.sqrt(sigma[-1, -1]) - sum(mcr)
  
  result = np.concatenate((np.sqrt(np.matrix(sigma[-1, -1])),
                           np.matrix(mcr).T,
                           np.matrix(mcr_eps)), axis = 1)
  
  return np.ravel(result)
lm_mcr(overlap_x_df, overlap_y_df, weights, intercept) * np.sqrt(scale["periods"] * scale["overlap"])
array([0.05487256, 0.        , 0.01922233, 0.0034666 , 0.00656938,
       0.00298036, 0.0226339 ])

Scenario analysis

Implied shocks

\[ \begin{aligned} \hat{\beta}&=(Z^\mathrm{T}WZ)^{-1}Z^\mathrm{T}WX \end{aligned} \]

def implied_shocks(shocks, x, z, weights):

  beta = np.linalg.lstsq(np.multiply(weights, z), np.multiply(weights, x), rcond = None)[0]
                   
  result = np.dot(shocks, beta)
  
  return result
shocks = np.array([-0.1, 0.1])
overlap_z_df = overlap_x_df.iloc[:, [0, 1]]
implied_shocks(shocks, overlap_x_df, overlap_z_df, weights)
array([-0.1       ,  0.1       , -0.00329699, -0.00594851])

Stress P&L

def pnl_stress(shocks, x, y, z, weights, intercept):
  
  coef = lm_coef(x, y, weights, intercept)
  
  if (intercept): x = sm.add_constant(x)
  
  result = np.multiply(coef.T, implied_shocks(shocks, x, z, weights))
  
  return result
pnl_stress(shocks, overlap_x_df, overlap_y_df, overlap_z_df, weights, intercept)
array([-0.00180301, -0.01944268, -0.01738388, -0.00796663, -0.00321806])