Risk

analysis
finance
python
Author

Jason Foster

Published

August 7, 2026

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

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).

import math
import statsmodels.api as sm

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

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.87415931e-04,  2.08120413e-01, -1.51095218e-01,  2.60153176e+00,
        1.25777282e+00])
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.87415931e-04,  2.08120413e-01, -1.51095218e-01,  2.60153176e+00,
        1.25777282e+00])

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

def lm_rsq(x, y, weights, intercept):
          
  coef = 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.6752057381748476)
fit.rsquared
np.float64(0.6752057381748475)

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

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.02241748e-05, 2.30934491e-02, 4.33090175e-02, 3.34593361e-01,
       3.49378212e-01])
np.array(fit.bse)
array([6.02241748e-05, 2.30934491e-02, 4.33090175e-02, 3.34593361e-01,
       3.49378212e-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 combinations 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)) \]

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 pd.Series(result, index = x.columns)
lm_shap(overlap_x_df, overlap_y_df, weights, intercept)
SP500           0.306750
DTWEXAFEGS      0.108573
DGS10           0.122555
BAMLH0A0HYM2    0.137328
dtype: float64

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:

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)
  _, V = np.linalg.eigh(np.cov(x.T, ddof = 1))
  V = V[:, ::-1]
  
  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([ 0.00055299, -0.0004465 ,  0.00018054,  0.00048137])
pcr_coef(overlap_x_df, overlap_y_df, comps)
array([ 0.3049027 , -0.07029764,  0.00323122,  0.01239269])
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([ 0.00055299, -0.0004465 ,  0.00018054,  0.00048137])

R-squared

def pcr_rsq(x, y, comps):
  
  coef = 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.6161576642566288)
pcr_rsq(overlap_x_df, overlap_y_df, comps)
np.float64(0.5809505963812225)
fit = LinearRegression().fit(pca_x_df[:, :comps], overlap_y_df)
fit.score(pca_x_df[:, :comps], overlap_y_df)
0.6161576642566291

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

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:

def pcr_se(x, y, comps):

  n_rows = x.shape[0]

  rsq = pcr_rsq(x, y, comps)

  y = y - np.average(y, axis = 0)

  df_resid = n_rows - comps - 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([2.76042396e-05, 2.22887554e-05, 9.01213758e-06, 2.40291373e-05])
pcr_se(overlap_x_df, overlap_y_df, comps)
array([0.01637775, 0.00377602, 0.00017356, 0.00066567])
fit = sm.OLS(overlap_y_df.values, sm.add_constant(pca_x_df[:, :comps])).fit()

V = pca.components_.T[:, :comps]
np.sqrt(np.dot(V, np.dot(fit.cov_params()[1:, 1:], V.T)).diagonal())
array([2.76042396e-05, 2.22887554e-05, 9.01213758e-06, 2.40291373e-05])

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

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(([sigma[-1, -1]], sar, [sar_eps])))

  return pd.Series(result, index = ["total"] + list(x.columns) + ["eps"])
lm_sar(overlap_x_df, overlap_y_df, weights, intercept) * np.sqrt(scale["periods"] * scale["overlap"])
total           0.056788
const           0.000000
SP500           0.028092
DTWEXAFEGS      0.008350
DGS10           0.017092
BAMLH0A0HYM2    0.010314
eps             0.032364
dtype: float64

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

def lm_mcr(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)
  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(sigma[-1, -1])], mcr, [mcr_eps]))

  return pd.Series(result, index = ["total"] + list(x.columns) + ["eps"])
lm_mcr(overlap_x_df, overlap_y_df, weights, intercept) * np.sqrt(scale["periods"] * scale["overlap"])
total           0.056788
const           0.000000
SP500           0.021106
DTWEXAFEGS      0.004312
DGS10           0.007299
BAMLH0A0HYM2    0.005627
eps             0.018444
dtype: float64

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

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

  beta = np.linalg.lstsq(np.multiply(np.sqrt(weights), z), np.multiply(np.sqrt(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.00246393, -0.00390347])

Stress 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:

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.00018201, -0.02081204, -0.01510952, -0.00640998, -0.00490967])