Risk

analysis
finance
python
Author

Jason Foster

Published

September 1, 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([ 2.08748558e-04,  2.23680226e-01, -1.24298709e-01,  3.28341494e+00,
        1.24750998e+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([ 2.08748558e-04,  2.23680226e-01, -1.24298709e-01,  3.28341494e+00,
        1.24750998e+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.7123821960605949)
fit.rsquared
np.float64(0.7123821960605949)

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([5.71925778e-05, 2.15649088e-02, 4.23118152e-02, 3.44422436e-01,
       3.35937116e-01])
np.array(fit.bse)
array([5.71925778e-05, 2.15649088e-02, 4.23118152e-02, 3.44422436e-01,
       3.35937116e-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)) \]

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.328151
DTWEXAFEGS      0.093712
DGS10           0.151418
BAMLH0A0HYM2    0.139101
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.00056532, -0.00045043,  0.00016783,  0.00050129])
pcr_coef(overlap_x_df, overlap_y_df, comps)
array([ 0.31474314, -0.06562135,  0.00280488,  0.01264038])
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.00056532, -0.00045043,  0.00016783,  0.00050129])

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.6190134302078931)
pcr_rsq(overlap_x_df, overlap_y_df, comps)
np.float64(0.5958854799924306)
fit = LinearRegression().fit(pca_x_df[:, :comps], overlap_y_df)
fit.score(pca_x_df[:, :comps], overlap_y_df)
0.6190134302078932

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.80498615e-05, 2.23494149e-05, 8.32747970e-06, 2.48727797e-05])
pcr_se(overlap_x_df, overlap_y_df, comps)
array([0.01639295, 0.0034178 , 0.00014609, 0.00065836])
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.80498615e-05, 2.23494149e-05, 8.32747970e-06, 2.48727797e-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.057470
const           0.000000
SP500           0.030238
DTWEXAFEGS      0.006592
DGS10           0.019856
BAMLH0A0HYM2    0.010114
eps             0.030821
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.057470
const           0.000000
SP500           0.023083
DTWEXAFEGS      0.003223
DGS10           0.009120
BAMLH0A0HYM2    0.005515
eps             0.016529
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.00224874, -0.00415565])

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.00041162, -0.02236802, -0.01242987, -0.00738353, -0.00518422])