Merge branch 'master' into master

arch/bootstrap/base.py

@@ -461,6 +461,7 @@ def conf_int(self, func, reps=1000, method='basic', size=0.95, tail='two',
                                   'too small to use BCa'
                         raise RuntimeError(message.format(jk_var=denom))
                     a = numer / denom
+                    a = np.atleast_1d(a)
                     a = a[:, None]
                 else:
                     a = 0.0

arch/tests/bootstrap/test_bootstrap.py

@@ -457,6 +457,19 @@ def func(y):
         ci = ci.T
         assert_allclose(ci_db, ci)
 
+    def test_conf_int_bca_scaler(self):
+        num_bootstrap = 100
+        bs = IIDBootstrap(self.y)
+        bs.seed(23456)
+
+        try:
+            ci = bs.conf_int(np.mean, reps=num_bootstrap, method='bca')
+            assert(ci.shape == (2, 1))
+        except IndexError:
+            pytest.fail('conf_int(method=\'bca\') scaler input regression. '
+                        'Ensure output is at least 1D with '
+                        'numpy.atleast_1d().')
+
     def test_conf_int_parametric(self):
         def param_func(x, params=None, state=None):
             if state is not None:

arch/tests/univariate/test_mean.py

@@ -395,7 +395,12 @@ def test_arch_arx(self):
                       'omega', 'alpha[1]', 'alpha[2]'])
 
         am = ARX(y=y, lags=2, x=x)
-        am.fit(disp=DISPLAY).summary()
+        res = am.fit(disp=DISPLAY)
+        summ = res.summary().as_text()
+        repr = res.__repr__()
+        assert str(hex(id(res))) in repr
+        assert summ[:10] == repr[:10]
+
         am.volatility = ARCH(p=2)
         results = am.fit(update_freq=0, disp='off')
         assert isinstance(results.pvalues, pd.Series)

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
@@ -776,6 +775,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/base.py

@@ -747,7 +747,23 @@ def forecast(self, params, horizon=1, start=None, align='origin', method='analyt
         raise NotImplementedError('Subclasses must implement')
 
 
-class ARCHModelFixedResult(object):
+class _SummaryRepr(object):
+    """Base class for returning summary as repr and str"""
+
+    def summary(self):
+        return Summary()
+
+    def __repr__(self):
+        out = self.__str__() + '\n'
+        out += self.__class__.__name__
+        out += ', id: {0}'.format(hex(id(self)))
+        return out
+
+    def __str__(self):
+        return self.summary().as_text()
+
+
+class ARCHModelFixedResult(_SummaryRepr):
     """
     Results for fixed parameters for an ARCHModel model
 
@@ -1653,12 +1669,12 @@ class ARCHModelForecast(object):
     simulations : ARCHModelForecastSimulation
         Object containing detailed simulation results if using a simulation-based method
     """
+
     def __init__(self, index, mean, variance, residual_variance,
                  lterm_residual_variance=None,
                  simulated_paths=None, simulated_variances=None,
                  simulated_residual_variances=None, simulated_residuals=None,
                  align='origin'):
-
         mean = _format_forecasts(mean, index)
         variance = _format_forecasts(variance, index)
         residual_variance = _format_forecasts(residual_variance, index)

arch/univariate/recursions.pyx

@@ -14,6 +14,8 @@ cdef extern from 'math.h':
 cdef extern from 'float.h':
     double DBL_MAX
 
+cdef double LNSIGMA_MAX = log(DBL_MAX)
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 @cython.cdivision(True)
@@ -265,14 +267,13 @@ def egarch_recursion(double[:] parameters,
             else:
                 lnsigma2[t] += parameters[loc] * lnsigma2[t - 1 - j]
             loc += 1
+        if lnsigma2[t] > LNSIGMA_MAX:
+            lnsigma2[t] = LNSIGMA_MAX
         sigma2[t] = exp(lnsigma2[t])
         if sigma2[t] < var_bounds[t, 0]:
             sigma2[t] = var_bounds[t, 0]
         elif sigma2[t] > var_bounds[t, 1]:
-            if sigma2[t] > DBL_MAX:
-                sigma2[t] = var_bounds[t, 1] + 1000
-            else:
-                sigma2[t] = var_bounds[t, 1] + log(sigma2[t] / var_bounds[t, 1])
+            sigma2[t] = var_bounds[t, 1] + log(sigma2[t]) - log(var_bounds[t, 1])
         std_resids[t] = resids[t] / sqrt(sigma2[t])
         abs_std_resids[t] = fabs(std_resids[t])
 

arch/univariate/recursions_python.py

@@ -13,6 +13,8 @@
 __all__ = ['harch_recursion', 'arch_recursion', 'garch_recursion',
            'egarch_recursion', 'cgarch_recursion']
 
+LNSIGMA_MAX = np.log(np.finfo(np.double).max) - .1
+
 
 def harch_recursion_python(parameters, resids, sigma2, lags, nobs, backcast,
                            var_bounds):
@@ -225,14 +227,13 @@ def egarch_recursion_python(parameters, resids, sigma2, p, o, q, nobs,
             else:
                 lnsigma2[t] += parameters[loc] * lnsigma2[t - 1 - j]
             loc += 1
+        if lnsigma2[t] > LNSIGMA_MAX:
+            lnsigma2[t] = LNSIGMA_MAX
         sigma2[t] = np.exp(lnsigma2[t])
         if sigma2[t] < var_bounds[t, 0]:
             sigma2[t] = var_bounds[t, 0]
         elif sigma2[t] > var_bounds[t, 1]:
-            if not np.isinf(sigma2[t]):
-                sigma2[t] = var_bounds[t, 1] + log(sigma2[t] / var_bounds[t, 1])
-            else:
-                sigma2[t] = var_bounds[t, 1] + 1000
+            sigma2[t] = var_bounds[t, 1] + log(sigma2[t]) - log(var_bounds[t, 1])
         std_resids[t] = resids[t] / np.sqrt(sigma2[t])
         abs_std_resids[t] = np.abs(std_resids[t])
 

arch/univariate/volatility.py

@@ -5,7 +5,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).')
@@ -1318,8 +1318,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
     ----------
@@ -1341,38 +1342,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):
@@ -1386,7 +1401,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]
@@ -1417,7 +1436,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,5 +8,6 @@ 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
 - Added Component GARCH process which decomposes volatility into long run and short run components.
   The model comes with an extra ``longterm_var`` attribute when fit.

examples/univariate_volatility_modeling.ipynb

@@ -543,6 +543,7 @@
     "    * CGARCH (`CGARCH`)\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",