Long-Run Marginal Expected Shortfall (LRMES)¶
Introduction¶
LRMES is used in estimating SRISK by Brownlees and Engle (2017). It measures the expected firm return conditional on a systemic event.
A systemic event is a market decline below a threshold C over a time horizon h.
If the multiperiod arithmetic market return between t+1 to t+h is R_{mt+1:t+h}, then the systemic event is \{R_{mt+1:t+h}<C\}.
LRMES_{it} for firm i at time t is then defined as
Brownlees and Engle (2017) use a GARCH-DCC model to construct the LRMES predictions.
GARCH-DCC¶
Let firm and market log returns be r_{it}=\log(1+R_{it}) and r_{mt}=\log(1+R_{mt}). Conditional on the information set \mathcal{F}_{t−1} available at time t−1, the return pair has an (unspecified) distribution \mathcal{D} with zero mean and time–varying covariance,
To specify the evolution of the time varying volatilities and correlation, GJR-GARCH volatility model and the standard DCC correlation model are used.
GARCH for volatility¶
Specifically, the GJR-GARCH model equations for the volatility dynamics are:
where I^-_{it} = 1 if r_{it}<0 and I^-_{mt} = 1 if r_{mt}<0.
DCC for correlation¶
The DCC models correlation through the volatility-adjusted returns \epsilon=r/\sigma.
where \mathbf Q_{it} is the so-called pseudo correlation matrix.
The DCC then models the dynamics of \mathbf Q_{it} as
where,
- \bar{\mathbf{Q}}_i is the unconditional correlation matrix of the firm and market adjusted returns, and
- \mathbf e_{t-1}=\begin{bmatrix} \epsilon_{it-1} \\ \epsilon_{mt-1} \end{bmatrix}
Estimating GARCH-DCC¶
The above model is typically estimated by a two–step QML estimation procedure. More extensive details on this modeling approach and estimation are provided in Engle (2009).
Note by Mingze
Equation 4.33 in Engle (2009) shows that the log likelihood can be additively divided into two parts, one concerns the variance and the other concerns correlation. Therefore, we can solve for the variance and correlation parameters in two separate steps, hence "two-step" QML.
Specifically, we first estimate a GJR-GARCH(1,1) for each firm (and market)'s log return series to obtain the conditional volatilities \sigma and hence \epsilon=r/\sigma. In the second step, we use the estimated coefficients to estimate the DCC model for \epsilon for each pair of firm returns and market returns.
Note
Refer to algorithms/dcc_garch for the estimation process.
Predicting LRMES¶
Appendix A. Simulation Algorithm for LRMES in Engle (2009) describes the exact steps to construct LRMES forecasts.
Summary by Mingze
The general idea is to simulate market returns and use the estimated GARCH-DCC model to derive the corresponding firm returns. We then use the distribution of returns to estimate LRMES.
- Step 1. Construct GARCH-DCC standardized innovations for the training sample t=1,...,T, where \xi_{it} is standardized, linearly orthogonal shocks of the firm to the market on day t,
-
Step 2. Sample with replacement S\times h pairs of [\xi_{it}, \epsilon_{mt}]', which are used as the simulated innovations from time T+1 to T+h. Notice that in the algorithm, the innovations are simulated by resampling the standardized residuals of the GARCH-DCC rather than relying on parametric assumptions.
-
Step 3. Use the pseudo sample of innovations as inputs of the DCC and GARCH filters, respectively. Initial conditions are the last values of the conditional correlation \rho_{iT} and variances \sigma^2_{iT} and \sigma^2_{mT}. This step delivers S pseudo samples of GARCH-DCC logarithmic returns from period T+1 to period T+h, conditional on the realized process up to time T, that is
Note by Mingze
Suppose we have a simulated volatility-adjusted market return \epsilon^s_{mT+h} at time T+h, then the corresponding volatility-adjusted firm return at T+h is computed as
Therefore, we need to predict \rho_{iT+h}, i.e., the off-diagnal element of \mathbf R_{T+h} = \begin{bmatrix} 1 & \rho_{iT+h} \\ \rho_{iT+h} & 1 \end{bmatrix}.
Note that \mathbf R_{T+h} = \text{diag}(\mathbf Q_{iT+h})^{-1/2} \mathbf Q_{iT+h} \text{diag}(\mathbf Q_{iT+h})^{-1/2}, and
We therefore need to make assumptions about E_T[\mathbf{e}_{T+h-1} \mathbf{e}'_{T+h-1}] because these are the volatility-adjusted returns in the future, but we don't have future returns.
Since, E_T[\mathbf e_{T+h-1} \mathbf e'_{T+h-1}]=E_T[R_{T+h-1}], the assumption we make here is that
- \bar{\mathbf R} \approx \bar{\mathbf Q}, and
- E_T[\mathbf R_{T+h}] \approx E_T[\mathbf Q_{T+h}]
According to Engle and Sheppard (2001), this assumption seems to provide better bias properties.
So,
where,
- \bar{\mathbf R}_i = \text{diag}(\bar{\mathbf{Q}}_{i})^{-1/2} \bar{\mathbf{Q}}_{i} \text{diag}(\bar{\mathbf{Q}}_{i})^{-1/2}
- E_T[\mathbf R_{T+1}]= \text{diag}(\hat{\mathbf{Q}}_{iT+1})^{-1/2} \hat{\mathbf{Q}}_{iT+1} \text{diag}(\hat{\mathbf{Q}}_{iT+1})^{-1/2},
Further,
and that \mathbf{e}_T = [\epsilon_{iT}, \epsilon_{mT}]' and \mathbf Q_{iT} are known.
- Step 4. Construct the multiperiod arithmetic firm (market) return of each pseudo sample,
- Step 5. Compute LRMES as the Monte Carlo average of the simulated multiperiod arithmetic returns conditional on the systemic event,
API¶
estimate(firm_returns, market_returns, h=22, S=10000, C=-0.1, random_seed=42)
¶
h-step-ahead LRMES forecasts conditional on a systemic event of market decline C
Parameters:
Name | Type | Description | Default |
---|---|---|---|
firm_returns | np.ndarray | (n_days,n_firms) array of firm log returns. | required |
market_returns | np.ndarray | (n_days,) array of market log returns. | required |
h | int | h-period-ahead prediction horizon. Defaults to 22. | 22 |
S | int | sample size used in simulation to generate LRMES forecasts. Defaults to 10000. | 10000 |
C | float | market decline used to define systemic event. Defaults to -0.1, i.e. -10%. | -0.1 |
random_seed | int | random seed. Defaults to 42. | 42 |
Returns:
Type | Description |
---|---|
np.ndarray | np.ndarray: (n_firms,) array of LRMES forecasts |
Examples:
>>> import frds.measures.long_run_mes as lrmes
>>> import yfinance as yf
>>> import numpy as np
>>> df = yf.download(tickers="SPY JPM GS", start="2017-01-01", end="2022-12-31", progress=False)
>>> df = df[['Adj Close']]
>>> df.columns = df.columns.droplevel(0)
>>> df
GS JPM SPY
Date
2017-01-03 214.042953 72.665413 202.085266
2017-01-04 215.425156 72.799454 203.287506
2017-01-05 213.821426 72.129333 203.126007
2017-01-06 216.993469 72.137695 203.852783
2017-01-09 215.212524 72.187935 203.179840
... ... ... ...
2022-12-23 343.053650 129.302628 381.454193
2022-12-27 339.538818 129.755707 379.949921
2022-12-28 338.446625 130.464844 375.227936
2022-12-29 340.988434 131.213409 381.982178
2022-12-30 340.938812 132.080154 380.975983
[1510 rows x 3 columns]
>>> mkt_returns = np.log(df.SPY.pct_change()+1).dropna().to_numpy()
>>> firm_returns = np.array([np.log(df.JPM.pct_change()+1).dropna().to_numpy(),
... np.log(df.GS.pct_change()+1).dropna().to_numpy()]).T
>>> lrmes.estimate(firm_returns, mkt_returns)
array([ 0.12958814, -0.09460028])
Note
yfinance.download
may lead to different data at each run, so the results can vary. However, given a fixed input dataset, the function will yield the same estimate conditional on the same random seed.
Hence, the 22-day LRMES for JPM is 12.96% and -9.5% for GS. The latter would be ignored in computing SRISK, though.
Source code in src/frds/measures/long_run_mes.py
def estimate(
firm_returns: np.ndarray,
market_returns: np.ndarray,
h=22,
S=10_000,
C=-0.1,
random_seed=42,
) -> np.ndarray:
"""h-step-ahead LRMES forecasts conditional on a systemic event of market decline C
Args:
firm_returns (np.ndarray): (n_days,n_firms) array of firm log returns.
market_returns (np.ndarray): (n_days,) array of market log returns.
h (int, optional): h-period-ahead prediction horizon. Defaults to 22.
S (int, optional): sample size used in simulation to generate LRMES forecasts. Defaults to 10000.
C (float, optional): market decline used to define systemic event. Defaults to -0.1, i.e. -10%.
random_seed (int, optional): random seed. Defaults to 42.
Returns:
np.ndarray: (n_firms,) array of LRMES forecasts
Examples:
>>> import frds.measures.long_run_mes as lrmes
>>> import yfinance as yf
>>> import numpy as np
>>> df = yf.download(tickers="SPY JPM GS", start="2017-01-01", end="2022-12-31", progress=False)
>>> df = df[['Adj Close']]
>>> df.columns = df.columns.droplevel(0)
>>> df
GS JPM SPY
Date
2017-01-03 214.042953 72.665413 202.085266
2017-01-04 215.425156 72.799454 203.287506
2017-01-05 213.821426 72.129333 203.126007
2017-01-06 216.993469 72.137695 203.852783
2017-01-09 215.212524 72.187935 203.179840
... ... ... ...
2022-12-23 343.053650 129.302628 381.454193
2022-12-27 339.538818 129.755707 379.949921
2022-12-28 338.446625 130.464844 375.227936
2022-12-29 340.988434 131.213409 381.982178
2022-12-30 340.938812 132.080154 380.975983
[1510 rows x 3 columns]
>>> mkt_returns = np.log(df.SPY.pct_change()+1).dropna().to_numpy()
>>> firm_returns = np.array([np.log(df.JPM.pct_change()+1).dropna().to_numpy(),
... np.log(df.GS.pct_change()+1).dropna().to_numpy()]).T
>>> lrmes.estimate(firm_returns, mkt_returns)
array([ 0.12958814, -0.09460028])
!!! note
`yfinance.download` may lead to different data at each run, so the results can vary.
However, given a fixed input dataset, the function will yield the same estimate conditional on the same random seed.
Hence, the 22-day LRMES for JPM is 12.96% and -9.5% for GS.
The latter would be ignored in computing SRISK, though.
"""
if len(firm_returns.shape) == 1:
firm_returns = firm_returns.reshape((firm_returns.shape[0], 1))
n_days, n_firms = firm_returns.shape
(_n_days,) = market_returns.shape
assert n_days == _n_days
# Fit GJR-GARCH for the market returns
mkt_am = arch_model(market_returns, p=1, o=1, q=1, rescale=True)
mkt_res = mkt_am.fit(update_freq=0, disp=False)
epsilon_mkt = market_returns / mkt_res.conditional_volatility
# Forecasted volatility
mkt_vol_hat = np.sqrt(
np.squeeze(
mkt_res.forecast(
mkt_res.params, horizon=h, reindex=False
).variance.T.to_numpy()
)
) # (h,1) array of volatility forecasts
lrmes = np.zeros((n_firms,))
for i in range(n_firms):
# Fit GJR-GARCH for each firm's returns
firm_am = arch_model(firm_returns[:, i], p=1, o=1, q=1, rescale=True)
firm_res = firm_am.fit(update_freq=0, disp=False)
epsilon_firm = firm_returns[:, i] / firm_res.conditional_volatility
# Forecasted volatility
firm_vol_hat = np.sqrt(
np.squeeze(
firm_res.forecast(
firm_res.params, horizon=h, reindex=False
).variance.T.to_numpy()
)
) # (h,1) array of volatility forecasts
# Estimate DCC for each firm-market pair
epsilon = np.array([epsilon_firm, epsilon_mkt])
a, b = dcc(epsilon) # params in DCC model
Q_avg = calc_Q_avg(epsilon)
Q = calc_Q(epsilon, a, b) # Qt for training data
R = calc_R(epsilon, a, b) # Rt for training data
# DCC predictions for correlations
et = epsilon[:, -1]
Q_Tplus1 = (1 - a - b) * Q_avg + a * np.outer(et, et) + b * Q[-1]
diag_q = 1.0 / np.sqrt(np.abs(Q_Tplus1))
diag_q = diag_q * np.eye(2)
R_Tplus1 = np.dot(np.dot(diag_q, Q_Tplus1), diag_q)
diag_q = 1.0 / np.sqrt(np.abs(Q_avg))
diag_q = diag_q * np.eye(2)
R_avg = np.dot(np.dot(diag_q, Q_avg), diag_q)
Rhat = []
for _h in range(1, h + 1):
_ab = np.power(a + b, _h - 1)
# R_Tplus_h is correlation matrix for T+h, symmetric 2*2 matrix
R_Tplus_h = (1 - _ab) * R_avg + _ab * R_Tplus1
# In case predicted correlation is larger than 1
if abs(R_Tplus_h[0, 1]) >= 1:
R_Tplus_h[0, 1] = 0.9999 * (1 if R_Tplus_h[0, 1] >= 0 else -1)
R_Tplus_h[1, 0] = R_Tplus_h[0, 1]
Rhat.append(R_Tplus_h)
# Sample innovations
innov = np.zeros((n_days,))
for t in range(n_days):
rho = R[t][0, 1]
innov[t] = (epsilon_firm[t] - epsilon_mkt[t] * rho) / np.sqrt(1 - rho**2)
sample = np.array([epsilon_mkt, innov]) # shape=(2,n_days)
# Sample S*h pairs of standardized innovations
rng = np.random.RandomState(random_seed)
# sample.shape=(S,h,2)
sample = sample.T[rng.choice(sample.shape[1], (S, h), replace=True)]
# list of simulated firm total returns when there're systemic events
firm_total_return = []
for s in range(S):
# Each simulation
inv = sample[s, :, :] # (h,2)
# mkt log return = mkt innovation * predicted mkt volatility
mktrets = np.multiply(inv[:, 0], mkt_vol_hat)
# systemic event if over the prediction horizon,
# the market falls by more than C
systemic_event = np.exp(np.sum(mktrets)) - 1 < C
# no need to simulate firm returns if there is no systemic event
if not systemic_event:
continue
# when there is a systemic event
firmrets = np.zeros((h,))
for _h in range(h):
mktinv, firminv = inv[_h][0], inv[_h][1]
# Simulated firm return at T+h
rho = Rhat[_h][0, 1]
firmret = firminv * np.sqrt(1 - rho**2) + rho * mktinv
firmret = firm_vol_hat[_h] * firmret
firmrets[_h] = firmret
# firm return over the horizon
firmret = np.exp(np.sum(firmrets)) - 1
firm_total_return.append(firmret)
# Store result
if len(firm_total_return):
lrmes[i] = np.mean(firm_total_return)
else:
lrmes[i] = np.nan
return lrmes
TODO¶
- Use
multiprocessing
for execution speed. - Use C++ to rewrite for execution speed.
References¶
- Brownlees and Engle (2017), SRISK: A Conditional Capital Shortfall Measure of Systemic Risk, Review of Financial Studies, 30 (1), 48–79.
- Duan and Zhang (2015), Non-Gaussian Bridge Sampling with an Application, SSRN.
- Orskaug (2009), Multivariate dcc-garch model with various error distributions.
See Also¶
- SRISK
- Absorption Ratio
- Contingent Claim Analysis
- Distress Insurance Premium
- Marginal Expected Shortfall (MES)
- Systemic Expected Shortfall (SES)