Merge pull request #192 from bashtage/estimate-ewma

ENH: Allow estimation of EWMA Variance parameters

Author: bashtage

arch/tests/univariate/test_volatility.py

@@ -1,10 +1,9 @@
-import warnings
-from unittest import TestCase
-
 import numpy as np
 import pytest
+import warnings
 from numpy.testing import assert_almost_equal, assert_equal, assert_allclose, \
     assert_array_equal
+from unittest import TestCase
 
 try:
     from arch.univariate import _recursions as rec
@@ -697,6 +696,55 @@ def test_ewma(self):
         txt = ewma.__repr__()
         assert str(hex(id(ewma))) in txt
 
+    def test_ewma_estimated(self):
+        ewma = EWMAVariance(lam=None)
+
+        sv = ewma.starting_values(self.resids)
+        assert sv == 0.94
+        assert sv.shape[0] == ewma.num_params
+
+        bounds = ewma.bounds(self.resids)
+        assert len(bounds) == 1
+        assert bounds == [(0, 1)]
+
+        var_bounds = ewma.variance_bounds(self.resids)
+
+        backcast = ewma.backcast(self.resids)
+        w = 0.94 ** np.arange(75)
+        assert_almost_equal(backcast,
+                            np.sum((self.resids[:75] ** 2) * (w / w.sum())))
+
+        parameters = np.array([0.9234])
+
+        var_bounds = ewma.variance_bounds(self.resids)
+        ewma.compute_variance(parameters, self.resids, self.sigma2, backcast, var_bounds)
+        cond_var_direct = np.zeros_like(self.sigma2)
+        cond_var_direct[0] = backcast
+        parameters = np.array([0, 1-parameters[0], parameters[0]])
+        rec.garch_recursion(parameters,
+                            self.resids ** 2.0,
+                            np.sign(self.resids),
+                            cond_var_direct,
+                            1, 0, 1, self.T, backcast, var_bounds)
+        assert_allclose(self.sigma2, cond_var_direct)
+        assert_allclose(self.sigma2 / cond_var_direct, np.ones_like(self.sigma2))
+
+        names = ewma.parameter_names()
+        names_target = ['lam']
+        assert_equal(names, names_target)
+
+        a, b = ewma.constraints()
+        a_target = np.ones((1, 1))
+        b_target = np.zeros((1,))
+        assert_array_equal(a, a_target)
+        assert_array_equal(b, b_target)
+
+        assert_equal(ewma.num_params, 1)
+        assert_equal(ewma.name, 'EWMA/RiskMetrics')
+        assert isinstance(ewma.__str__(), str)
+        txt = ewma.__repr__()
+        assert str(hex(id(ewma))) in txt
+
     def test_riskmetrics(self):
         rm06 = RiskMetrics2006()
 

arch/univariate/volatility.py

@@ -6,7 +6,6 @@
 from __future__ import absolute_import, division
 
 import itertools
-
 import numpy as np
 from numpy import (sqrt, ones, zeros, isscalar, sign, ones_like, arange, empty, abs, array, finfo,
                    float64, log, exp, floor)
@@ -37,6 +36,7 @@ class BootstrapRng(object):
     start : int
         Location of first forecast
     """
+
     def __init__(self, std_resid, start):
         if start <= 0 or start > std_resid.shape[0]:
             raise ValueError('start must be > 0 and <= len(std_resid).')
@@ -1294,8 +1294,9 @@ class EWMAVariance(VolatilityProcess):
 
     Parameters
     ----------
-    lam : float, optional
-        Smoothing parameter. Default is 0.94
+    lam : {float, None}, optional
+        Smoothing parameter. Default is 0.94. Set to None to estimate lam
+        jointly with other model parameters
 
     Attributes
     ----------
@@ -1317,38 +1318,52 @@ class EWMAVariance(VolatilityProcess):
 
         \sigma_t^{2}=\lambda\sigma_{t-1}^2 + (1-\lambda)\epsilon^2_{t-1}
 
-    This model has no parameters since the smoothing parameter is fixed.
+    When lam is provided, this model has no parameters since the smoothing
+    parameter is treated as fixed. Sel lam to `None` to jointly estimate this
+    parameter when fitting the model.
     """
 
     def __init__(self, lam=0.94):
         super(EWMAVariance, self).__init__()
         self.lam = lam
-        self.num_params = 0
-        if not 0.0 < lam < 1.0:
+        self._estimate_lam = lam is None
+        self.num_params = 1 if self._estimate_lam else 0
+        if lam is not None and not 0.0 < lam < 1.0:
             raise ValueError('lam must be strictly between 0 and 1')
         self.name = 'EWMA/RiskMetrics'
 
     def __str__(self):
-        descr = self.name + '(lam: ' + '{0:0.2f}'.format(self.lam) + ')'
+        if self._estimate_lam:
+            descr = self.name + '(lam: Estimated)'
+        else:
+            descr = self.name + '(lam: ' + '{0:0.2f}'.format(self.lam) + ')'
         return descr
 
     def starting_values(self, resids):
-        return np.empty((0,))
+        if self._estimate_lam:
+            return np.array([0.94])
+        return np.array([])
 
     def parameter_names(self):
+        if self._estimate_lam:
+            return ['lam']
         return []
 
-    def variance_bounds(self, resids, power=2.0):
-        return ones((resids.shape[0], 1)) * array([-1.0, finfo(float64).max])
-
     def bounds(self, resids):
+        if self._estimate_lam:
+            return [(0, 1)]
         return []
 
     def compute_variance(self, parameters, resids, sigma2, backcast,
                          var_bounds):
-        return ewma_recursion(self.lam, resids, sigma2, resids.shape[0], backcast)
+        lam = parameters[0] if self._estimate_lam else self.lam
+        return ewma_recursion(lam, resids, sigma2, resids.shape[0], backcast)
 
     def constraints(self):
+        if self._estimate_lam:
+            a = ones((1, 1))
+            b = zeros((1,))
+            return a, b
         return np.empty((0, 0)), np.empty((0,))
 
     def simulate(self, parameters, nobs, rng, burn=500, initial_value=None):
@@ -1362,7 +1377,11 @@ def simulate(self, parameters, nobs, rng, burn=500, initial_value=None):
 
         sigma2[0] = initial_value
         data[0] = sqrt(sigma2[0])
-        lam, one_m_lam = self.lam, 1.0 - self.lam
+        if self._estimate_lam:
+            lam = parameters[0]
+        else:
+            lam = self.lam
+        one_m_lam = 1.0 - lam
         for t in range(1, nobs + burn):
             sigma2[t] = lam * sigma2[t - 1] + one_m_lam * data[t - 1] ** 2.0
             data[t] = np.sqrt(sigma2[t]) * errors[t]
@@ -1393,7 +1412,10 @@ def _simulation_forecast(self, parameters, resids, backcast, var_bounds, start,
         paths.fill(np.nan)
         shocks = np.empty((t, simulations, horizon))
         shocks.fill(np.nan)
-        lam = self.lam
+        if self._estimate_lam:
+            lam = parameters[0]
+        else:
+            lam = self.lam
 
         for i in range(start, t):
             std_shocks = rng((simulations, horizon))

doc/source/changes/4.0.txt

@@ -8,3 +8,4 @@ Changes since 4.0
   a ``ZeroMean`` variance process.
 - Fixed a bug that prevented ``fix`` from being used with a new model (:issue:`156`)
 - Added ``first_obs`` and ``last_obs`` parameters to ``fix`` to mimic ``fit``
+- Added ability to jointly estimate smoothing parameter in EWMA variance when fitting the model

examples/univariate_volatility_modeling.ipynb

@@ -542,6 +542,7 @@
     "    * GJR-GARCH (`GARCH` using `o` argument) \n",
     "    * TARCH/ZARCH (`GARCH` using `power` argument set to `1`)\n",
     "    * Power GARCH and Asymmetric Power GARCH (`GARCH` using `power`)\n",
+    "    * Exponentially Weighted Moving Average Variance with estimated coefficient (`EWMAVariance`)\n",
     "    * Heterogeneous ARCH (`HARCH`)\n",
     "    * Parameterless Models\n",
     "        * Exponentially Weighted Moving Average Variance, known as RiskMetrics (`EWMAVariance`)\n",