Merge 353f89f into 475ef9a

arch/tests/univariate/test_variance_forecasting.py

@@ -2,20 +2,23 @@
 
 import numpy as np
 import pytest
+
 try:
     from arch.univariate.recursions import garch_recursion
 except ImportError:
     from arch.univariate.recursions_python import garch_recursion
 from numpy.random import RandomState
 from numpy.testing import assert_allclose
+import pandas as pd
 
+from arch.univariate import arch_model
 from arch.univariate.distribution import Normal, StudentsT
+from arch.univariate.mean import ConstantMean
 from arch.univariate.volatility import GARCH, ConstantVariance, HARCH, EWMAVariance, \
     RiskMetrics2006, BootstrapRng, EGARCH, FixedVariance
 
 
 def _compare_truncated_forecasts(full, trunc, start):
-
     assert np.all(np.isnan(trunc.forecasts[:start]))
     assert_allclose(trunc.forecasts[start:], full.forecasts[start:])
 
@@ -30,6 +33,7 @@ class PreservedState(object):
     Context manager that will save NumPy's random generator's state when entering and restore
     the original state when exiting.
     """
+
     def __init__(self, random_state):
         self._random_state = random_state
         self._state = None
@@ -307,7 +311,7 @@ def test_arch_1_forecast_bootstrap(self):
         shocks.fill(np.nan)
         for i in range(333, 1000):
             locs = self.rng.random_sample((2000, 5))
-            int_locs = np.floor((i+1) * locs).astype(np.int64)
+            int_locs = np.floor((i + 1) * locs).astype(np.int64)
             std_shocks = std_resids[int_locs]
             paths[i, :, 0] = params[0] + params[1] * self.resid[i] ** 2.0
             shocks[i, :, 0] = np.sqrt(paths[i, :, 0]) * std_shocks[:, 0]
@@ -599,7 +603,7 @@ def test_harch_forecast_bootstrap(self):
         arch[:22] += 0.2 / 22
         for i in range(100, t):
             locs = self.rng.random_sample((1000, 10))
-            locs *= (i+1)
+            locs *= (i + 1)
             int_locs = np.floor(locs).astype(np.int64)
             std_shocks = std_resid[int_locs]
             temp = np.empty((1000, 32))
@@ -835,7 +839,7 @@ def test_egarch_211_forecast(self):
         sqrt2pi = np.sqrt(2 / np.pi)
         for i in range(731, t):
             locs = self.rng.random_sample((1111, 3))
-            int_locs = np.floor((i+1) * locs).astype(np.int64)
+            int_locs = np.floor((i + 1) * locs).astype(np.int64)
             std_shocks = std_resids[int_locs]
             paths[i, :, 0] = np.log(one_step[i])
 
@@ -934,7 +938,7 @@ def test_egarch_212_forecast_smoke(self):
         sqrt2pi = np.sqrt(2 / np.pi)
         for i in range(731, t):
             locs = self.rng.random_sample((1111, 3))
-            int_locs = np.floor((i+1) * locs).astype(np.int64)
+            int_locs = np.floor((i + 1) * locs).astype(np.int64)
             std_shocks = std_resids[int_locs]
             paths[i, :, 0] = np.log(one_step[i])
 
@@ -1252,11 +1256,11 @@ def test_ewma_forecast(self):
 
         forecast = vol.forecast(params, resids, backcast, var_bounds, horizon=10, start=0)
 
-        expected = np.zeros((t+1))
+        expected = np.zeros((t + 1))
         expected[0] = backcast
         lam = 0.94
-        for i in range(1, t+1):
-            expected[i] = lam * expected[i-1] + (1-lam) * resids[i-1] ** 2
+        for i in range(1, t + 1):
+            expected[i] = lam * expected[i - 1] + (1 - lam) * resids[i - 1] ** 2
         for i in range(10):
             assert_allclose(forecast.forecasts[:, i], expected[1:])
 
@@ -1289,7 +1293,7 @@ def test_ewma_simulation(self):
             shocks[i, :, 0] = np.sqrt(paths[i, :, 0]) * std_shocks[:, 0]
             for j in range(1, 10):
                 paths[i, :, j] = vol.lam * paths[i, :, j - 1] + \
-                    (1 - vol.lam) * shocks[i, :, j - 1] ** 2
+                                 (1 - vol.lam) * shocks[i, :, j - 1] ** 2
                 shocks[i, :, j] = np.sqrt(paths[i, :, j]) * std_shocks[:, j]
 
         assert_allclose(forecasts.forecasts, paths.mean(1))
@@ -1469,3 +1473,25 @@ def test_bs_rng_errors(self):
         y = self.rng.rand(1000)
         with pytest.raises(ValueError):
             BootstrapRng(y, 0)
+
+
+def test_external_rng():
+
+    arch_mod = arch_model(None, mean='Constant', vol='GARCH', p=1, q=1)
+    data = arch_mod.simulate(np.array([0.1, 0.1, 0.1, 0.88]), 1000)
+    data.index = pd.date_range('2000-01-01', periods=data.index.shape[0])
+
+    rand_state = np.random.RandomState(1)
+    state = rand_state.get_state()
+    volatility = GARCH(1, 0, 1)
+    distribution = Normal(rand_state)
+    mod = ConstantMean(data.data, volatility=volatility, distribution=distribution)
+    res = mod.fit(disp='off')
+
+    rand_state.set_state(state)
+    fcast_1 = res.forecast(res.params, horizon=5, method='simulation', start=900, simulations=250)
+    rand_state.set_state(state)
+    rng = rand_state.standard_normal
+    fcast_2 = res.forecast(res.params, horizon=5, method='simulation', start=900,
+                           simulations=250, rng=rng)
+    assert_allclose(fcast_1.residual_variance, fcast_2.residual_variance)

arch/univariate/base.py

@@ -652,7 +652,7 @@ def compute_param_cov(self, params, backcast=None, robust=True):
             return inv_hess / nobs
 
     def forecast(self, params, horizon=1, start=None, align='origin', method='analytic',
-                 simulations=1000):
+                 simulations=1000, rng=None):
         """
         Construct forecasts from estimated model
 
@@ -685,12 +685,16 @@ def forecast(self, params, horizon=1, start=None, align='origin', method='analyt
         simulations : int
             Number of simulations to run when computing the forecast using
             either simulation or bootstrap.
+        rng : callable, optional
+            Custom random number generator to use in simulation-based forecasts.
+            Must produce random samples using the syntax `rng(size)` where size
+            is a tuple containing the dimension of the random values.
 
         Returns
         -------
         forecasts : DataFrame
-            t by h data frame containing the forecasts.  The alignment of the forecasts is
-            controlled by `align`.
+            t by h data frame containing the forecasts.  The alignment of the
+            forecasts is controlled by `align`.
 
         Examples
         --------
@@ -1048,7 +1052,7 @@ def plot(self, annualize=None, scale=None):
         return fig
 
     def forecast(self, params=None, horizon=1, start=None, align='origin', method='analytic',
-                 simulations=1000):
+                 simulations=1000, rng=None):
         """
         Construct forecasts from estimated model
 
@@ -1072,21 +1076,26 @@ def forecast(self, params=None, horizon=1, start=None, align='origin', method='a
             from time t-1, the 2 step from time t-2, ..., and the h-step from
             time t-h.  'target' simplified computing forecast errors since the
             realization and h-step forecast are aligned.
-        method : {'analytic', 'simulation', 'bootstrap'}
+        method : {'analytic', 'simulation', 'bootstrap'}, optional
             Method to use when producing the forecast. The default is analytic.
             The method only affects the variance forecast generation.  Not all
             volatility models support all methods. In particular, volatility
             models that do not evolve in squares such as EGARCH or TARCH do not
             support the 'analytic' method for horizons > 1.
-        simulations : int
+        simulations : int, optional
             Number of simulations to run when computing the forecast using
             either simulation or bootstrap.
+        rng : callable, optional
+            Custom random number generator to use in simulation-based forecasts.
+            Must produce random samples using the syntax `rng(size)` where size
+            is either a scalar or a tuple containing the dimension of the
+            random values.
 
         Returns
         -------
         forecasts : DataFrame
-            t by h data frame containing the forecasts.  The alignment of the forecasts is
-            controlled by `align`.
+            t by h data frame containing the forecasts.  The alignment of the
+            forecasts is controlled by `align`.
 
         Notes
         -----
@@ -1115,7 +1124,7 @@ def forecast(self, params=None, horizon=1, start=None, align='origin', method='a
             if (params.size != np.array(self._params).size or
                     params.ndim != self._params.ndim):
                 raise ValueError('params have incorrect dimensions')
-        return self.model.forecast(params, horizon, start, align, method, simulations)
+        return self.model.forecast(params, horizon, start, align, method, simulations, rng)
 
     def hedgehog_plot(self, params=None, horizon=10, step=10, start=None,
                       type='volatility', method='analytic', simulations=1000):

arch/univariate/mean.py

@@ -629,7 +629,7 @@ def _fit_no_arch_normal_errors(self, cov_type='robust'):
                                copy.deepcopy(self))
 
     def forecast(self, parameters, horizon=1, start=None, align='origin',
-                 method='analytic', simulations=1000):
+                 method='analytic', simulations=1000, rng=None):
         # Check start
         earliest, default_start = self._fit_indices
         default_start = max(0, default_start - 1)
@@ -651,7 +651,8 @@ def forecast(self, parameters, horizon=1, start=None, align='origin',
         backcast = self._volatility.backcast(resids)
         full_resids = self.resids(mp, self._y[earliest:], self.regressors[earliest:])
         vb = self._volatility.variance_bounds(full_resids, 2.0)
-        rng = self._distribution.simulate(dp)
+        if rng is None:
+            rng = self._distribution.simulate(dp)
         variance_start = max(0, start_index - earliest)
         vfcast = self._volatility.forecast(vp, full_resids, backcast, vb,
                                            start=variance_start,

arch/univariate/volatility.py

@@ -22,7 +22,7 @@
                                                    egarch_recursion)
 
 __all__ = ['GARCH', 'ARCH', 'HARCH', 'ConstantVariance', 'EWMAVariance', 'RiskMetrics2006',
-           'EGARCH', 'FixedVariance']
+           'EGARCH', 'FixedVariance', 'BootstrapRng']
 
 
 class BootstrapRng(object):

doc/source/changes/4.0.txt

@@ -1,5 +1,7 @@
 Changes since 4.0
 =================
+- Added a parameter to forecast that allows a user-provided callable random
+  generator to be used in place of the model random generator. (:issue:`225`)
 - Added a low memory automatic lag selection method that can be used with very
   large time-series.
 - Improved performance of automatic lag selection in ADF and related tests.
@@ -15,7 +17,10 @@ Changes since 4.0
   so-called zig-zag estimation where a mean model is estimated with a fixed
   variance, and then a variance model is estimated on the residuals using
   a ``ZeroMean`` variance process.
-- Fixed a bug that prevented ``fix`` from being used with a new model (:issue:`156`)
+- 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 ability to pass optimization options to ARCH model estimation (:issue:`195`)
+- Added ability to jointly estimate smoothing parameter in EWMA variance when
+  fitting the model
+- Added ability to pass optimization options to ARCH model estimation
+  (:issue:`195`)

doc/source/univariate/univariate.rst

@@ -10,6 +10,7 @@ Univariate Volatility Models
     Examples <univariate_volatility_modeling.ipynb>
     Forecasting <forecasting>
     Forecasting Examples <univariate_volatility_forecasting.ipynb>
+    Forecasting Scenarios <univariate_volatility_scenarios.ipynb>
     Mean Models <mean>
     Volatility Processes <volatility>
     Using the Fixed Variance Process <univariate_using_fixed_variance.ipynb>