From 9fd1e724ca4e94906b3cb02a4205b1d81c33250b Mon Sep 17 00:00:00 2001 From: Andrew Scott <84710236+andrewscottm@users.noreply.github.com> Date: Fri, 22 Nov 2024 19:07:17 +0000 Subject: [PATCH] fix: seasonal naive confidence interval bug (#932) --- nbs/src/core/models.ipynb | 24 ++++++++++-- python/statsforecast/models.py | 69 +++++++++++++++++----------------- 2 files changed, 55 insertions(+), 38 deletions(-) diff --git a/nbs/src/core/models.ipynb b/nbs/src/core/models.ipynb index 947cd11a4..42f3ddb4f 100644 --- a/nbs/src/core/models.ipynb +++ b/nbs/src/core/models.ipynb @@ -6212,7 +6212,7 @@ " if self.prediction_intervals is not None:\n", " res = self._add_predict_conformal_intervals(res, level)\n", " else:\n", - " k = np.floor((h - 1) / self.season_length)\n", + " k = np.floor(np.arange(h) / self.season_length)\n", " sigma = self.model_['sigma']\n", " sigmah = sigma * np.sqrt(k + 1)\n", " pred_int = _calculate_intervals(res, level, h, sigmah)\n", @@ -6287,14 +6287,13 @@ " if self.prediction_intervals is not None:\n", " res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)\n", " else:\n", - " k = np.floor((h - 1) / self.season_length)\n", + " k = np.floor(np.arange(h) / self.season_length)\n", " residuals = y - out[\"fitted\"]\n", " sigma = _calculate_sigma(residuals, len(y) - self.season_length)\n", " sigmah = sigma * np.sqrt(k + 1)\n", " pred_int = _calculate_intervals(out, level, h, sigmah)\n", " res = {**res, **pred_int}\n", " if fitted:\n", - " k = np.floor((h - 1) / self.season_length)\n", " residuals = y - out[\"fitted\"]\n", " sigma = _calculate_sigma(residuals, len(y) - self.season_length)\n", " res = _add_fitted_pi(res=res, se=sigma, level=level)\n", @@ -6360,6 +6359,25 @@ "_plot_insample_pi(fcst_seas_naive)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#| hide\n", + "# test h > season_length\n", + "seas_naive_bigh = SeasonalNaive(season_length=12)\n", + "fcst_seas_naive_bigh = seas_naive_bigh.forecast(ap, 13, None, None, (80,95), True)\n", + "np.testing.assert_almost_equal(\n", + " fcst_seas_naive_bigh['lo-80'][:12],\n", + " np.array([370.4595, 344.4595, 372.4595, 414.4595, 425.4595, 488.4595, \n", + " 575.4595, 559.4595, 461.4595, 414.4595, 343.4595, 385.4595]),\n", + " decimal=4\n", + ")\n", + "_plot_fcst(fcst_seas_naive_bigh)" + ] + }, { "cell_type": "code", "execution_count": null, diff --git a/python/statsforecast/models.py b/python/statsforecast/models.py index dad5285ce..9f7cf0994 100644 --- a/python/statsforecast/models.py +++ b/python/statsforecast/models.py @@ -3334,7 +3334,7 @@ def predict( if self.prediction_intervals is not None: res = self._add_predict_conformal_intervals(res, level) else: - k = np.floor((h - 1) / self.season_length) + k = np.floor(np.arange(h) / self.season_length) sigma = self.model_["sigma"] sigmah = sigma * np.sqrt(k + 1) pred_int = _calculate_intervals(res, level, h, sigmah) @@ -3410,20 +3410,19 @@ def forecast( if self.prediction_intervals is not None: res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level) else: - k = np.floor((h - 1) / self.season_length) + k = np.floor(np.arange(h) / self.season_length) residuals = y - out["fitted"] sigma = _calculate_sigma(residuals, len(y) - self.season_length) sigmah = sigma * np.sqrt(k + 1) pred_int = _calculate_intervals(out, level, h, sigmah) res = {**res, **pred_int} if fitted: - k = np.floor((h - 1) / self.season_length) residuals = y - out["fitted"] sigma = _calculate_sigma(residuals, len(y) - self.season_length) res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../../nbs/src/core/models.ipynb 262 +# %% ../../nbs/src/core/models.ipynb 263 def _window_average( y: np.ndarray, # time series h: int, # forecasting horizon @@ -3438,7 +3437,7 @@ def _window_average( mean = _repeat_val(val=wavg, h=h) return {"mean": mean} -# %% ../../nbs/src/core/models.ipynb 263 +# %% ../../nbs/src/core/models.ipynb 264 class WindowAverage(_TS): def __init__( @@ -3596,7 +3595,7 @@ def forecast( raise Exception("You must pass `prediction_intervals` to compute them.") return res -# %% ../../nbs/src/core/models.ipynb 274 +# %% ../../nbs/src/core/models.ipynb 275 def _seasonal_window_average( y: np.ndarray, h: int, @@ -3613,7 +3612,7 @@ def _seasonal_window_average( out = _repeat_val_seas(season_vals=season_avgs, h=h) return {"mean": out} -# %% ../../nbs/src/core/models.ipynb 275 +# %% ../../nbs/src/core/models.ipynb 276 class SeasonalWindowAverage(_TS): def __init__( @@ -3785,7 +3784,7 @@ def forecast( raise Exception("You must pass `prediction_intervals` to compute them.") return res -# %% ../../nbs/src/core/models.ipynb 287 +# %% ../../nbs/src/core/models.ipynb 288 def _chunk_forecast(y, aggregation_level): lost_remainder_data = len(y) % aggregation_level y_cut = y[lost_remainder_data:] @@ -3868,7 +3867,7 @@ def _adida( res["fitted"] = np.append(np.nan, sums_fitted / fitted_aggregation_levels) return res -# %% ../../nbs/src/core/models.ipynb 288 +# %% ../../nbs/src/core/models.ipynb 289 class ADIDA(_TS): def __init__( @@ -4037,7 +4036,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../../nbs/src/core/models.ipynb 300 +# %% ../../nbs/src/core/models.ipynb 301 def _croston_classic( y: np.ndarray, # time series h: int, # forecasting horizon @@ -4065,7 +4064,7 @@ def _croston_classic( out["fitted"] = ydf / yif return out -# %% ../../nbs/src/core/models.ipynb 301 +# %% ../../nbs/src/core/models.ipynb 302 class CrostonClassic(_TS): def __init__( @@ -4229,7 +4228,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../../nbs/src/core/models.ipynb 312 +# %% ../../nbs/src/core/models.ipynb 313 def _croston_optimized( y: np.ndarray, # time series h: int, # forecasting horizon @@ -4271,7 +4270,7 @@ def _croston_optimized( out["fitted"] = ydf / yif return out -# %% ../../nbs/src/core/models.ipynb 313 +# %% ../../nbs/src/core/models.ipynb 314 class CrostonOptimized(_TS): def __init__( @@ -4432,7 +4431,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../../nbs/src/core/models.ipynb 324 +# %% ../../nbs/src/core/models.ipynb 325 def _croston_sba( y: np.ndarray, # time series h: int, # forecasting horizon @@ -4444,7 +4443,7 @@ def _croston_sba( out["fitted"] *= 0.95 return out -# %% ../../nbs/src/core/models.ipynb 325 +# %% ../../nbs/src/core/models.ipynb 326 class CrostonSBA(_TS): def __init__( @@ -4610,7 +4609,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../../nbs/src/core/models.ipynb 336 +# %% ../../nbs/src/core/models.ipynb 337 def _imapa( y: np.ndarray, # time series h: int, # forecasting horizon @@ -4644,7 +4643,7 @@ def _imapa( res["fitted"] = fitted_vals return res -# %% ../../nbs/src/core/models.ipynb 337 +# %% ../../nbs/src/core/models.ipynb 338 class IMAPA(_TS): def __init__( @@ -4809,7 +4808,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../../nbs/src/core/models.ipynb 348 +# %% ../../nbs/src/core/models.ipynb 349 def _tsb( y: np.ndarray, # time series h: int, # forecasting horizon @@ -4834,7 +4833,7 @@ def _tsb( res["fitted"] = ypft * ydft return res -# %% ../../nbs/src/core/models.ipynb 349 +# %% ../../nbs/src/core/models.ipynb 350 class TSB(_TS): def __init__( @@ -5009,7 +5008,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../../nbs/src/core/models.ipynb 361 +# %% ../../nbs/src/core/models.ipynb 362 def _predict_mstl_components(mstl_ob, h, season_length): seasoncolumns = mstl_ob.filter(regex="seasonal*").columns nseasons = len(seasoncolumns) @@ -5030,7 +5029,7 @@ def _predict_mstl_seas(mstl_ob, h, season_length): seascomp = _predict_mstl_components(mstl_ob, h, season_length) return seascomp.sum(axis=1) -# %% ../../nbs/src/core/models.ipynb 362 +# %% ../../nbs/src/core/models.ipynb 363 class MSTL(_TS): r"""MSTL model. @@ -5309,7 +5308,7 @@ def forward( } return res -# %% ../../nbs/src/core/models.ipynb 378 +# %% ../../nbs/src/core/models.ipynb 379 class TBATS(_TS): r"""Trigonometric Box-Cox transform, ARMA errors, Trend and Seasonal components (TBATS) model. @@ -5523,7 +5522,7 @@ def forecast( res_trans = res return res_trans -# %% ../../nbs/src/core/models.ipynb 386 +# %% ../../nbs/src/core/models.ipynb 387 class AutoTBATS(TBATS): r"""AutoTBATS model. @@ -5582,7 +5581,7 @@ def __init__( alias=alias, ) -# %% ../../nbs/src/core/models.ipynb 396 +# %% ../../nbs/src/core/models.ipynb 397 class Theta(AutoTheta): r"""Standard Theta Method. @@ -5619,7 +5618,7 @@ def __init__( prediction_intervals=prediction_intervals, ) -# %% ../../nbs/src/core/models.ipynb 410 +# %% ../../nbs/src/core/models.ipynb 411 class OptimizedTheta(AutoTheta): r"""Optimized Theta Method. @@ -5656,7 +5655,7 @@ def __init__( prediction_intervals=prediction_intervals, ) -# %% ../../nbs/src/core/models.ipynb 424 +# %% ../../nbs/src/core/models.ipynb 425 class DynamicTheta(AutoTheta): r"""Dynamic Standard Theta Method. @@ -5693,7 +5692,7 @@ def __init__( prediction_intervals=prediction_intervals, ) -# %% ../../nbs/src/core/models.ipynb 438 +# %% ../../nbs/src/core/models.ipynb 439 class DynamicOptimizedTheta(AutoTheta): r"""Dynamic Optimized Theta Method. @@ -5730,7 +5729,7 @@ def __init__( prediction_intervals=prediction_intervals, ) -# %% ../../nbs/src/core/models.ipynb 453 +# %% ../../nbs/src/core/models.ipynb 454 class GARCH(_TS): r"""Generalized Autoregressive Conditional Heteroskedasticity (GARCH) model. @@ -5924,7 +5923,7 @@ def forecast( res = _add_fitted_pi(res=res, se=se, level=level) return res -# %% ../../nbs/src/core/models.ipynb 466 +# %% ../../nbs/src/core/models.ipynb 467 class ARCH(GARCH): r"""Autoregressive Conditional Heteroskedasticity (ARCH) model. @@ -5968,7 +5967,7 @@ def __init__( self.alias = alias super().__init__(p, q=0, alias=alias) -# %% ../../nbs/src/core/models.ipynb 477 +# %% ../../nbs/src/core/models.ipynb 478 class SklearnModel(_TS): r"""scikit-learn model wrapper @@ -6176,7 +6175,7 @@ def forward( res = _add_fitted_pi(res=res, se=se, level=level) return res -# %% ../../nbs/src/core/models.ipynb 487 +# %% ../../nbs/src/core/models.ipynb 488 class MFLES(_TS): r"""MFLES model. @@ -6458,7 +6457,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../../nbs/src/core/models.ipynb 495 +# %% ../../nbs/src/core/models.ipynb 496 class AutoMFLES(_TS): r"""AutoMFLES @@ -6659,7 +6658,7 @@ def forecast( res = _add_fitted_pi(res=res, se=sigma, level=level) return res -# %% ../../nbs/src/core/models.ipynb 499 +# %% ../../nbs/src/core/models.ipynb 500 class ConstantModel(_TS): def __init__(self, constant: float, alias: str = "ConstantModel"): @@ -6845,7 +6844,7 @@ def forward( ) return res -# %% ../../nbs/src/core/models.ipynb 513 +# %% ../../nbs/src/core/models.ipynb 514 class ZeroModel(ConstantModel): def __init__(self, alias: str = "ZeroModel"): @@ -6860,7 +6859,7 @@ def __init__(self, alias: str = "ZeroModel"): """ super().__init__(constant=0, alias=alias) -# %% ../../nbs/src/core/models.ipynb 527 +# %% ../../nbs/src/core/models.ipynb 528 class NaNModel(ConstantModel): def __init__(self, alias: str = "NaNModel"):