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

$LRMES_{it}=-E_i[R_{it+1:t+h} | R_{mt+1:t+h} < C]$

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,

$\begin{bmatrix}r_{it} \\ r_{mt}\end{bmatrix} | \mathcal F_{t-1} \sim \mathcal D\left(\mathbf 0, \begin{bmatrix}\sigma_{it}^2 & \rho_{it}\sigma_{it}\sigma_{mt} \\ \rho_{it}\sigma_{it}\sigma_{mt} & \sigma_{mt}^2 \end{bmatrix}\right)$

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:

$\begin{align} \sigma_{it}^2 &= \omega_{i} + \alpha_{i} r^2_{it-1} + \gamma_{i} r^2_{it-1} I^-_{it-1} + \beta_{i} \sigma^2_{it-1} \\ \sigma_{mt}^2 &= \omega_{m} + \alpha_{m} r^2_{mt-1} + \gamma_{m} r^2_{mt-1} I^-_{mt-1} + \beta_{m} \sigma^2_{mt-1} \end{align}$

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

$\mathbf R_t = \text{Corr} \begin{pmatrix} \epsilon_{it} \\ \epsilon_{mt} \end{pmatrix} = \begin{bmatrix} 1 & \rho_{it} \\ \rho_{it} & 1 \end{bmatrix} = \text{diag}(\mathbf Q_{it})^{-1/2} \mathbf Q_{it} \text{diag}(\mathbf Q_{it})^{-1/2}$

where $\mathbf Q_{it}$ is the so-called pseudo correlation matrix.

The DCC then models the dynamics of $\mathbf Q_{it}$ as

$\mathbf Q_{it} = (1-a-b) \bar{\mathbf Q}_i + a \mathbf{e}_{t-1} \mathbf{e}'_{t-1} + b \mathbf Q_{it-1}$

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$ ,

$\epsilon_{mt} = \frac{r_{mt}}{\sigma_{mt}} \text{ and } \xi_{it} = \left(\frac{r_{it}}{\sigma_{it}} - \rho_{it} \epsilon_{mt}\right) / \sqrt{1-\rho^2_{it}}$

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

$\begin{bmatrix} r^s_{iT+t} \\ r^s_{mT+t} \end{bmatrix}_{t=1,...,h} | \mathcal{F}_{T}$

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

$\left(\sqrt{1-\rho^2_{iT+h}} \times \xi^s_{iT+h} + \rho_{iT+h} \epsilon^s_{mt}\right)$

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

$E_T[\mathbf Q_{iT+h}] = (1-a-b)\bar{\mathbf Q}_i + a E_T[\mathbf{e}_{T+h-1} \mathbf{e}'_{T+h-1}] + b E_T[\mathbf Q_{iT+h-1}]$

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,

$\begin{align} E_T[\mathbf{R}_{T+h}] &\approx E_T[\mathbf Q_{T+h}] \\ &= (1-a-b)\bar{\mathbf Q}_i + a E_T[\mathbf R_{T+h-1}] + b E_T[\mathbf Q_{iT+h-1}] \\ &\approx (1-a-b)\bar{\mathbf R}_i + (a + b) \mathbf R_{T+h-1} \\ &= \dots \\ &= (1-(a+b)^{h-1}) \bar{\mathbf R}_i + (a+b)^{h-1} E_T[\mathbf R_{T+1}] \end{align}$

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,

$\hat{\mathbf{Q}}_{T+1} = (1-a-b)\bar{\mathbf Q}_i + a \mathbf{e}_T \mathbf{e}'_T+b \mathbf Q_{iT}$

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,

$R^s_{iT+1:T+h} = \exp \left(\sum_{t=1}^{h} r^s_{iT+t} \right) -1$

Step 5. Compute LRMES as the Monte Carlo average of the simulated multiperiod arithmetic returns conditional on the systemic event,

$LRMES_{iT} = - \frac{\sum_{s=1}^S R^s_{iT+1:T+h} I(R^s_{mT+1:T+h}<C)}{\sum_{s=1}^S I(R^s_{mT+1:T+h}<C)}$

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.