diff --git a/arch/tests/univariate/test_variance_forecasting.py b/arch/tests/univariate/test_variance_forecasting.py index 7fb5ab9b41..65c4af2709 100644 --- a/arch/tests/univariate/test_variance_forecasting.py +++ b/arch/tests/univariate/test_variance_forecasting.py @@ -2,20 +2,22 @@ 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.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 +32,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 +310,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 +602,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 +838,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 +937,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 +1255,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 +1292,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 +1472,26 @@ def test_bs_rng_errors(self): y = self.rng.rand(1000) with pytest.raises(ValueError): BootstrapRng(y, 0) + + +def test_external_rng(): + from arch.univariate import ConstantMean, GARCH, Normal + + am = arch_model(None, mean='Constant', vol='GARCH', p=1, q=1) + data = am.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]) + + rs = np.random.RandomState(1) + state = rs.get_state() + volatility = GARCH(1, 0, 1) + distribution = Normal(rs) + mod = ConstantMean(data.data, volatility=volatility, distribution=distribution) + res = mod.fit(disp='off') + + rs.set_state(state) + f1 = res.forecast(res.params, horizon=5, method='simulation', start=900, simulations=250) + rs.set_state(state) + rng = rs.standard_normal + f2 = res.forecast(res.params, horizon=5, method='simulation', start=900, simulations=250, + rng=rng) + assert_allclose(f1.residual_variance, f2.residual_variance) diff --git a/arch/univariate/base.py b/arch/univariate/base.py index 208db07a05..182deafefb 100644 --- a/arch/univariate/base.py +++ b/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): diff --git a/arch/univariate/mean.py b/arch/univariate/mean.py index 56b484cc99..585f190010 100644 --- a/arch/univariate/mean.py +++ b/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, diff --git a/doc/source/changes/4.0.txt b/doc/source/changes/4.0.txt index 8bc9eacb9e..547e1d9539 100644 --- a/doc/source/changes/4.0.txt +++ b/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`)