LRMES updated. C++ backend for (GJR)GARCH DCC

docs/source/measures/long_run_mes.rst

@@ -92,9 +92,10 @@ where,
    volatility at time :math:`t`.
 
 -  :math:`z_{it}, z_{mt}` follows an unknown bivariate distribution with
-   zero mean and some covariance structure.
+   zero mean and maybe some covariance structure.
 
-However, generally we can assume :math:`z_{it}, z_{mt}\sim N(0,1)`.
+However, generally we can assume :math:`z_{it}` and :math:`z_{mt}` are i.i.d. 
+standard normal
 
 .. important::
 
@@ -144,6 +145,14 @@ Specifically, the GJR-GARCH models the conditional variances as:
 where :math:`I^-_{it} = 1` if :math:`r_{it} < 0` and :math:`I^-_{mt} =
 1` if :math:`r_{mt} < 0`.
 
+.. note::
+
+   `Brownlees and Engle (2017) <https://doi.org/10.1093/rfs/hhw060>`_ use a 
+   zero-mean assumption, so that :math:`r` is used in the above equation 
+   instead of the residual :math:`\epsilon`. They use :math:`\epsilon` to denote
+   the standardized residual, which in this note is :math:`z` to be consistent
+   with other notes of ``frds``.
+
 Dynamic correlation
 ===================
 
@@ -175,8 +184,8 @@ the time-varying :math:`\mathbf{R}_t` via a proxy process
    :label: eq:dcc_corr
 
    \mathbf{R}_t = \text{Corr}\begin{pmatrix}
-   \epsilon_{it} \\\\
-   \epsilon_{mt}
+   z_{it} \\\\
+   z_{mt}
    \end{pmatrix}
    = \begin{bmatrix}
    1 & \rho_{it} \\\\
@@ -217,7 +226,7 @@ The Q process is updated as follows:
 .. math::
    :label: q11
 
-   q_{it} = (1 - a - b) \overline{q}_{i} + a z^2_{1,t-1} + b q_{i,t-1}
+   q_{it} = (1 - a - b) \overline{q}_{i} + a z^2_{i,t-1} + b q_{i,t-1}
 
 .. math::
    :label: q22
@@ -234,7 +243,7 @@ The dynamic conditional correlation :math:`\rho_t` is given by:
 .. math::
    :label: rho_t
 
-   \rho_{it} = \frac{q_{imt}}{\sqrt{q_{i} q_{mt}}}
+   \rho_{it} = \frac{q_{imt}}{\sqrt{q_{it} q_{mt}}}
 
 Estimating GJR-GARCH-DCC
 ========================
@@ -319,11 +328,13 @@ assumptions.
 
 .. note::
 
-   Basically, we need innovations :math:`z_{it}` and :math:`z_{mt}` in order to forecast.
-   To do so, after estimating the GJR-GARCH-DCC model, we use the standardized residuals
-   from the market time series and the estimated conditinal correlations to compute the "innovation" 
-   to the firm series *such that* its standarized residuals conforms the conditional correlations.
+   Basically, we need innovations :math:`z_{it}` and :math:`z_{mt}` for :math:`t>T` in order to forecast.
 
+   For the market returns' :math:`z_{mt}` in the prediction horizon, we sample the standardized residuals 
+   from the past with replacement. For the firm return innovations :math:`z_{it}` in the prediction horizon, 
+   we use the computed :math:`\xi_{it}`, which is *independent* from :math:`z_{mt}`. This guarantees that
+   the market and firm return innovations are i.i.d. standard normal, while the return covariances conform 
+   to the DCC model.
 
 **Step 3**. Use the pseudo sample of innovations as inputs of the DCC
 and (GJR)GARCH filters, respectively. Initial conditions are the last values
@@ -341,6 +352,44 @@ up to time :math:`T`, that is
       r^s_{mT+t}
    \end{bmatrix}_{t=1,...,h} | \mathcal{F}_{T}
 
+Specifically, in each simulation :math:`s`, for the :math:`h`-step-ahead prediction, 
+compute :math:`\hat{\sigma}^2_{iT+h}` and :math:`\hat{\sigma}^2_{mT+h}`,
+
+.. math:: 
+
+   \begin{align}
+   \hat{\sigma}^2_{iT+h} &= \omega_i + \left[\alpha_i+\gamma_i I(\xi_{ih}<0)\right] \xi_{ih}^2 + \beta_i \hat{\sigma}^2_{iT+h-1} \\\\
+   \hat{\sigma}^2_{mT+h} &= \omega_m + \left[\alpha_m+\gamma_m I(z_{mh}<0)\right] z_{mh}^2 + \beta_m \hat{\sigma}^2_{mT+h-1} 
+   \end{align}
+
+Then, update DCC coefficients,
+
+.. math::
+
+   \begin{align}
+   \hat{q}_{iT+h} &= (1 - a - b) \overline{q}_{i} + a \xi^2_{i,h-1} + b \hat{q}_{i,T+h-1} \\\\
+   \hat{q}_{mT+h} &= (1 - a - b) \overline{q}_{m} + a z^2_{m,h-1} + b \hat{q}_{m,T+h-1} \\\\
+   \hat{q}_{imT+h} &= (1 - a - b) \overline{q}_{im} + a \xi_{i,h-1} z_{m,h-1} + b \hat{q}_{im,T+h-1}
+   \end{align}
+
+The dynamic conditional correlation :math:`\hat{\rho}_{iT+1}` is given by:
+
+.. math::
+
+   \hat{\rho}_{iT+h} = \frac{\hat{q}_{imT+h}}{\sqrt{\hat{q}_{iT+h} \hat{q}_{mT+h}}}
+
+This conditional correlation :math:`\hat{\rho}_{iT+1}` is then used to compute the 1-step-ahead returns 
+given the conditional variances and innovations,
+
+.. math::
+
+   \begin{align}
+   \hat{r}_{iT+h} &= \mu_i + \hat{\rho}_{iT+h} \xi_{ih} \\\\
+   \hat{r}_{mT+h} &= \mu_m + \hat{\rho}_{iT+h} z_{mh}
+   \end{align}
+
+This process is performed for :math:`T+2,\dots,T+h` forecasts, so tha we have in this simulation
+:math:`s` a set of market and firm (log) returns, :math:`r^s_{iT+t}` and :math:`r^s_{mT+t}`, :math:`t=1,\dots,h`.
 
 **Step 4**. Construct the multiperiod arithmetic firm (market) return of
 each pseudo sample,
@@ -385,6 +434,30 @@ multiperiod arithmetic returns conditional on the systemic event,
 
 .. autoclass:: frds.measures.LRMES
 
+.. autoclass:: frds.measures.LongRunMarginalExpectedShortfall
+
 **********
  Examples
-**********
+**********
+
+Import built-in dataset of stock returns.
+
+>>> from frds.datasets import StockReturns
+>>> returns = StockReturns.stocks_us
+>>> returns.head()
+               GOOGL        GS       JPM     ^GSPC
+Date                                              
+2010-01-05 -0.004404  0.017680  0.019370  0.003116
+2010-01-06 -0.025209 -0.010673  0.005494  0.000546
+2010-01-07 -0.023280  0.019569  0.019809  0.004001
+2010-01-08  0.013331 -0.018912 -0.002456  0.002882
+2010-01-11 -0.001512 -0.015776 -0.003357  0.001747
+
+Estimate LRMES using the last 600 days' observtions.
+
+>>> from frds.measures import LRMES
+>>> gs = returns["GS"].to_numpy()[-600:]
+>>> sp500 = returns["^GSPC"].to_numpy()[-600:]
+>>> lrmes = LRMES(gs, sp500)
+>>> lrmes.estimate(S=10000, h=22, C=-0.1)
+-0.02442517334543748

setup.py

@@ -21,17 +21,17 @@
     language="c++",
 )
 
-# mod_dcc = Extension(
-#     "frds.algorithms.dcc.dcc_ext",
-#     include_dirs=[get_path("platinclude"), numpy.get_include()],
-#     sources=["src/frds/algorithms/dcc/dcc_module.cpp"],
-#     extra_compile_args=extra_compile_args,
-#     language="c++",
-# )
+mod_algo_utils = Extension(
+    "frds.algorithms.utils.utils_ext",
+    include_dirs=[get_path("platinclude"), numpy.get_include()],
+    sources=["src/frds/algorithms/utils/utils_module.cpp"],
+    extra_compile_args=extra_compile_args,
+    language="c++",
+)
 
 ext_modules = [
     mod_isolation_forest,
-    # mod_dcc,
+    mod_algo_utils,
 ]
 
 setup(

src/frds/algorithms/_garch.py

@@ -5,6 +5,12 @@
 import numpy as np
 from scipy.optimize import minimize, OptimizeResult
 
+USE_CPP_EXTENSION = True
+try:
+    import frds.algorithms.utils.utils_ext as ext
+except ImportError:
+    USE_CPP_EXTENSION = False
+
 
 class GARCHModel:
     """:doc:`/algorithms/garch` model with constant mean and Normal noise
@@ -162,6 +168,8 @@ def compute_variance(
         Returns:
             np.ndarray: conditional variance
         """
+        if USE_CPP_EXTENSION:
+            return ext.compute_garch_variance(params, resids, backcast, var_bounds)
         omega, alpha, beta = params
         sigma2 = np.zeros_like(resids)
         sigma2[0] = omega + (alpha + beta) * backcast
@@ -252,6 +260,8 @@ def bounds_check(sigma2: float, var_bounds: np.ndarray) -> float:
         Returns:
             float: adjusted conditional variance
         """
+        if USE_CPP_EXTENSION:
+            return ext.bounds_check(sigma2, var_bounds)
         lower, upper = var_bounds[0], var_bounds[1]
         sigma2 = max(lower, sigma2)
         if sigma2 > upper:
@@ -292,6 +302,8 @@ def ewma(resids: np.ndarray, initial_value: float, lam=0.94) -> np.ndarray:
         Returns:
             np.ndarray: Array containing the conditional variance estimates.
         """
+        if USE_CPP_EXTENSION:
+            return ext.ewma(resids, initial_value, lam)
         T = len(resids)
         variance = np.empty(T)
         variance[0] = initial_value  # Set the initial value
@@ -384,6 +396,8 @@ def compute_variance(
         Returns:
             np.ndarray: conditional variance
         """
+        if USE_CPP_EXTENSION:
+            return ext.compute_gjrgarch_variance(params, resids, backcast, var_bounds)
         # fmt: off
         omega, alpha, gamma, beta = params
         sigma2 = np.zeros_like(resids)
@@ -394,6 +408,31 @@ def compute_variance(
             sigma2[t] = GJRGARCHModel.bounds_check(sigma2[t], var_bounds[t])
         return sigma2
 
+    @staticmethod
+    def forecast_variance(
+        params: Parameters,
+        resids: np.ndarray,
+        initial_variance: float,
+    ) -> np.ndarray:
+        """Forecast the variances conditional on given parameters and residuals.
+
+        Args:
+            params (Parameters): :class:`frds.algorithms.GJRGARCHModel.Parameters`
+            resids (np.ndarray): residuals to use
+            initial_variance (float): starting value of variance forecasts
+
+        Returns:
+            np.ndarray: conditional variance
+        """
+        # fmt: off
+        omega, alpha, gamma, beta = params.omega, params.alpha, params.gamma, params.beta
+        sigma2 = np.zeros_like(resids)
+        sigma2[0] = initial_variance
+        for t in range(1, len(resids)):
+            sigma2[t] = omega + alpha * (resids[t - 1] ** 2) + beta * sigma2[t - 1]
+            sigma2[t] += gamma * resids[t - 1] ** 2 if resids[t - 1] < 0 else 0
+        return sigma2[1:]
+
     def starting_values(self, resids: np.ndarray) -> List[float]:
         """Finds the optimal initial values for the volatility model via a grid
         search. For varying target persistence and alpha values, return the

src/frds/algorithms/_mgarch.py

@@ -320,11 +320,44 @@ def loglikelihood_model(self, params: np.ndarray) -> float:
         resids2 = self.model2.resids
         sigma2_1 = self.model1.sigma2
         sigma2_2 = self.model2.sigma2
+        # z1 and z2 are standardized residuals
+        z1 = resids1 / np.sqrt(sigma2_1)
+        z2 = resids2 / np.sqrt(sigma2_2)
 
         # The loglikelihood of the variance component (Step 1)
         l1 = self.model1.parameters.loglikelihood + self.model2.parameters.loglikelihood
 
         # The loglikelihood of the correlation component (Step 2)
+        rho = self.conditional_correlations(a, b)
+
+        log_likelihood_terms = -0.5 * (
+            - (z1**2 + z2**2)
+            + np.log(1 - rho ** 2)
+            + (z1 ** 2  + z2 ** 2  - 2 * rho * z1 * z2) / (1 - rho ** 2)
+        )
+        l2 = np.sum(log_likelihood_terms)
+
+        negative_loglikelihood = - (l1 + l2)
+        return negative_loglikelihood
+
+    def conditional_correlations(self, a: float, b: float) -> np.ndarray:
+        """Computes the conditional correlations based on given a and b.
+        Other parameters are
+
+        Args:
+            a (float): DCC parameter
+            b (float): DCC parameter
+
+        Returns:
+            np.ndarray: array of conditional correlations
+        """
+        self.model1.fit()  # in case it was not estimated, no performance loss
+        self.model2.fit()  # in case it was not estimated
+        resids1 = self.model1.resids
+        resids2 = self.model2.resids
+        sigma2_1 = self.model1.sigma2
+        sigma2_2 = self.model2.sigma2
+
         # z1 and z2 are standardized residuals
         z1 = resids1 / np.sqrt(sigma2_1)
         z2 = resids2 / np.sqrt(sigma2_2)
@@ -335,28 +368,18 @@ def loglikelihood_model(self, params: np.ndarray) -> float:
         q12 = np.empty_like(z1)
         q22 = np.empty_like(z1)
         rho = np.zeros_like(z1)
-        # TODO: initial values for Q
         q11[0] = q_11_bar
         q22[0] = q_22_bar
         q12[0] = q_12_bar
         rho[0] = q12[0] / np.sqrt(q11[0] * q22[0])
 
-        # TODO: fix RuntimeWarning about invalid values
         for t in range(1, T):
             q11[t] = (1 - a - b) * q_11_bar + a * z1[t - 1] ** 2 + b * q11[t - 1]
             q22[t] = (1 - a - b) * q_22_bar + a * z2[t - 1] ** 2 + b * q22[t - 1]
             q12[t] = (1 - a - b) * q_12_bar + a * z1[t - 1] * z2[t - 1] + b * q12[t - 1]
             rho[t] = q12[t] / np.sqrt(q11[t] * q22[t])
 
-        log_likelihood_terms = -0.5 * (
-            - (z1**2 + z2**2)
-            + np.log(1 - rho ** 2)
-            + (z1 ** 2  + z2 ** 2  - 2 * rho * z1 * z2) / (1 - rho ** 2)
-        )
-        l2 = np.sum(log_likelihood_terms)
-
-        negative_loglikelihood = - (l1 + l2)
-        return negative_loglikelihood
+        return rho
 
 
 class GJRGARCHModel_DCC(GARCHModel_DCC):