Update GARCH-CCC starting value

docs/source/algorithms/garch-ccc.rst

@@ -186,10 +186,14 @@ Estimation techniques
 My implementation of :class:`frds.algorithms.GARCHModel_CCC` fits the GARCH-DCC model 
 by simultaneously estimating all parameters via maxmimizing the log-likelihood :math:numref:`log_likelihood`.
 
+General steps are:
+
 1. Use :class:`frds.algorithms.GARCHModel` to estimate the :doc:`/algorithms/garch` model for each of the returns.
 
-2. Use as starting vaues the estimates and a correlation found to maximize loglikelihood.
-`
+2. Use the standardized residuals from the estimated GARCH models to compute correlation coefficient.
+
+3. Use as starting vaues the estimated parameters from above in optimizing the loglikelihood function.
+
 
 References
 ==========
@@ -228,16 +232,16 @@ Use :class:`frds.algorithms.GARCHModel_CCC` to estimate a GARCH(1,1)-CCC.
 >>> res = model_ccc.fit()
 >>> from pprint import pprint
 >>> pprint(res)
-Parameters(mu1=0.02746938560178172,
-           omega1=0.034013928248590015,
-           alpha1=0.06593333053443852,
-           beta1=0.9219586517939616,
-           mu2=0.009412154182986385,
-           omega2=0.05869540397594033,
-           alpha2=0.08305499592375533,
-           beta2=0.9040973431326755,
-           rho=0.6506782161411894,
-           loglikelihood=-7281.321453325305)
+Parameters(mu1=0.02745814255283541,
+           omega1=0.03401400758840226,
+           alpha1=0.06593379740524756,
+           beta1=0.9219575443861723,
+           mu2=0.009390068254041505,
+           omega2=0.058694325049554734,
+           alpha2=0.0830561828957614,
+           beta2=0.9040961791372522,
+           rho=0.6506770477876749,
+           loglikelihood=-7281.321453218112)
 
 These results are comparable to the ones obtained in Stata, and even marginally 
 better based on log-likelihood. In Stata, we can estimate the same model as below:

docs/source/index.rst

@@ -9,11 +9,17 @@ of ready-to-use methods for computing a wide array of measures in the literature
 It is developed by Dr. `Mingze Gao <http://mingze-gao.com>`_  from the University of Sydney, as 
 a personal project during his postdoctoral research fellowship.
 
-|GitHub license|
+|GitHub license| |PyPI Downloads| |Tests|
 
 .. |GitHub license| image:: https://img.shields.io/github/license/mgao6767/frds?color=blue
    :target: https://github.com/mgao6767/frds/blob/master/LICENSE
 
+.. |PyPI Downloads| image:: https://img.shields.io/pypi/dm/frds?label=PyPI%20downloads 
+   :target: https://pypi.org/project/frds/
+
+.. |Tests| image:: https://github.com/mgao6767/frds/actions/workflows/test.yml/badge.svg
+   :target: https://github.com/mgao6767/frds/actions/workflows/test.yml
+
 .. |frds| replace:: :code:`frds`
 
 .. important:: 
@@ -36,16 +42,15 @@ installed via ``pip``.
 
 The structure of |frds| is simple:
 
-* :mod:`frds.measures` provides a collection of measures.
 * :mod:`frds.algorithms` provides a collection of algorithms.
+* :mod:`frds.measures` provides a collection of measures.
 * :mod:`frds.datasets` provides example datasets.
 
 --------
 Examples
 --------
 
-|frds| aims to provide out-of-the-box.
-
+Some simple examples.
 
 Measure
 -------
@@ -74,8 +79,8 @@ liabilities.
 ...         [ 0.283,  0.053, -0.085, 0.508, -0.370,  1.000],
 ...     ]
 ... )
->>> dip = DistressInsurancePremium()
->>> dip.estimate(default_probabilities, correlations)       
+>>> dip = DistressInsurancePremium(default_probabilities, correlations)
+>>> dip.estimate()
 0.2865733550799999
 
 Algorithm
@@ -107,16 +112,16 @@ The results are as good as those otained in Stata, if not better.
 >>> model_ccc = GARCHModel_CCC(toyota, nissan)
 >>> res = model_ccc.fit()
 >>> pprint(res)
-Parameters(mu1=0.02746938560178172,
-           omega1=0.034013928248590015,
-           alpha1=0.06593333053443852,
-           beta1=0.9219586517939616,
-           mu2=0.009412154182986385,
-           omega2=0.05869540397594033,
-           alpha2=0.08305499592375533,
-           beta2=0.9040973431326755,
-           rho=0.6506782161411894,
-           loglikelihood=-7281.321453325305)
+Parameters(mu1=0.02745814255283541,
+           omega1=0.03401400758840226,
+           alpha1=0.06593379740524756,
+           beta1=0.9219575443861723,
+           mu2=0.009390068254041505,
+           omega2=0.058694325049554734,
+           alpha2=0.0830561828957614,
+           beta2=0.9040961791372522,
+           rho=0.6506770477876749,
+           loglikelihood=-7281.321453218112)
 
 ---------
 Read more

src/frds/algorithms/_mgarch.py

@@ -47,15 +47,20 @@ def __init__(self, returns1: np.ndarray, returns2: np.ndarray) -> None:
         self.parameters = type(self).Parameters()
 
     def fit(self) -> Parameters:
-        """Estimates the Multivariate GARCH(1,1) parameters via MLE
+        """Estimates the Multivariate GARCH(1,1)-CCC parameters via MLE
 
         Returns:
             params: :class:`frds.algorithms.GARCHModel_CCC.Parameters`
         """
         m1, m2 = self.model1, self.model2
         m1.fit()
         m2.fit()
-        starting_vals = self.starting_values(m1.resids, m2.resids)
+        z1 = m1.resids / np.sqrt(m1.sigma2)
+        z2 = m2.resids / np.sqrt(m2.sigma2)
+        rho = np.corrcoef(z1, z2)[1, 0]
+        m1_params = list(asdict(m1.parameters).values())[:-1]
+        m2_params = list(asdict(m2.parameters).values())[:-1]
+        starting_vals = [*m1_params, *m2_params, rho]
 
         # Step 4. Set bounds for parameters
         bounds = [
@@ -174,41 +179,6 @@ def loglikelihood(
         log_likelihood = np.sum(log_likelihood_terms)
         return log_likelihood
 
-    def starting_values(self, resids1: np.ndarray, resids2: 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
-        combination of alpha and beta that gives the highest loglikelihood.
-
-        Args:
-            resids1 (np.ndarray): Array of residuals for the first return series.
-            resids2 (np.ndarray): Array of residuals for the second return series.
-
-        Returns:
-            List[float]: [mu1, omega1, alpha1, beta1, mu2, omega2, alpha2, beta2, rho]
-        """
-        m1, m2 = self.model1, self.model2
-        # Constant correlation
-        rho_grid = np.linspace(-0.9, 0.9, 10)
-
-        initial_params = []
-        max_likelihood = -np.inf
-        for rho in rho_grid:
-            # last one is loglikelihood
-            m1_params = list(asdict(m1.parameters).values())[:-1]
-            m2_params = list(asdict(m2.parameters).values())[:-1]
-            params = [*m1_params, *m2_params, rho]
-            ll = -self.loglikelihood_model(
-                params,
-                m1.backcast_value,
-                m2.backcast_value,
-                m1.var_bounds,
-                m2.var_bounds,
-            )
-            if ll > max_likelihood:
-                initial_params = params
-
-        return initial_params
-
 
 if __name__ == "__main__":
     import pandas as pd
@@ -223,6 +193,6 @@ def starting_values(self, resids1: np.ndarray, resids2: np.ndarray) -> List[floa
     nissan = df["nissan"].to_numpy() * 100
     honda = df["honda"].to_numpy() * 100
 
-    model = GARCHModel_CCC(toyota, honda)
+    model = GARCHModel_CCC(toyota, nissan)
     res = model.fit()
     pprint(res)

tests/algorithms/test_garch_ccc.py

@@ -1,43 +1,28 @@
 import pytest
-import numpy as np
-
 from frds.algorithms import GARCHModel_CCC
+from frds.datasets import StockReturns
 
 
 def test_garch_ccc():
-    rng = np.random.default_rng(42)
-
-    # Initialize parameters
-    T = 500  # Number of observations
-    omega1, alpha1, beta1 = 0.1, 0.2, 0.7  # Parameters for first series
-    omega2, alpha2, beta2 = 0.1, 0.3, 0.6  # Parameters for second series
-    rho = 0.5  # Constant correlation
-
-    # Initialize variables
-    y1 = np.zeros(T)
-    y2 = np.zeros(T)
-    h1 = np.zeros(T)
-    h2 = np.zeros(T)
-    z1 = rng.normal(size=T)
-    z2 = rng.normal(size=T)
-
-    # Simulate conditional variances and returns
-    for t in range(1, T):
-        h1[t] = omega1 + alpha1 * y1[t - 1] ** 2 + beta1 * h1[t - 1]
-        h2[t] = omega2 + alpha2 * y2[t - 1] ** 2 + beta2 * h2[t - 1]
+    returns = StockReturns.stocks_us
 
-        # Apply constant correlation
-        e1 = np.sqrt(h1[t]) * z1[t]
-        e2 = np.sqrt(h2[t]) * (rho * z1[t] + np.sqrt(1 - rho**2) * z2[t])
+    sp500 = returns["^GSPC"].to_numpy() * 100
+    googl = returns["GOOGL"].to_numpy() * 100
 
-        y1[t] = e1
-        y2[t] = e2
+    model_ccc = GARCHModel_CCC(sp500, googl)
+    res = model_ccc.fit()
 
-    model = GARCHModel_CCC(y1, y2)
-    res = model.fit()
-    print(res)
+    tol = 0.01
+    assert res.mu1 == pytest.approx(0.0699378, rel=tol)
+    assert res.omega1 == pytest.approx(0.0585878, rel=tol)
+    assert res.alpha1 == pytest.approx(0.1477404, rel=tol)
+    assert res.beta1 == pytest.approx(0.7866691, rel=tol)
+    assert res.mu2 == pytest.approx(0.0940275, rel=tol)
+    assert res.omega2 == pytest.approx(0.4842512, rel=tol)
+    assert res.alpha2 == pytest.approx(0.12166, rel=tol)
+    assert res.beta2 == pytest.approx(0.7113389, rel=tol)
+    assert res.rho == pytest.approx(0.6646705, rel=tol)
 
 
 if __name__ == "__main__":
-    # pytest.main([__file__])
-    test_garch_ccc()
+    pytest.main([__file__])