diff --git a/Forecasting.md b/Forecasting.md new file mode 100644 index 0000000000..62790c8849 --- /dev/null +++ b/Forecasting.md @@ -0,0 +1,4 @@ +1. Compute the 1-step ahead forecast in sample +2. If parameters are not common, then + 1. Update the 1-step ahead forecast out-of-sample. +3. \ No newline at end of file diff --git a/arch/univariate/base.py b/arch/univariate/base.py index f1cebb976f..b1ce16a37c 100644 --- a/arch/univariate/base.py +++ b/arch/univariate/base.py @@ -712,9 +712,9 @@ def compute_param_cov(self, params, backcast=None, robust=True): @abstractmethod def forecast(self, params, horizon=1, start=None, align='origin', method='analytic', - simulations=1000, rng=None): + simulations=1000, rng=None, insample_params=None): """ - Construct forecasts from estimated model + Construct forecasts from a model Parameters ---------- @@ -749,6 +749,11 @@ def forecast(self, params, horizon=1, start=None, align='origin', method='analyt Custom random number generator to use in simulation-based forecasts. Must produce random samples using the syntax `rng(size)` where size the 2-element tuple (simulations, horizon). + insample_params : {ndarray, Series}, optional + Parameters to use when computing in-sample residuals. If not + provided, insample_params is set to params. Setting this input + allows for alternative parameters to be used to produce forecasts + while keeping the in-sample parameters fixed. Returns ------- @@ -765,7 +770,7 @@ def forecast(self, params, horizon=1, start=None, align='origin', method='analyt >>> sim_data.index = pd.date_range('2000-01-01',periods=250) >>> am = arch_model(sim_data['data'],mean='HAR',lags=[1,5,22], vol='Constant') >>> res = am.fit() - >>> fig = res.hedgehog_plot() + >>> forecasts = res.forecast(horizon=5) Notes ----- @@ -1114,13 +1119,13 @@ def _set_tight_x(axis, index): return fig def forecast(self, params=None, horizon=1, start=None, align='origin', method='analytic', - simulations=1000, rng=None): + simulations=1000, rng=None, insample_params=None): """ - Construct forecasts from estimated model + Construct forecasts from a model Parameters ---------- - params : ndarray, optional + params : {ndarray, Series}, optional Alternative parameters to use. If not provided, the parameters estimated when fitting the model are used. Must be identical in shape to the parameters computed by fitting the model. @@ -1151,6 +1156,11 @@ def forecast(self, params=None, horizon=1, start=None, align='origin', method='a Custom random number generator to use in simulation-based forecasts. Must produce random samples using the syntax `rng(size)` where size the 2-element tuple (simulations, horizon). + insample_params : {ndarray, Series}, optional + Parameters to use when computing in-sample residuals. If not + provided, insample_params is set to params. Setting this input + allows for alternative parameters to be used to produce forecasts + while keeping the in-sample parameters fixed. Returns ------- @@ -1182,19 +1192,25 @@ def forecast(self, params=None, horizon=1, start=None, align='origin', method='a if params is None: params = self._params else: - if (params.size != np.array(self._params).size or + if (params.size != np.asarray(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, rng) + if insample_params is not None: + if (insample_params.size != np.asarray(params).size or + insample_params.ndim != params.ndim): + raise ValueError('model_params have incorrect dimensions') - def hedgehog_plot(self, params=None, horizon=10, step=10, start=None, - type='volatility', method='analytic', simulations=1000): + return self.model.forecast(params, horizon, start, align, method, simulations, rng, + insample_params) + + def hedgehog_plot(self, params=None, horizon=10, step=10, start=None, type='volatility', + method='analytic', simulations=1000, rng=None, insample_params=None): """ Plot forecasts from estimated model Parameters ---------- - params : {ndarray, Series} + params : {ndarray, Series}, optional Alternative parameters to use. If not provided, the parameters computed by fitting the model are used. Must be 1-d and identical in shape to the parameters computed by fitting the model. @@ -1219,6 +1235,15 @@ def hedgehog_plot(self, params=None, horizon=10, step=10, start=None, 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 + the 2-element tuple (simulations, horizon). + insample_params : {ndarray, Series}, optional + Parameters to use when computing in-sample residuals. If not + provided, insample_params is set to params. Setting this input + allows for alternative parameters to be used to produce forecasts + while keeping the in-sample parameters fixed. Returns ------- @@ -1244,14 +1269,16 @@ def hedgehog_plot(self, params=None, horizon=10, step=10, start=None, start = 0 while invalid_start: try: - forecasts = self.forecast(params, horizon, start, - method=method, simulations=simulations) + forecasts = self.forecast(params, horizon, start, method=method, + simulations=simulations, rng=rng, + insample_params=insample_params) invalid_start = False except ValueError: start += 1 else: forecasts = self.forecast(params, horizon, start, method=method, - simulations=simulations) + simulations=simulations, rng=rng, + insample_params=insample_params) fig, ax = plt.subplots(1, 1) use_date = isinstance(self._dep_var.index, pd.DatetimeIndex) diff --git a/arch/univariate/mean.py b/arch/univariate/mean.py index 32481503f8..78482ebaa1 100644 --- a/arch/univariate/mean.py +++ b/arch/univariate/mean.py @@ -610,7 +610,7 @@ def _fit_no_arch_normal_errors(self, cov_type='robust'): copy.deepcopy(self)) def forecast(self, params, horizon=1, start=None, align='origin', - method='analytic', simulations=1000, rng=None): + method='analytic', simulations=1000, rng=None, insample_params=None): # Check start earliest, default_start = self._fit_indices default_start = max(0, default_start - 1) @@ -623,22 +623,26 @@ def forecast(self, params, horizon=1, start=None, align='origin', # Parse params params = np.asarray(params) mp, vp, dp = self._parse_parameters(params) - + insample_mp, insample_vp, insample_db = mp, vp, dp + if insample_params is not None: + insample_mp, insample_vp, insample_db = self._parse_parameters(insample_params) ##################################### # Compute residual variance forecasts ##################################### # Back cast should use only the sample used in fitting - resids = self.resids(mp) + resids = self.resids(insample_mp) backcast = self._volatility.backcast(resids) - full_resids = self.resids(mp, self._y[earliest:], self.regressors[earliest:]) + full_resids = self.resids(insample_mp, self._y[earliest:], self.regressors[earliest:]) vb = self._volatility.variance_bounds(full_resids, 2.0) if rng is None: rng = self._distribution.simulate(dp) variance_start = max(0, start_index - earliest) + # TODO: Adapt for in-sample parameters vfcast = self._volatility.forecast(vp, full_resids, backcast, vb, start=variance_start, horizon=horizon, method=method, - simulations=simulations, rng=rng) + simulations=simulations, rng=rng, + insample_params=insample_vp) var_fcasts = vfcast.forecasts var_fcasts = _forecast_pad(earliest, var_fcasts) diff --git a/arch/univariate/volatility.py b/arch/univariate/volatility.py index 6117ecfb41..33888ddb04 100644 --- a/arch/univariate/volatility.py +++ b/arch/univariate/volatility.py @@ -192,7 +192,8 @@ def _check_forecasting_method(self, method, horizon): """ pass - def _one_step_forecast(self, parameters, resids, backcast, var_bounds, horizon): + def _one_step_forecast(self, parameters, insample_parameters, resids, backcast, var_bounds, + horizon, start): """ One-step ahead forecast @@ -200,6 +201,10 @@ def _one_step_forecast(self, parameters, resids, backcast, var_bounds, horizon): ---------- parameters : ndarray Parameters required to forecast the volatility model + insample_parameters ndarray + Parameters used to compute in-sample volatility. Set equal to + parameters to use the same values to fit the model and to + forecast from the model. resids : ndarray Residuals to use in the recursion backcast : float @@ -209,6 +214,8 @@ def _one_step_forecast(self, parameters, resids, backcast, var_bounds, horizon): horizon : int Forecast horizon. Must be 1 or larger. Forecasts are produced for horizons in [1, horizon]. + start : int + Observation index of first out-of-sample forecast Returns ------- @@ -221,15 +228,43 @@ def _one_step_forecast(self, parameters, resids, backcast, var_bounds, horizon): _resids = np.concatenate((resids, [0])) _var_bounds = np.concatenate((var_bounds, [[0, np.inf]])) sigma2 = np.zeros(t + 1) - self.compute_variance(parameters, _resids, sigma2, backcast, _var_bounds) - forecasts = np.zeros((t, horizon)) - forecasts[:, 0] = sigma2[1:] + self.compute_variance(insample_parameters, _resids, sigma2, backcast, _var_bounds) + forecasts = np.full((t, horizon), np.nan) + forecasts[start:, 0] = sigma2[start + 1:] sigma2 = sigma2[:-1] - + self._update_one_step_forecast(parameters, resids, backcast, sigma2, forecasts, start) return sigma2, forecasts @abstractmethod - def _analytic_forecast(self, parameters, resids, backcast, var_bounds, start, horizon): + def _update_one_step_forecast(self, parameters, resids, backcast, sigma2, forecasts, start): + """ + Update the out-of-sample forecasts using the correct parameters + + Parameters + ---------- + parameters : ndarray + Parameters required to forecast the volatility model + resids : ndarray + Residuals to use in the recursion + backcast : float + Value to use when initializing the recursion + sigma2 : ndarray + t element array containing the one-step ahead forecasts + forecasts : ndarray + t by horizon array containing the one-step ahead forecasts in the first location + start : int + Observation index of first out-of-sample forecast + + Returns + ------- + forecasts : ndarray + t by horizon array containing the updates one-step ahead forecasts in column 0 + """ + pass + + @abstractmethod + def _analytic_forecast(self, parameters, insample_parameters, resids, backcast, var_bounds, + start, horizon): """ Analytic multi-step volatility forecasts from the model @@ -237,6 +272,10 @@ def _analytic_forecast(self, parameters, resids, backcast, var_bounds, start, ho ---------- parameters : ndarray Parameters required to forecast the volatility model + insample_parameters ndarray + Parameters used to compute in-sample volatility. Set equal to + parameters to use the same values to fit the model and to + forecast from the model. resids : ndarray Residuals to use in the recursion backcast : float @@ -259,8 +298,8 @@ def _analytic_forecast(self, parameters, resids, backcast, var_bounds, start, ho pass @abstractmethod - def _simulation_forecast(self, parameters, resids, backcast, var_bounds, start, horizon, - simulations, rng): + def _simulation_forecast(self, parameters, insample_parameters, resids, backcast, var_bounds, + start, horizon, simulations, rng): """ Simulation-based volatility forecasts from the model @@ -268,6 +307,10 @@ def _simulation_forecast(self, parameters, resids, backcast, var_bounds, start, ---------- parameters : ndarray Parameters required to forecast the volatility model + insample_parameters ndarray + Parameters used to compute in-sample volatility. Set equal to + parameters to use the same values to fit the model and to + forecast from the model. resids : ndarray Residuals to use in the recursion backcast : float @@ -296,8 +339,8 @@ def _simulation_forecast(self, parameters, resids, backcast, var_bounds, start, """ pass - def _bootstrap_forecast(self, parameters, resids, backcast, var_bounds, start, horizon, - simulations, random_state): + def _bootstrap_forecast(self, parameters, insample_parameters, resids, backcast, var_bounds, + start, horizon, simulations, random_state): """ Simulation-based volatility forecasts using model residuals @@ -305,6 +348,10 @@ def _bootstrap_forecast(self, parameters, resids, backcast, var_bounds, start, h ---------- parameters : ndarray Parameters required to forecast the volatility model + insample_parameters : ndarray + Parameters used to compute in-sample volatility. Set equal to + parameters to use the same values to fit the model and to + forecast from the model. resids : ndarray Residuals to use in the recursion backcast : float @@ -336,8 +383,8 @@ def _bootstrap_forecast(self, parameters, resids, backcast, var_bounds, start, h raise ValueError('start must include more than {0} ' 'observations'.format(self._min_bootstrap_obs)) rng = BootstrapRng(std_resid, start, random_state=random_state).rng() - return self._simulation_forecast(parameters, resids, backcast, var_bounds, - start, horizon, simulations, rng) + return self._simulation_forecast(parameters, insample_parameters, resids, backcast, + var_bounds, start, horizon, simulations, rng) def variance_bounds(self, resids, power=2.0): """ @@ -494,7 +541,8 @@ def constraints(self): pass def forecast(self, parameters, resids, backcast, var_bounds, start=None, horizon=1, - method='analytic', simulations=1000, rng=None, random_state=None): + method='analytic', simulations=1000, rng=None, random_state=None, + insample_parameters=None): """ Forecast volatility from the model @@ -525,6 +573,9 @@ def forecast(self, parameters, resids, backcast, var_bounds, start=None, horizon samples numbers with that shape. random_state : RandomState, optional NumPy RandomState instance to use when method is 'bootstrap' + insample_parameters : ndarray, optional + Parameters to use for in-sample volatility fitting. It not + provided, insample_parameters is set to parameters. Returns ------- @@ -545,6 +596,9 @@ def forecast(self, parameters, resids, backcast, var_bounds, start=None, horizon to use this method when not available will raise a ValueError. """ parameters = np.asarray(parameters) + if insample_parameters is None: + insample_parameters = parameters + insample_parameters = np.asarray(insample_parameters) method = method.lower() if method not in ('analytic', 'simulation', 'bootstrap'): raise ValueError('{0} is not a known forecasting method'.format(method)) @@ -553,18 +607,18 @@ def forecast(self, parameters, resids, backcast, var_bounds, start=None, horizon start = len(resids) - 1 if start is None else start if method == 'analytic': - return self._analytic_forecast(parameters, resids, backcast, var_bounds, start, - horizon) + return self._analytic_forecast(parameters, insample_parameters, resids, backcast, + var_bounds, start, horizon) elif method == 'simulation': - return self._simulation_forecast(parameters, resids, backcast, var_bounds, start, - horizon, simulations, rng) + return self._simulation_forecast(parameters, insample_parameters, resids, backcast, + var_bounds, start, horizon, simulations, rng) else: if start < 10 or (horizon / start) >= .2: raise ValueError('Bootstrap forecasting requires at least 10 initial ' 'observations, and the ratio of horizon-to-start < 20%.') - return self._bootstrap_forecast(parameters, resids, backcast, var_bounds, start, - horizon, simulations, random_state) + return self._bootstrap_forecast(parameters, insample_parameters, resids, backcast, + var_bounds, start, horizon, simulations, random_state) @abstractmethod def simulate(self, parameters, nobs, rng, burn=500, initial_value=None): @@ -667,7 +721,11 @@ def parameter_names(self): def _check_forecasting_method(self, method, horizon): return - def _analytic_forecast(self, parameters, resids, backcast, var_bounds, start, horizon): + def _update_one_step_forecast(self, parameters, resids, backcast, sigma2, forecasts, start): + return forecasts + + def _analytic_forecast(self, parameters, insample_parameters, resids, backcast, var_bounds, + start, horizon): t = resids.shape[0] forecasts = np.full((t, horizon), np.nan) @@ -675,8 +733,8 @@ def _analytic_forecast(self, parameters, resids, backcast, var_bounds, start, ho forecast_paths = None return VarianceForecast(forecasts, forecast_paths) - def _simulation_forecast(self, parameters, resids, backcast, var_bounds, start, horizon, - simulations, rng): + def _simulation_forecast(self, parameters, insample_parameters, resids, backcast, var_bounds, + start, horizon, simulations, rng): t = resids.shape[0] forecasts = np.full((t, horizon), np.nan) forecast_paths = np.full((t, simulations, horizon), np.nan) @@ -965,12 +1023,37 @@ def _check_forecasting_method(self, method, horizon): raise ValueError('Analytic forecasts not available for horizon > 1 when power != 2') return - def _analytic_forecast(self, parameters, resids, backcast, var_bounds, start, horizon): + def _update_one_step_forecast(self, parameters, resids, backcast, sigma2, forecasts, start): + p, o, q = self.p, self.o, self.q + t = resids.shape[0] - sigma2, forecasts = self._one_step_forecast(parameters, resids, backcast, - var_bounds, horizon) + sresids = resids < 0 + power = self.power + power_resids = np.abs(resids) ** power + half_power = power / 2.0 + for i in range(start, t): + # 1-step ahead using time i information + fcast = parameters[0] + for j in range(p): + shock = power_resids[i - j] if (i - j) >= 0 else backcast ** half_power + fcast += parameters[1 + j] * shock + for j in range(o): + if (i - j) >= 0: + shock = power_resids[i - j] * sresids[i - j] + else: + shock = 0.5 * backcast ** half_power + fcast += parameters[1 + p + j] * shock + for j in range(q): + lag = sigma2 ** half_power if (i - j) >= 0 else backcast ** half_power + fcast += parameters[1 + p + o + j] * lag + forecasts[i, 0] = fcast ** (2.0 / power) + return forecasts + + def _analytic_forecast(self, parameters, insample_parameters, resids, backcast, var_bounds, start, horizon): + + sigma2, forecasts = self._one_step_forecast(parameters, insample_parameters, resids, + backcast, var_bounds, horizon, start) if horizon == 1: - forecasts[:start] = np.nan return VarianceForecast(forecasts) t = resids.shape[0] @@ -1051,11 +1134,11 @@ def _simulate_paths(self, m, parameters, horizon, std_shocks, return np.mean(forecast_paths, 0), forecast_paths, shock[:, m:] - def _simulation_forecast(self, parameters, resids, backcast, var_bounds, start, horizon, - simulations, rng): + def _simulation_forecast(self, parameters, insample_parameters, resids, backcast, var_bounds, + start, horizon, simulations, rng): - sigma2, forecasts = self._one_step_forecast(parameters, resids, backcast, - var_bounds, horizon) + sigma2, forecasts = self._one_step_forecast(parameters, insample_parameters, resids, + backcast, var_bounds, horizon, start) t = resids.shape[0] paths = np.full((t, simulations, horizon), np.nan) shocks = np.full((t, simulations, horizon), np.nan)