models.py

# AUTOGENERATED! DO NOT EDIT! File to edit: ../nbs/src/core/models.ipynb.

# %% auto 0
__all__ = ['AutoARIMA', 'AutoETS', 'ETS', 'AutoCES', 'AutoTheta', 'ARIMA', 'AutoRegressive', 'SimpleExponentialSmoothing',
           'SimpleExponentialSmoothingOptimized', 'SeasonalExponentialSmoothing',
           'SeasonalExponentialSmoothingOptimized', 'Holt', 'HoltWinters', 'HistoricAverage', 'Naive',
           'RandomWalkWithDrift', 'SeasonalNaive', 'WindowAverage', 'SeasonalWindowAverage', 'ADIDA', 'CrostonClassic',
           'CrostonOptimized', 'CrostonSBA', 'IMAPA', 'TSB', 'MSTL', 'TBATS', 'AutoTBATS', 'Theta', 'OptimizedTheta',
           'DynamicTheta', 'DynamicOptimizedTheta', 'GARCH', 'ARCH', 'SklearnModel', 'ConstantModel', 'ZeroModel',
           'NaNModel']

# %% ../nbs/src/core/models.ipynb 5
import warnings
from math import trunc
from typing import Any, Dict, List, Optional, Sequence, Tuple, Union

import numpy as np
from numba import njit
from scipy.optimize import minimize
from scipy.special import inv_boxcox

from statsforecast.arima import (
    Arima,
    auto_arima_f,
    fitted_arima,
    forecast_arima,
    forward_arima,
    is_constant,
)
from .ces import auto_ces, forecast_ces, forward_ces
from statsforecast.ets import (
    _PHI_LOWER,
    _PHI_UPPER,
    ets_f,
    forecast_ets,
    forward_ets,
)
from .mstl import mstl
from .theta import auto_theta, forecast_theta, forward_theta
from .garch import garch_model, garch_forecast
from .tbats import tbats_selection, tbats_forecast, _compute_sigmah
from statsforecast.utils import (
    _calculate_sigma,
    _calculate_intervals,
    _ensure_float,
    _naive,
    _old_kw_to_pos,
    _quantiles,
    _repeat_val,
    _repeat_val_seas,
    _seasonal_naive,
    CACHE,
    ConformalIntervals,
    NOGIL,
)

# %% ../nbs/src/core/models.ipynb 9
def _add_fitted_pi(res, se, level):
    level = sorted(level)
    level = np.asarray(level)
    quantiles = _quantiles(level=level)
    lo = res["fitted"].reshape(-1, 1) - quantiles * se.reshape(-1, 1)
    hi = res["fitted"].reshape(-1, 1) + quantiles * se.reshape(-1, 1)
    lo = lo[:, ::-1]
    lo = {f"fitted-lo-{l}": lo[:, i] for i, l in enumerate(reversed(level))}
    hi = {f"fitted-hi-{l}": hi[:, i] for i, l in enumerate(level)}
    res = {**res, **lo, **hi}
    return res

# %% ../nbs/src/core/models.ipynb 10
def _add_conformal_distribution_intervals(
    fcst: Dict,
    cs: np.ndarray,
    level: List[Union[int, float]],
) -> Dict:
    """
    Adds conformal intervals to the `fcst` dict based on conformal scores `cs`.
    `level` should be already sorted. This strategy creates forecasts paths
    based on errors and calculate quantiles using those paths.
    """
    alphas = [100 - lv for lv in level]
    cuts = [alpha / 200 for alpha in reversed(alphas)]
    cuts.extend(1 - alpha / 200 for alpha in alphas)
    mean = fcst["mean"].reshape(1, -1)
    scores = np.vstack([mean - cs, mean + cs])
    quantiles = np.quantile(
        scores,
        cuts,
        axis=0,
    )
    quantiles = quantiles.reshape(len(cuts), -1)
    lo_cols = [f"lo-{lv}" for lv in reversed(level)]
    hi_cols = [f"hi-{lv}" for lv in level]
    out_cols = lo_cols + hi_cols
    for i, col in enumerate(out_cols):
        fcst[col] = quantiles[i]
    return fcst

# %% ../nbs/src/core/models.ipynb 11
def _get_conformal_method(method: str):
    available_methods = {
        "conformal_distribution": _add_conformal_distribution_intervals,
        # "conformal_error": _add_conformal_error_intervals,
    }
    if method not in available_methods.keys():
        raise ValueError(
            f"prediction intervals method {method} not supported "
            f'please choose one of {", ".join(available_methods.keys())}'
        )
    return available_methods[method]

# %% ../nbs/src/core/models.ipynb 12
class _TS:
    uses_exog = False

    def new(self):
        b = type(self).__new__(type(self))
        b.__dict__.update(self.__dict__)
        return b

    def _conformity_scores(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ) -> np.ndarray:
        n_windows = self.prediction_intervals.n_windows  # type: ignore[attr-defined]
        h = self.prediction_intervals.h  # type: ignore[attr-defined]
        n_samples = y.size
        # use as many windows as possible for short series
        # subtract 1 for the training set
        n_windows = min(n_windows, (n_samples - 1) // h)
        if n_windows < 2:
            raise ValueError(
                f"Prediction intervals settings require at least {2 * h + 1:,} samples, serie has {n_samples:,}."
            )
        test_size = n_windows * h
        cs = np.empty((n_windows, h), dtype=np.float32)
        for i_window in range(n_windows):
            train_end = n_samples - test_size + i_window * h
            y_train = y[:train_end]
            y_test = y[train_end : train_end + h]
            if X is not None:
                X_train = X[:train_end]
                X_test = X[train_end : train_end + h]
            else:
                X_train = None
                X_test = None
            fcst_window = self.forecast(h=h, y=y_train, X=X_train, X_future=X_test)  # type: ignore[attr-defined]
            cs[i_window] = np.abs(fcst_window["mean"] - y_test)
        return cs

    @property
    def _conformal_method(self):
        return _get_conformal_method(self.prediction_intervals.method)

    def _store_cs(self, y, X):
        if self.prediction_intervals is not None:
            self._cs = self._conformity_scores(y, X)

    def _add_conformal_intervals(self, fcst, y, X, level):
        if self.prediction_intervals is not None and level is not None:
            cs = self._conformity_scores(y, X) if y is not None else self._cs
            res = self._conformal_method(fcst=fcst, cs=cs, level=level)
            return res
        return fcst

    def _add_predict_conformal_intervals(self, fcst, level):
        return self._add_conformal_intervals(fcst=fcst, y=None, X=None, level=level)

# %% ../nbs/src/core/models.ipynb 17
class AutoARIMA(_TS):
    """AutoARIMA model.

    Automatically selects the best ARIMA (AutoRegressive Integrated Moving Average)
    model using an information criterion. Default is Akaike Information Criterion (AICc).

    Notes
    -----
    This implementation is a mirror of Hyndman's [forecast::auto.arima](https://github.com/robjhyndman/forecast).

    References
    ----------
    [Rob J. Hyndman, Yeasmin Khandakar (2008). "Automatic Time Series Forecasting: The forecast package for R"](https://www.jstatsoft.org/article/view/v027i03).

    Parameters
    ----------
    d : Optional[int]
        Order of first-differencing.
    D : Optional[int]
        Order of seasonal-differencing.
    max_p : int
        Max autorregresives p.
    max_q : int
        Max moving averages q.
    max_P : int
        Max seasonal autorregresives P.
    max_Q : int
        Max seasonal moving averages Q.
    max_order : int
        Max p+q+P+Q value if not stepwise selection.
    max_d : int
        Max non-seasonal differences.
    max_D : int
        Max seasonal differences.
    start_p : int
        Starting value of p in stepwise procedure.
    start_q : int
        Starting value of q in stepwise procedure.
    start_P : int
        Starting value of P in stepwise procedure.
    start_Q : int
        Starting value of Q in stepwise procedure.
    stationary : bool
        If True, restricts search to stationary models.
    seasonal : bool
        If False, restricts search to non-seasonal models.
    ic : str
        Information criterion to be used in model selection.
    stepwise : bool
        If True, will do stepwise selection (faster).
    nmodels : int
        Number of models considered in stepwise search.
    trace : bool
        If True, the searched ARIMA models is reported.
    approximation : Optional[bool]
        If True, conditional sums-of-squares estimation, final MLE.
    method : Optional[str]
        Fitting method between maximum likelihood or sums-of-squares.
    truncate : Optional[int]
        Observations truncated series used in model selection.
    test : str
        Unit root test to use. See `ndiffs` for details.
    test_kwargs : Optional[str]
        Unit root test additional arguments.
    seasonal_test : str
        Selection method for seasonal differences.
    seasonal_test_kwargs : Optional[dict]
        Seasonal unit root test arguments.
    allowdrift : bool (default True)
        If True, drift models terms considered.
    allowmean : bool (default True)
        If True, non-zero mean models considered.
    blambda : Optional[float]
        Box-Cox transformation parameter.
    biasadj : bool
        Use adjusted back-transformed mean Box-Cox.
    season_length : int
        Number of observations per unit of time. Ex: 24 Hourly data.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    uses_exog = True

    def __init__(
        self,
        d: Optional[int] = None,
        D: Optional[int] = None,
        max_p: int = 5,
        max_q: int = 5,
        max_P: int = 2,
        max_Q: int = 2,
        max_order: int = 5,
        max_d: int = 2,
        max_D: int = 1,
        start_p: int = 2,
        start_q: int = 2,
        start_P: int = 1,
        start_Q: int = 1,
        stationary: bool = False,
        seasonal: bool = True,
        ic: str = "aicc",
        stepwise: bool = True,
        nmodels: int = 94,
        trace: bool = False,
        approximation: Optional[bool] = False,
        method: Optional[str] = None,
        truncate: Optional[bool] = None,
        test: str = "kpss",
        test_kwargs: Optional[str] = None,
        seasonal_test: str = "seas",
        seasonal_test_kwargs: Optional[Dict] = None,
        allowdrift: bool = False,
        allowmean: bool = False,
        blambda: Optional[float] = None,
        biasadj: bool = False,
        season_length: int = 1,
        alias: str = "AutoARIMA",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.d = d
        self.D = D
        self.max_p = max_p
        self.max_q = max_q
        self.max_P = max_P
        self.max_Q = max_Q
        self.max_order = max_order
        self.max_d = max_d
        self.max_D = max_D
        self.start_p = start_p
        self.start_q = start_q
        self.start_P = start_P
        self.start_Q = start_Q
        self.stationary = stationary
        self.seasonal = seasonal
        self.ic = ic
        self.stepwise = stepwise
        self.nmodels = nmodels
        self.trace = trace
        self.approximation = approximation
        self.method = method
        self.truncate = truncate
        self.test = test
        self.test_kwargs = test_kwargs
        self.seasonal_test = seasonal_test
        self.seasonal_test_kwargs = seasonal_test_kwargs
        self.allowdrift = allowdrift
        self.allowmean = allowmean
        self.blambda = blambda
        self.biasadj = biasadj
        self.season_length = season_length
        self.alias = alias
        self.prediction_intervals = prediction_intervals

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the AutoARIMA model.

        Fit an AutoARIMA to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            AutoARIMA fitted model.
        """
        with np.errstate(invalid="ignore"):
            self.model_ = auto_arima_f(
                x=y,
                d=self.d,
                D=self.D,
                max_p=self.max_p,
                max_q=self.max_q,
                max_P=self.max_P,
                max_Q=self.max_Q,
                max_order=self.max_order,
                max_d=self.max_d,
                max_D=self.max_D,
                start_p=self.start_p,
                start_q=self.start_q,
                start_P=self.start_P,
                start_Q=self.start_Q,
                stationary=self.stationary,
                seasonal=self.seasonal,
                ic=self.ic,
                stepwise=self.stepwise,
                nmodels=self.nmodels,
                trace=self.trace,
                approximation=self.approximation,
                method=self.method,
                truncate=self.truncate,
                xreg=X,
                test=self.test,
                test_kwargs=self.test_kwargs,
                seasonal_test=self.seasonal_test,
                seasonal_test_kwargs=self.seasonal_test_kwargs,
                allowdrift=self.allowdrift,
                allowmean=self.allowmean,
                blambda=self.blambda,
                biasadj=self.biasadj,
                period=self.season_length,
            )

        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted AutoArima.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        fcst = forecast_arima(self.model_, h=h, xreg=X, level=level)
        mean = fcst["mean"]
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            res = {
                "mean": mean,
                **{f"lo-{l}": fcst["lower"][f"{l}%"] for l in reversed(level)},
                **{f"hi-{l}": fcst["upper"][f"{l}%"] for l in level},
            }
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted AutoArima insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = fitted_arima(self.model_)
        res = {"fitted": mean}
        if level is not None:
            se = np.sqrt(self.model_["sigma2"])
            res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient AutoARIMA predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenpus of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x) optional exogenous.
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        with np.errstate(invalid="ignore"):
            mod = auto_arima_f(
                x=y,
                d=self.d,
                D=self.D,
                max_p=self.max_p,
                max_q=self.max_q,
                max_P=self.max_P,
                max_Q=self.max_Q,
                max_order=self.max_order,
                max_d=self.max_d,
                max_D=self.max_D,
                start_p=self.start_p,
                start_q=self.start_q,
                start_P=self.start_P,
                start_Q=self.start_Q,
                stationary=self.stationary,
                seasonal=self.seasonal,
                ic=self.ic,
                stepwise=self.stepwise,
                nmodels=self.nmodels,
                trace=self.trace,
                approximation=self.approximation,
                method=self.method,
                truncate=self.truncate,
                xreg=X,
                test=self.test,
                test_kwargs=self.test_kwargs,
                seasonal_test=self.seasonal_test,
                seasonal_test_kwargs=self.seasonal_test_kwargs,
                allowdrift=self.allowdrift,
                allowmean=self.allowmean,
                blambda=self.blambda,
                biasadj=self.biasadj,
                period=self.season_length,
            )
        fcst = forecast_arima(mod, h, xreg=X_future, level=level)
        res = {"mean": fcst["mean"]}
        if fitted:
            res["fitted"] = fitted_arima(mod)
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                res = {
                    **res,
                    **{f"lo-{l}": fcst["lower"][f"{l}%"] for l in reversed(level)},
                    **{f"hi-{l}": fcst["upper"][f"{l}%"] for l in level},
                }
            if fitted:
                # add prediction intervals for fitted values
                se = np.sqrt(mod["sigma2"])
                res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forward(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Apply fitted ARIMA model to a new time series.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        if not hasattr(self, "model_"):
            raise Exception("You have to use the `fit` method first")
        with np.errstate(invalid="ignore"):
            mod = forward_arima(self.model_, y=y, xreg=X, method=self.method)
        fcst = forecast_arima(mod, h, xreg=X_future, level=level)
        res = {"mean": fcst["mean"]}
        if fitted:
            res["fitted"] = fitted_arima(mod)
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                res = {
                    **res,
                    **{f"lo-{l}": fcst["lower"][f"{l}%"] for l in reversed(level)},
                    **{f"hi-{l}": fcst["upper"][f"{l}%"] for l in level},
                }
            if fitted:
                # add prediction intervals for fitted values
                se = np.sqrt(mod["sigma2"])
                res = _add_fitted_pi(res=res, se=se, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 33
class AutoETS(_TS):
    """Automatic Exponential Smoothing model.

    Automatically selects the best ETS (Error, Trend, Seasonality)
    model using an information criterion. Default is Akaike Information Criterion (AICc), while particular models are estimated using maximum likelihood.
    The state-space equations can be determined based on their $M$ multiplicative, $A$ additive,
    $Z$ optimized or $N$ ommited components. The `model` string parameter defines the ETS equations:
    E in [$M, A, Z$], T in [$N, A, M, Z$], and S in [$N, A, M, Z$].

    For example when model='ANN' (additive error, no trend, and no seasonality), ETS will
    explore only a simple exponential smoothing.

    If the component is selected as 'Z', it operates as a placeholder to ask the AutoETS model
    to figure out the best parameter.

    Notes
    -----
    This implementation is a mirror of Hyndman's [forecast::ets](https://github.com/robjhyndman/forecast).

    References
    ----------
    [Rob J. Hyndman, Yeasmin Khandakar (2008). "Automatic Time Series Forecasting: The forecast package for R"](https://www.jstatsoft.org/article/view/v027i03).

    [Hyndman, Rob, et al (2008). "Forecasting with exponential smoothing: the state space approach"](https://robjhyndman.com/expsmooth/).

    Parameters
    ----------
    model : str
        Controlling state-space-equations.
    season_length : int
        Number of observations per unit of time. Ex: 24 Hourly data.
    damped : bool
        A parameter that 'dampens' the trend.
    phi : float, optional (default=None)
        Smoothing parameter for trend damping. Only used when `damped=True`.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals],
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        season_length: int = 1,
        model: str = "ZZZ",
        damped: Optional[bool] = None,
        phi: Optional[float] = None,
        alias: str = "AutoETS",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.season_length = season_length
        self.model = model
        self.damped = damped
        if phi is not None:
            if not isinstance(phi, float):
                raise ValueError("phi must be `None` or float.")
            if not _PHI_LOWER <= phi <= _PHI_UPPER:
                raise ValueError(f"Valid range for phi is [{_PHI_LOWER}, {_PHI_UPPER}]")
        self.phi = phi
        self.alias = alias
        self.prediction_intervals = prediction_intervals

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the Exponential Smoothing model.

        Fit an Exponential Smoothing model to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            Exponential Smoothing fitted model.
        """
        self.model_ = ets_f(
            y, m=self.season_length, model=self.model, damped=self.damped, phi=self.phi
        )
        self.model_["actual_residuals"] = y - self.model_["fitted"]
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self, h: int, X: Optional[np.ndarray] = None, level: Optional[List[int]] = None
    ):
        """Predict with fitted Exponential Smoothing.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenpus of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        fcst = forecast_ets(self.model_, h=h, level=level)
        res = {"mean": fcst["mean"]}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            res = {
                **res,
                **{f"lo-{l}": fcst[f"lo-{l}"] for l in reversed(level)},
                **{f"hi-{l}": fcst[f"hi-{l}"] for l in level},
            }
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted Exponential Smoothing insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            residuals = self.model_["actual_residuals"]
            se = _calculate_sigma(residuals, len(residuals) - self.model_["n_params"])
            res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient Exponential Smoothing predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenpus of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mod = ets_f(
            y, m=self.season_length, model=self.model, damped=self.damped, phi=self.phi
        )
        fcst = forecast_ets(mod, h=h, level=level)
        keys = ["mean"]
        if fitted:
            keys.append("fitted")
        res = {key: fcst[key] for key in keys}
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                res = {
                    **res,
                    **{f"lo-{l}": fcst[f"lo-{l}"] for l in reversed(level)},
                    **{f"hi-{l}": fcst[f"hi-{l}"] for l in level},
                }
            if fitted:
                # add prediction intervals for fitted values
                se = _calculate_sigma(y - mod["fitted"], len(y) - mod["n_params"])
                res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forward(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Apply fitted Exponential Smoothing model to a new time series.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenpus of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        if not hasattr(self, "model_"):
            raise Exception("You have to use the `fit` method first")
        mod = forward_ets(self.model_, y=y)
        fcst = forecast_ets(mod, h=h, level=level)
        keys = ["mean"]
        if fitted:
            keys.append("fitted")
        res = {key: fcst[key] for key in keys}
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                res = {
                    **res,
                    **{f"lo-{l}": fcst[f"lo-{l}"] for l in reversed(level)},
                    **{f"hi-{l}": fcst[f"hi-{l}"] for l in level},
                }
            if fitted:
                # add prediction intervals for fitted values
                se = _calculate_sigma(y - mod["fitted"], len(y) - mod["n_params"])
                res = _add_fitted_pi(res=res, se=se, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 48
class ETS(AutoETS):
    @classmethod
    def _warn(cls):
        warnings.warn(
            "`ETS` will be deprecated in future versions of `StatsForecast`. Please use `AutoETS` instead.",
            category=FutureWarning,
            stacklevel=2,
        )

    def __init__(
        self,
        season_length: int = 1,
        model: str = "ZZZ",
        damped: Optional[bool] = None,
        phi: Optional[float] = None,
        alias: str = "ETS",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        ETS._warn()
        super().__init__(
            season_length=season_length,
            model=model,
            damped=damped,
            phi=phi,
            alias=alias,
            prediction_intervals=prediction_intervals,
        )

    def __repr__(self):
        return self.alias

# %% ../nbs/src/core/models.ipynb 53
class AutoCES(_TS):
    """Complex Exponential Smoothing model.

    Automatically selects the best Complex Exponential Smoothing
    model using an information criterion. Default is Akaike Information Criterion (AICc), while particular
    models are estimated using maximum likelihood.
    The state-space equations can be determined based on their $S$ simple, $P$ parial,
    $Z$ optimized or $N$ ommited components. The `model` string parameter defines the
    kind of CES model: $N$ for simple CES (withous seasonality), $S$ for simple seasonality (lagged CES),
    $P$ for partial seasonality (without complex part), $F$ for full seasonality (lagged CES
    with real and complex seasonal parts).

    If the component is selected as 'Z', it operates as a placeholder to ask the AutoCES model
    to figure out the best parameter.

    References
    ----------
    [Svetunkov, Ivan & Kourentzes, Nikolaos. (2015). "Complex Exponential Smoothing". 10.13140/RG.2.1.3757.2562. ](https://onlinelibrary.wiley.com/doi/full/10.1002/nav.22074).

    Parameters
    ----------
    model : str
        Controlling state-space-equations.
    season_length : int
        Number of observations per unit of time. Ex: 24 Hourly data.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        season_length: int = 1,
        model: str = "Z",
        alias: str = "CES",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.season_length = season_length
        self.model = model
        self.alias = alias
        self.prediction_intervals = prediction_intervals

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the Complex Exponential Smoothing model.

        Fit the Complex Exponential Smoothing model to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            Complex Exponential Smoothing fitted model.
        """
        if is_constant(y):
            model = Naive(
                alias=self.alias, prediction_intervals=self.prediction_intervals
            )
            model.fit(y=y, X=X)
            return model
        self.model_ = auto_ces(y, m=self.season_length, model=self.model)
        self.model_["actual_residuals"] = y - self.model_["fitted"]
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self, h: int, X: Optional[np.ndarray] = None, level: Optional[List[int]] = None
    ):
        """Predict with fitted Exponential Smoothing.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level: List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        fcst = forecast_ces(self.model_, h=h, level=level)
        res = {"mean": fcst["mean"]}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            res = {
                **res,
                **{f"lo-{l}": fcst[f"lo-{l}"] for l in reversed(level)},
                **{f"hi-{l}": fcst[f"hi-{l}"] for l in level},
            }
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted Exponential Smoothing insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            residuals = self.model_["actual_residuals"]
            se = _calculate_sigma(residuals, self.model_["n"])
            res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient Complex Exponential Smoothing predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenpus of shape (h, n_x).
        level: List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        if is_constant(y):
            model = Naive(
                alias=self.alias, prediction_intervals=self.prediction_intervals
            )
            return model.forecast(
                y=y, h=h, X=X, X_future=X_future, level=level, fitted=fitted
            )
        mod = auto_ces(y, m=self.season_length, model=self.model)
        fcst = forecast_ces(mod, h, level=level)
        keys = ["mean"]
        if fitted:
            keys.append("fitted")
        res = {key: fcst[key] for key in keys}
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                res = {
                    **res,
                    **{f"lo-{l}": fcst[f"lo-{l}"] for l in reversed(level)},
                    **{f"hi-{l}": fcst[f"hi-{l}"] for l in level},
                }
            if fitted:
                # add prediction intervals for fitted values
                se = _calculate_sigma(y - mod["fitted"], len(y))
                res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forward(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Apply fitted Complex Exponential Smoothing to a new time series.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenpus of shape (h, n_x).
        level: List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        if not hasattr(self, "model_"):
            raise Exception("You have to use the `fit` method first")
        mod = forward_ces(self.model_, y=y)
        fcst = forecast_ces(mod, h, level=level)
        keys = ["mean"]
        if fitted:
            keys.append("fitted")
        res = {key: fcst[key] for key in keys}
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                res = {
                    **res,
                    **{f"lo-{l}": fcst[f"lo-{l}"] for l in reversed(level)},
                    **{f"hi-{l}": fcst[f"hi-{l}"] for l in level},
                }
            if fitted:
                # add prediction intervals for fitted values
                se = _calculate_sigma(y - mod["fitted"], len(y))
                res = _add_fitted_pi(res=res, se=se, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 71
class AutoTheta(_TS):
    """AutoTheta model.

    Automatically selects the best Theta (Standard Theta Model ('STM'),
    Optimized Theta Model ('OTM'), Dynamic Standard Theta Model ('DSTM'),
    Dynamic Optimized Theta Model ('DOTM')) model using mse.

    References
    ----------
    [Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016). "Models for optimising the theta method and their relationship to state space models". International Journal of Forecasting](https://www.sciencedirect.com/science/article/pii/S0169207016300243)

    Parameters
    ----------
    season_length : int
        Number of observations per unit of time. Ex: 24 Hourly data.
    decomposition_type : str
        Sesonal decomposition type, 'multiplicative' (default) or 'additive'.
    model : str
        Controlling Theta Model. By default searchs the best model.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        season_length: int = 1,
        decomposition_type: str = "multiplicative",
        model: Optional[str] = None,
        alias: str = "AutoTheta",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.season_length = season_length
        self.decomposition_type = decomposition_type
        self.model = model
        self.alias = alias
        self.prediction_intervals = prediction_intervals

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the AutoTheta model.

        Fit an AutoTheta model to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            AutoTheta fitted model.
        """
        self.model_ = auto_theta(
            y=y,
            m=self.season_length,
            model=self.model,
            decomposition_type=self.decomposition_type,
        )
        self.model_["fitted"] = y - self.model_["residuals"]
        self._store_cs(y, X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted AutoTheta.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        fcst = forecast_theta(self.model_, h=h, level=level)
        if self.prediction_intervals is not None and level is not None:
            fcst = self._add_predict_conformal_intervals(fcst, level)
        return fcst

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted AutoTheta insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            se = np.std(self.model_["residuals"][3:], ddof=1)
            res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient AutoTheta predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mod = auto_theta(
            y=y,
            m=self.season_length,
            model=self.model,
            decomposition_type=self.decomposition_type,
        )
        res = forecast_theta(mod, h, level=level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        if fitted:
            res["fitted"] = y - mod["residuals"]
        if level is not None and fitted:
            # add prediction intervals for fitted values
            se = np.std(mod["residuals"][3:], ddof=1)
            res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forward(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Apply fitted AutoTheta to a new time series.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        if not hasattr(self, "model_"):
            raise Exception("You have to use the `fit` method first")
        mod = forward_theta(self.model_, y=y)
        res = forecast_theta(mod, h, level=level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        if fitted:
            res["fitted"] = y - mod["residuals"]
        if level is not None and fitted:
            # add prediction intervals for fitted values
            se = np.std(mod["residuals"][3:], ddof=1)
            res = _add_fitted_pi(res=res, se=se, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 87
class ARIMA(_TS):
    """ARIMA model.

    AutoRegressive Integrated Moving Average model.

    References
    ----------
    [Rob J. Hyndman, Yeasmin Khandakar (2008). "Automatic Time Series Forecasting: The forecast package for R"](https://www.jstatsoft.org/article/view/v027i03).

    Parameters
    ----------
    order : tuple (default=(0, 0, 0))
        A specification of the non-seasonal part of the ARIMA model: the three components (p, d, q) are the AR order, the degree of differencing, and the MA order.
    season_length : int (default=1)
        Number of observations per unit of time. Ex: 24 Hourly data.
    seasonal_order : tuple (default=(0, 0, 0))
        A specification of the seasonal part of the ARIMA model.
        (P, D, Q) for the  AR order, the degree of differencing, the MA order.
    include_mean : bool (default=True)
        Should the ARIMA model include a mean term?
        The default is True for undifferenced series, False for differenced ones (where a mean would not affect the fit nor predictions).
    include_drift : bool (default=False)
        Should the ARIMA model include a linear drift term?
        (i.e., a linear regression with ARIMA errors is fitted.)
    include_constant : bool, optional (default=None)
        If True, then includ_mean is set to be True for undifferenced series and include_drift is set to be True for differenced series.
        Note that if there is more than one difference taken, no constant is included regardless of the value of this argument.
        This is deliberate as otherwise quadratic and higher order polynomial trends would be induced.
    blambda : float, optional (default=None)
        Box-Cox transformation parameter.
    biasadj : bool (default=False)
        Use adjusted back-transformed mean Box-Cox.
    method : str (default='CSS-ML')
        Fitting method: maximum likelihood or minimize conditional sum-of-squares.
        The default (unless there are missing values) is to use conditional-sum-of-squares to find starting values, then maximum likelihood.
    fixed : dict, optional (default=None)
        Dictionary containing fixed coefficients for the arima model. Example: `{'ar1': 0.5, 'ma2': 0.75}`.
        For autoregressive terms use the `ar{i}` keys. For its seasonal version use `sar{i}`.
        For moving average terms use the `ma{i}` keys. For its seasonal version use `sma{i}`.
        For intercept and drift use the `intercept` and `drift` keys.
        For exogenous variables use the `ex_{i}` keys.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    uses_exog = True

    def __init__(
        self,
        order: Tuple[int, int, int] = (0, 0, 0),
        season_length: int = 1,
        seasonal_order: Tuple[int, int, int] = (0, 0, 0),
        include_mean: bool = True,
        include_drift: bool = False,
        include_constant: Optional[bool] = None,
        blambda: Optional[float] = None,
        biasadj: bool = False,
        method: str = "CSS-ML",
        fixed: Optional[dict] = None,
        alias: str = "ARIMA",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.order = order
        self.season_length = season_length
        self.seasonal_order = seasonal_order
        self.include_mean = include_mean
        self.include_drift = include_drift
        self.include_constant = include_constant
        self.blambda = blambda
        self.biasadj = biasadj
        self.method = method
        self.fixed = fixed
        self.alias = alias
        self.prediction_intervals = prediction_intervals

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """
        Fit the model to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            Fitted model.
        """
        with np.errstate(invalid="ignore"):
            self.model_ = Arima(
                x=y,
                order=self.order,
                seasonal={"order": self.seasonal_order, "period": self.season_length},
                xreg=X,
                include_mean=self.include_mean,
                include_constant=self.include_constant,
                include_drift=self.include_drift,
                blambda=self.blambda,
                biasadj=self.biasadj,
                method=self.method,
                fixed=self.fixed,
            )
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted model.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        fcst = forecast_arima(self.model_, h=h, xreg=X, level=level)
        mean = fcst["mean"]
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            res = {
                "mean": mean,
                **{f"lo-{l}": fcst["lower"][f"{l}%"] for l in reversed(level)},
                **{f"hi-{l}": fcst["upper"][f"{l}%"] for l in level},
            }
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = fitted_arima(self.model_)
        res = {"fitted": mean}
        if level is not None:
            se = np.sqrt(self.model_["sigma2"])
            res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory efficient predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x) optional exogenous.
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        with np.errstate(invalid="ignore"):
            mod = Arima(
                x=y,
                order=self.order,
                seasonal={"order": self.seasonal_order, "period": self.season_length},
                xreg=X,
                include_mean=self.include_mean,
                include_constant=self.include_constant,
                include_drift=self.include_drift,
                blambda=self.blambda,
                biasadj=self.biasadj,
                method=self.method,
                fixed=self.fixed,
            )
        fcst = forecast_arima(mod, h, xreg=X_future, level=level)
        res = {"mean": fcst["mean"]}
        if fitted:
            res["fitted"] = fitted_arima(mod)
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                res = {
                    **res,
                    **{f"lo-{l}": fcst["lower"][f"{l}%"] for l in reversed(level)},
                    **{f"hi-{l}": fcst["upper"][f"{l}%"] for l in level},
                }
            if fitted:
                # add prediction intervals for fitted values
                se = np.sqrt(mod["sigma2"])
                res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forward(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Apply fitted model to a new time series.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        if not hasattr(self, "model_"):
            raise Exception("You have to use the `fit` method first")
        with np.errstate(invalid="ignore"):
            mod = forward_arima(self.model_, y=y, xreg=X, method=self.method)
        fcst = forecast_arima(mod, h, xreg=X_future, level=level)
        res = {"mean": fcst["mean"]}
        if fitted:
            res["fitted"] = fitted_arima(mod)
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                res = {
                    **res,
                    **{f"lo-{l}": fcst["lower"][f"{l}%"] for l in reversed(level)},
                    **{f"hi-{l}": fcst["upper"][f"{l}%"] for l in level},
                }
            if fitted:
                # add prediction intervals for fitted values
                se = np.sqrt(mod["sigma2"])
                res = _add_fitted_pi(res=res, se=se, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 102
class AutoRegressive(ARIMA):
    """Simple Autoregressive model.

    Parameters
    ----------
    lags : int or list
        Number of lags to include in the model.
        If an int is passed then all lags up to `lags` are considered.
        If a list, only the elements of the list are considered as lags.
    include_mean : bool (default=True)
        Should the AutoRegressive model include a mean term?
        The default is True for undifferenced series, False for differenced ones (where a mean would not affect the fit nor predictions).
    include_drift : bool (default=False)
        Should the AutoRegressive model include a linear drift term?
        (i.e., a linear regression with AutoRegressive errors is fitted.)
    blambda : float, optional (default=None)
        Box-Cox transformation parameter.
    biasadj : bool (default=False)
        Use adjusted back-transformed mean Box-Cox.
    method : str (default='CSS-ML')
        Fitting method: maximum likelihood or minimize conditional sum-of-squares.
        The default (unless there are missing values) is to use conditional-sum-of-squares to find starting values, then maximum likelihood.
    fixed : dict, optional (default=None)
        Dictionary containing fixed coefficients for the AutoRegressive model. Example: `{'ar1': 0.5, 'ar5': 0.75}`.
        For autoregressive terms use the `ar{i}` keys.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        lags: Tuple[int, List],
        include_mean: bool = True,
        include_drift: bool = False,
        blambda: Optional[float] = None,
        biasadj: bool = False,
        method: str = "CSS-ML",
        fixed: Optional[dict] = None,
        alias: str = "AutoRegressive",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        if isinstance(lags, int):
            order = (lags, 0, 0)
        elif isinstance(lags, list):
            order = (max(lags), 0, 0)
            fixed_lags = {
                f"ar{i+1}": np.nan if (i + 1) in lags else 0 for i in range(max(lags))
            }
            if fixed is not None:
                fixed_lags.update(fixed)
            fixed = fixed_lags
        else:
            raise ValueError(
                "Please provide an int or a list specifying the lags to use."
            )
        super().__init__(
            order=order,
            include_mean=include_mean,
            include_drift=include_drift,
            blambda=blambda,
            biasadj=biasadj,
            method=method,
            alias=alias,
            fixed=fixed,
            prediction_intervals=prediction_intervals,
        )

    def __repr__(self):
        return self.alias

# %% ../nbs/src/core/models.ipynb 117
@njit(nogil=NOGIL, cache=CACHE)
def _ses_fcst_mse(x: np.ndarray, alpha: float) -> Tuple[float, float, np.ndarray]:
    """Perform simple exponential smoothing on a series.

    This function returns the one step ahead prediction
    as well as the mean squared error of the fit.
    """
    smoothed = x[0]
    n = x.size
    mse = 0.0
    fitted = np.full(n, np.nan, np.float32)

    for i in range(1, n):
        smoothed = (alpha * x[i - 1] + (1 - alpha) * smoothed).item()
        error = x[i] - smoothed
        mse += error * error
        fitted[i] = smoothed

    mse /= n
    forecast = alpha * x[-1] + (1 - alpha) * smoothed
    return forecast, mse, fitted


def _ses_mse(alpha: float, x: np.ndarray) -> float:
    """Compute the mean squared error of a simple exponential smoothing fit."""
    _, mse, _ = _ses_fcst_mse(x, alpha)
    return mse


def _ses_forecast(x: np.ndarray, alpha: float) -> Tuple[float, np.ndarray]:
    """One step ahead forecast with simple exponential smoothing."""
    forecast, _, fitted = _ses_fcst_mse(x, alpha)
    return forecast, fitted


def _demand(x: np.ndarray) -> np.ndarray:
    """Extract the positive elements of a vector."""
    return x[x > 0]


def _intervals(x: np.ndarray) -> np.ndarray:
    """Compute the intervals between non zero elements of a vector."""
    nonzero_idxs = np.where(x != 0)[0]
    return np.diff(nonzero_idxs + 1, prepend=0)


def _probability(x: np.ndarray) -> np.ndarray:
    """Compute the element probabilities of being non zero."""
    return (x != 0).astype(np.int32)


def _optimized_ses_forecast(
    x: np.ndarray, bounds: Sequence[Tuple[float, float]] = [(0.1, 0.3)]
) -> Tuple[float, np.ndarray]:
    """Searches for the optimal alpha and computes SES one step forecast."""
    alpha = minimize(
        fun=_ses_mse, x0=(0,), args=(x,), bounds=bounds, method="L-BFGS-B"
    ).x[0]
    forecast, fitted = _ses_forecast(x, alpha)
    return forecast, fitted


def _chunk_sums(array: np.ndarray, chunk_size: int) -> np.ndarray:
    """Splits an array into chunks and returns the sum of each chunk.

    Incomplete chunks are discarded"""
    n_chunks = array.size // chunk_size
    n_elems = n_chunks * chunk_size
    return array[:n_elems].reshape(n_chunks, chunk_size).sum(axis=1)

# %% ../nbs/src/core/models.ipynb 118
def _ses(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
    alpha: float,  # smoothing parameter
) -> Dict[str, np.ndarray]:
    fcst, _, fitted_vals = _ses_fcst_mse(y, alpha)
    fcst = {"mean": _repeat_val(val=fcst, h=h)}
    if fitted:
        fcst["fitted"] = fitted_vals
    return fcst

# %% ../nbs/src/core/models.ipynb 119
class SimpleExponentialSmoothing(_TS):
    """SimpleExponentialSmoothing model.

    Uses a weighted average of all past observations where the weights decrease exponentially into the past.
    Suitable for data with no clear trend or seasonality.
    Assuming there are $t$ observations, the one-step forecast is given by: $\hat{y}_{t+1} = \\alpha y_t + (1-\\alpha) \hat{y}_{t-1}$

    The rate $0 \leq \\alpha \leq 1$ at which the weights decrease is called the smoothing parameter. When $\\alpha = 1$, SES is equal to the naive method.

    References
    ----------
    [Charles C Holt (1957). “Forecasting seasonals and trends by exponentially weighted moving averages”](https://doi.org/10.1016/j.ijforecast).

    Parameters
    ----------
    alpha : float
        Smoothing parameter.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        alpha: float,
        alias: str = "SES",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.alpha = alpha
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the SimpleExponentialSmoothing model.

        Fit an SimpleExponentialSmoothing to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            SimpleExponentialSmoothing fitted model.
        """
        mod = _ses(y=y, alpha=self.alpha, h=1, fitted=True)
        self.model_ = dict(mod)
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted SimpleExponentialSmoothing.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val(val=self.model_["mean"][0], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception("You must pass `prediction_intervals` to " "compute them.")
        return res

    def predict_in_sample(self):
        """Access fitted SimpleExponentialSmoothing insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions.

        """
        res = {"fitted": self.model_["fitted"]}
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient SimpleExponentialSmoothing predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _ses(y=y, h=h, fitted=fitted, alpha=self.alpha)
        res = dict(res)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        else:
            raise Exception("You must pass `prediction_intervals` to " "compute them.")
        return res

# %% ../nbs/src/core/models.ipynb 131
def _ses_optimized(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
):
    fcst_, fitted_vals = _optimized_ses_forecast(y, [(0.01, 0.99)])
    mean = _repeat_val(val=fcst_, h=h)
    fcst = {"mean": mean}
    if fitted:
        fcst["fitted"] = fitted_vals
    return fcst

# %% ../nbs/src/core/models.ipynb 132
class SimpleExponentialSmoothingOptimized(_TS):
    """SimpleExponentialSmoothing model.

    Uses a weighted average of all past observations where the weights decrease exponentially into the past.
    Suitable for data with no clear trend or seasonality.
    Assuming there are $t$ observations, the one-step forecast is given by: $\hat{y}_{t+1} = \\alpha y_t + (1-\\alpha) \hat{y}_{t-1}$

    The smoothing parameter $\\alpha^*$ is optimized by square error minimization.

    References
    ----------
    [Charles C Holt (1957). “Forecasting seasonals and trends by exponentially weighted moving averages”](https://doi.org/10.1016/j.ijforecast).

    Parameters
    ----------
    alias: str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        This is required for generating future prediction intervals.
    """

    def __init__(
        self,
        alias: str = "SESOpt",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the SimpleExponentialSmoothingOptimized model.

        Fit an SimpleExponentialSmoothingOptimized to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            SimpleExponentialSmoothingOptimized fitted model.
        """
        mod = _ses_optimized(y=y, h=1, fitted=True)
        self.model_ = dict(mod)
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted SimpleExponentialSmoothingOptimized.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val(val=self.model_["mean"][0], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

    def predict_in_sample(self):
        """Access fitted SimpleExponentialSmoothingOptimized insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient SimpleExponentialSmoothingOptimized predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _ses_optimized(y=y, h=h, fitted=fitted)
        res = dict(res)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

# %% ../nbs/src/core/models.ipynb 144
def _seasonal_exponential_smoothing(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
    season_length: int,  # length of season
    alpha: float,  # smoothing parameter
) -> Dict[str, np.ndarray]:
    n = y.size
    if n < season_length:
        return {"mean": np.full(h, np.nan, np.float32)}
    season_vals = np.empty(season_length, np.float32)
    fitted_vals = np.full(y.size, np.nan, np.float32)
    for i in range(season_length):
        init_idx = i + n % season_length
        season_vals[i], fitted_vals[init_idx::season_length] = _ses_forecast(
            y[init_idx::season_length], alpha
        )
    out = _repeat_val_seas(season_vals=season_vals, h=h)
    fcst = {"mean": out}
    if fitted:
        fcst["fitted"] = fitted_vals
    return fcst

# %% ../nbs/src/core/models.ipynb 145
class SeasonalExponentialSmoothing(_TS):
    """SeasonalExponentialSmoothing model.

    Uses a weighted average of all past observations where the weights decrease exponentially into the past.
    Suitable for data with no clear trend or seasonality.
    Assuming there are $t$ observations and season $s$, the one-step forecast is given by:
    $\hat{y}_{t+1,s} = \\alpha y_t + (1-\\alpha) \hat{y}_{t-1,s}$

    Notes
    -----
    This method is an extremely simplified of Holt-Winter's method where the trend and level are set to zero.
    And a single seasonal smoothing parameter $\\alpha$ is shared across seasons.

    References
    ----------
    [Charles. C. Holt (1957). "Forecasting seasonals and trends by exponentially weighted moving averages", ONR Research Memorandum, Carnegie Institute of Technology 52.](https://www.sciencedirect.com/science/article/abs/pii/S0169207003001134).
    [Peter R. Winters (1960). "Forecasting sales by exponentially weighted moving averages". Management Science](https://pubsonline.informs.org/doi/abs/10.1287/mnsc.6.3.324).

    Parameters
    ----------
    alpha : float
        Smoothing parameter.
    season_length : int
        Number of observations per unit of time. Ex: 24 Hourly data.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        This is required for generating future prediction intervals.
    """

    def __init__(
        self,
        season_length: int,
        alpha: float,
        alias: str = "SeasonalES",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.season_length = season_length
        self.alpha = alpha
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the SeasonalExponentialSmoothing model.

        Fit an SeasonalExponentialSmoothing to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            SeasonalExponentialSmoothing fitted model.
        """
        mod = _seasonal_exponential_smoothing(
            y=y,
            season_length=self.season_length,
            alpha=self.alpha,
            fitted=True,
            h=self.season_length,
        )
        self.model_ = dict(mod)
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted SeasonalExponentialSmoothing.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val_seas(self.model_["mean"], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

    def predict_in_sample(self):
        """Access fitted SeasonalExponentialSmoothing insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient SeasonalExponentialSmoothing predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _seasonal_exponential_smoothing(
            y=y, h=h, fitted=fitted, alpha=self.alpha, season_length=self.season_length
        )
        res = dict(res)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

# %% ../nbs/src/core/models.ipynb 160
def _seasonal_ses_optimized(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
    season_length: int,  # season length
):
    n = y.size
    if n < season_length:
        return {"mean": np.full(h, np.nan, np.float32)}
    season_vals = np.empty(season_length, np.float32)
    fitted_vals = np.full(y.size, np.nan, np.float32)
    for i in range(season_length):
        init_idx = i + n % season_length
        season_vals[i], fitted_vals[init_idx::season_length] = _optimized_ses_forecast(
            y[init_idx::season_length], [(0.01, 0.99)]
        )
    out = _repeat_val_seas(season_vals=season_vals, h=h)
    fcst = {"mean": out}
    if fitted:
        fcst["fitted"] = fitted_vals
    return fcst

# %% ../nbs/src/core/models.ipynb 161
class SeasonalExponentialSmoothingOptimized(_TS):

    def __init__(
        self,
        season_length: int,
        alias: str = "SeasESOpt",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """SeasonalExponentialSmoothingOptimized model.

        Uses a weighted average of all past observations where the weights decrease exponentially into the past.
        Suitable for data with no clear trend or seasonality.
        Assuming there are $t$ observations and season $s$, the one-step forecast is given by:
        $\hat{y}_{t+1,s} = \\alpha y_t + (1-\\alpha) \hat{y}_{t-1,s}$

        The smoothing parameter $\\alpha^*$ is optimized by square error minimization.

        Notes
        -----
        This method is an extremely simplified of Holt-Winter's method where the trend and level are set to zero.
        And a single seasonal smoothing parameter $\\alpha$ is shared across seasons.

        References
        ----------
        [Charles. C. Holt (1957). "Forecasting seasonals and trends by exponentially weighted moving averages", ONR Research Memorandum, Carnegie Institute of Technology 52.](https://www.sciencedirect.com/science/article/abs/pii/S0169207003001134).
        [Peter R. Winters (1960). "Forecasting sales by exponentially weighted moving averages". Management Science](https://pubsonline.informs.org/doi/abs/10.1287/mnsc.6.3.324).

        Parameters
        ----------
        season_length : int
            Number of observations per unit of time. Ex: 24 Hourly data.
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            This is required for generating future prediction intervals.
        """
        self.season_length = season_length
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the SeasonalExponentialSmoothingOptimized model.

        Fit an SeasonalExponentialSmoothingOptimized to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            SeasonalExponentialSmoothingOptimized fitted model.
        """
        mod = _seasonal_ses_optimized(
            y=y,
            season_length=self.season_length,
            fitted=True,
            h=self.season_length,
        )
        self.model_ = dict(mod)
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted SeasonalExponentialSmoothingOptimized.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val_seas(self.model_["mean"], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

    def predict_in_sample(self):
        """Access fitted SeasonalExponentialSmoothingOptimized insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient SeasonalExponentialSmoothingOptimized predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _seasonal_ses_optimized(
            y=y, h=h, fitted=fitted, season_length=self.season_length
        )
        res = dict(res)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

# %% ../nbs/src/core/models.ipynb 174
class Holt(AutoETS):
    """Holt's method.

    Also known as double exponential smoothing, Holt's method is an extension of exponential smoothing for series with a trend.
    This implementation returns the corresponding `ETS` model with additive (A) or multiplicative (M) errors (so either 'AAN' or 'MAN').

    References
    ----------
    [Rob J. Hyndman and George Athanasopoulos (2018). "Forecasting principles and practice, Methods with trend"](https://otexts.com/fpp3/holt.html).

    Parameters
    ----------
    season_length : int
        Number of observations per unit of time. Ex: 12 Monthly data.
    error_type : str
        The type of error of the ETS model. Can be additive (A) or multiplicative (M).
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        season_length: int = 1,
        error_type: str = "A",
        alias: str = "Holt",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):

        self.season_length = season_length
        self.error_type = error_type
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        model = error_type + "AN"
        super().__init__(
            season_length, model, alias=alias, prediction_intervals=prediction_intervals
        )

    def __repr__(self):
        return self.alias

# %% ../nbs/src/core/models.ipynb 188
class HoltWinters(AutoETS):
    """Holt-Winters' method.

    Also known as triple exponential smoothing, Holt-Winters' method is an extension of exponential smoothing for series that contain both trend and seasonality.
    This implementation returns the corresponding `ETS` model with additive (A) or multiplicative (M) errors (so either 'AAA' or 'MAM').

    References
    ----------
    [Rob J. Hyndman and George Athanasopoulos (2018). "Forecasting principles and practice, Methods with seasonality"](https://otexts.com/fpp3/holt-winters.html).

    Parameters
    ----------
    season_length : int
        Number of observations per unit of time. Ex: 12 Monthly data.
    error_type : str
        The type of error of the ETS model. Can be additive (A) or multiplicative (M).
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        season_length: int = 1,  # season length
        error_type: str = "A",  # error type
        alias: str = "HoltWinters",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.season_length = season_length
        self.error_type = error_type
        self.alias = alias
        model = error_type + "A" + error_type
        super().__init__(
            season_length, model, alias=alias, prediction_intervals=prediction_intervals
        )

    def __repr__(self):
        return self.alias

# %% ../nbs/src/core/models.ipynb 203
def _historic_average(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
) -> Dict[str, np.ndarray]:
    fcst = {"mean": _repeat_val(val=y.mean(), h=h)}
    if fitted:
        fitted_vals = _repeat_val(val=y.mean(), h=len(y))
        fcst["fitted"] = fitted_vals
    return fcst

# %% ../nbs/src/core/models.ipynb 204
class HistoricAverage(_TS):

    def __init__(
        self,
        alias: str = "HistoricAverage",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """HistoricAverage model.

        Also known as mean method. Uses a simple average of all past observations.
        Assuming there are $t$ observations, the one-step forecast is given by:
        $$\hat{y}_{t+1} = \\frac{1}{t} \sum_{j=1}^t y_j$$

        References
        ----------
        [Rob J. Hyndman and George Athanasopoulos (2018). "Forecasting principles and practice, Simple Methods"](https://otexts.com/fpp3/simple-methods.html).

        Parameters
        ----------
        alias: str
              Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            By default, the model will compute the native prediction
            intervals.
        """
        self.alias = alias
        self.prediction_intervals = prediction_intervals

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the HistoricAverage model.

        Fit an HistoricAverage to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
         Clean time series of shape (t, ).
        X : array-like
         Optional exogenous of shape (t, n_x).

        Returns
        -------
        self
            HistoricAverage fitted model.
        """
        mod = _historic_average(y, h=1, fitted=True)
        mod = dict(mod)
        residuals = y - mod["fitted"]
        mod["sigma"] = _calculate_sigma(residuals, len(residuals) - 1)
        mod["n"] = len(y)
        self.model_ = mod
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted HistoricAverage.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : numpy.array
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val(val=self.model_["mean"][0], h=h)
        res = {"mean": mean}

        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            sigma = self.model_["sigma"]
            sigmah = sigma * np.sqrt(1 + (1 / self.model_["n"]))
            pred_int = _calculate_intervals(res, level, h, sigmah)
            res = {**res, **pred_int}

        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted HistoricAverage insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            sigma = self.model_["sigma"]
            sigmah = sigma * np.sqrt(1 + (1 / self.model_["n"]))
            res = _add_fitted_pi(res, se=sigmah, level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient HistoricAverage predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : list[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        out = _historic_average(y=y, h=h, fitted=fitted or (level is not None))
        res = {"mean": out["mean"]}

        if fitted:
            res["fitted"] = out["fitted"]
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                residuals = y - out["fitted"]
                sigma = _calculate_sigma(residuals, len(residuals) - 1)
                sigmah = sigma * np.sqrt(1 + (1 / len(y)))
                pred_int = _calculate_intervals(out, level, h, sigmah)
                res = {**res, **pred_int}
            if fitted:
                residuals = y - out["fitted"]
                sigma = _calculate_sigma(residuals, len(residuals) - 1)
                sigmah = sigma * np.sqrt(1 + (1 / len(y)))
                res = _add_fitted_pi(res=res, se=sigmah, level=level)

        return res

# %% ../nbs/src/core/models.ipynb 217
class Naive(_TS):

    def __init__(
        self,
        alias: str = "Naive",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """Naive model.

        All forecasts have the value of the last observation:
        $\hat{y}_{t+1} = y_t$ for all $t$

        References
        ----------
        [Rob J. Hyndman and George Athanasopoulos (2018). "forecasting principles and practice, Simple Methods"](https://otexts.com/fpp3/simple-methods.html).

        Parameters
        ----------
        alias: str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            By default, the model will compute the native prediction
            intervals.
        """
        self.alias = alias
        self.prediction_intervals = prediction_intervals

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the Naive model.

        Fit an Naive to a time series (numpy.array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self:
            Naive fitted model.
        """
        mod = _naive(y, h=1, fitted=True)
        mod = dict(mod)
        residuals = y - mod["fitted"]
        sigma = _calculate_sigma(residuals, len(residuals) - 1)
        mod["sigma"] = sigma
        self.model_ = mod
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,  # forecasting horizon
        X: Optional[np.ndarray] = None,  # exogenous regressors
        level: Optional[List[int]] = None,  # confidence level
    ):
        """Predict with fitted Naive.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val(self.model_["mean"][0], h=h)
        res = {"mean": mean}

        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            steps = np.arange(1, h + 1)
            sigma = self.model_["sigma"]
            sigmah = sigma * np.sqrt(steps)
            pred_int = _calculate_intervals(res, level, h, sigmah)
            res = {**res, **pred_int}
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted Naive insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            res = _add_fitted_pi(res=res, se=self.model_["sigma"], level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient Naive predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n,).
        h: int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        out = _naive(y=y, h=h, fitted=fitted or (level is not None))
        res = {"mean": out["mean"]}
        if fitted:
            res["fitted"] = out["fitted"]
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                steps = np.arange(1, h + 1)
                residuals = y - out["fitted"]
                sigma = _calculate_sigma(residuals, len(residuals) - 1)
                sigmah = sigma * np.sqrt(steps)
                pred_int = _calculate_intervals(out, level, h, sigmah)
                res = {**res, **pred_int}
            if fitted:
                residuals = y - out["fitted"]
                sigma = _calculate_sigma(residuals, len(residuals) - 1)
                res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

    def forward(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Apply fitted model to an new/updated series.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n,).
        h: int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = self.forecast(
            y=y, h=h, X=X, X_future=X_future, level=level, fitted=fitted
        )
        return res

# %% ../nbs/src/core/models.ipynb 233
def _random_walk_with_drift(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
) -> Dict[str, np.ndarray]:
    slope = (y[-1] - y[0]) / (y.size - 1)
    mean = slope * (1 + np.arange(h, dtype=np.float32)) + y[-1]
    fcst = {
        "mean": mean.astype(np.float32, copy=False),
        "slope": np.array([slope], dtype=np.float32),
        "last_y": np.array([y[-1]], dtype=np.float32),
    }
    if fitted:
        fitted_vals = np.full(y.size, np.nan, dtype=np.float32)
        fitted_vals[1:] = (slope + y[:-1]).astype(np.float32)
        fcst["fitted"] = fitted_vals
    return fcst

# %% ../nbs/src/core/models.ipynb 234
class RandomWalkWithDrift(_TS):

    def __init__(
        self,
        alias: str = "RWD",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """RandomWalkWithDrift model.

        A variation of the naive method allows the forecasts to change over time.
        The amout of change, called drift, is the average change seen in the historical data.

        $$\hat{y}_{t+1} = y_t+\\frac{1}{t-1}\sum_{j=1}^t (y_j-y_{j-1}) = y_t+ \\frac{y_t-y_1}{t-1}$$

        From the previous equation, we can see that this is equivalent to extrapolating a line between
        the first and the last observation.

        References
        ----------
        [Rob J. Hyndman and George Athanasopoulos (2018). "forecasting principles and practice, Simple Methods"](https://otexts.com/fpp3/simple-methods.html).

        Parameters
        ----------
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            By default, the model will compute the native prediction
            intervals.
        """
        self.alias = alias
        self.prediction_intervals = prediction_intervals

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the RandomWalkWithDrift model.

        Fit an RandomWalkWithDrift to a time series (numpy array) `y`.

        Parameters
        ----------
        y: numpy.array
            Clean time series of shape (t, ).

        Returns
        -------
        self :
            RandomWalkWithDrift fitted model.
        """
        mod = _random_walk_with_drift(y, h=1, fitted=True)
        mod = dict(mod)
        residuals = y - mod["fitted"]
        sigma = _calculate_sigma(residuals, len(residuals) - 1)
        mod["sigma"] = sigma
        mod["n"] = len(y)
        self.model_ = mod
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self, h: int, X: Optional[np.ndarray] = None, level: Optional[List[int]] = None
    ):
        """Predict with fitted RandomWalkWithDrift.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        hrange = np.arange(h, dtype=np.float32)
        mean = self.model_["slope"] * (1 + hrange) + self.model_["last_y"]
        res = {"mean": mean}

        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            steps = np.arange(1, h + 1)
            sigma = self.model_["sigma"]
            sigmah = sigma * np.sqrt(steps * (1 + steps / (self.model_["n"] - 1)))
            pred_int = _calculate_intervals(res, level, h, sigmah)
            res = {**res, **pred_int}
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted RandomWalkWithDrift insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            res = _add_fitted_pi(res=res, se=self.model_["sigma"], level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient RandomWalkWithDrift predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n,).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts: dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        out = _random_walk_with_drift(y=y, h=h, fitted=fitted or (level is not None))
        res = {"mean": out["mean"]}

        if fitted:
            res["fitted"] = out["fitted"]

        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                steps = np.arange(1, h + 1)
                residuals = y - out["fitted"]
                sigma = _calculate_sigma(residuals, len(residuals) - 1)
                sigmah = sigma * np.sqrt(steps * (1 + steps / (len(y) - 1)))
                pred_int = _calculate_intervals(out, level, h, sigmah)
                res = {**res, **pred_int}
            if fitted:
                residuals = y - out["fitted"]
                sigma = _calculate_sigma(residuals, len(residuals) - 1)
                res = _add_fitted_pi(res=res, se=sigma, level=level)

        return res

# %% ../nbs/src/core/models.ipynb 249
class SeasonalNaive(_TS):

    def __init__(
        self,
        season_length: int,
        alias: str = "SeasonalNaive",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """Seasonal naive model.

        A method similar to the naive, but uses the last known observation of the same period (e.g. the same month of the previous year) in order to capture seasonal variations.

        References
        ----------
        [Rob J. Hyndman and George Athanasopoulos (2018). "forecasting principles and practice, Simple Methods"](https://otexts.com/fpp3/simple-methods.html#seasonal-na%C3%AFve-method).

        Parameters
        ----------
        season_length : int
            Number of observations per unit of time. Ex: 24 Hourly data.
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            By default, the model will compute the native prediction
            intervals.
        """
        self.season_length = season_length
        self.alias = alias
        self.prediction_intervals = prediction_intervals

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the SeasonalNaive model.

        Fit an SeasonalNaive to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X: array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            SeasonalNaive fitted model.
        """
        mod = _seasonal_naive(
            y=y,
            season_length=self.season_length,
            h=self.season_length,
            fitted=True,
        )
        mod = dict(mod)
        residuals = y - mod["fitted"]
        mod["sigma"] = _calculate_sigma(residuals, len(y) - self.season_length)
        self.model_ = mod
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted Naive.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X: array-like
            Optional exogenous of shape (h, n_x).
        level: List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val_seas(season_vals=self.model_["mean"], h=h)
        res = {"mean": mean}

        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            k = np.floor((h - 1) / self.season_length)
            sigma = self.model_["sigma"]
            sigmah = sigma * np.sqrt(k + 1)
            pred_int = _calculate_intervals(res, level, h, sigmah)
            res = {**res, **pred_int}
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted SeasonalNaive insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            level = sorted(level)
            res = _add_fitted_pi(res=res, se=self.model_["sigma"], level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient SeasonalNaive predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        out = _seasonal_naive(
            y=y,
            h=h,
            fitted=fitted or (level is not None),
            season_length=self.season_length,
        )
        res = {"mean": out["mean"]}
        if fitted:
            res["fitted"] = out["fitted"]
        if level is not None:
            level = sorted(level)
            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)
                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 264
def _window_average(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
    window_size: int,  # window size
) -> Dict[str, np.ndarray]:
    if fitted:
        raise NotImplementedError("return fitted")
    if y.size < window_size:
        return {"mean": np.full(h, np.nan, np.float32)}
    wavg = y[-window_size:].mean()
    mean = _repeat_val(val=wavg, h=h)
    return {"mean": mean}

# %% ../nbs/src/core/models.ipynb 265
class WindowAverage(_TS):

    def __init__(
        self,
        window_size: int,
        alias: str = "WindowAverage",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """WindowAverage model.

        Uses the average of the last $k$ observations, with $k$ the length of the window.
        Wider windows will capture global trends, while narrow windows will reveal local trends.
        The length of the window selected should take into account the importance of past
        observations and how fast the series changes.

        References
        ----------
        [Rob J. Hyndman and George Athanasopoulos (2018). "forecasting principles and practice, Simple Methods"](https://otexts.com/fpp3/simple-methods.html).

        Parameters
        ----------
        window_size : int
            Size of truncated series on which average is estimated.
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            This is required for generating future prediction intervals.
        """
        self.window_size = window_size
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the WindowAverage model.

        Fit an WindowAverage to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            WindowAverage fitted model.
        """
        mod = _window_average(y=y, h=1, window_size=self.window_size, fitted=False)
        self.model_ = dict(mod)
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted WindowAverage.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : numpy.array
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val(self.model_["mean"][0], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

    def predict_in_sample(self):
        """Access fitted WindowAverage insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        raise NotImplementedError

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient WindowAverage predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _window_average(y=y, h=h, fitted=fitted, window_size=self.window_size)
        res = dict(res)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

# %% ../nbs/src/core/models.ipynb 276
def _seasonal_window_average(
    y: np.ndarray,
    h: int,
    fitted: bool,
    season_length: int,
    window_size: int,
) -> Dict[str, np.ndarray]:
    if fitted:
        raise NotImplementedError("return fitted")
    min_samples = season_length * window_size
    if y.size < min_samples:
        return {"mean": np.full(h, np.nan, np.float32)}
    season_avgs = y[-min_samples:].reshape(window_size, season_length).mean(axis=0)
    out = _repeat_val_seas(season_vals=season_avgs, h=h)
    return {"mean": out}

# %% ../nbs/src/core/models.ipynb 277
class SeasonalWindowAverage(_TS):

    def __init__(
        self,
        season_length: int,
        window_size: int,
        alias: str = "SeasWA",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """SeasonalWindowAverage model.

        An average of the last $k$ observations of the same period, with $k$ the length of the window.

        References
        ----------
        [Rob J. Hyndman and George Athanasopoulos (2018). "forecasting principles and practice, Simple Methods"](https://otexts.com/fpp3/simple-methods.html).

        Parameters
        ----------
        window_size : int
            Size of truncated series on which average is estimated.
        seasonal_length : int
            Number of observations per cycle.
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            This is required for generating future prediction intervals.
        """
        self.season_length = season_length
        self.window_size = window_size
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the SeasonalWindowAverage model.

        Fit an SeasonalWindowAverage to a time series (numpy array) `y`
        and optionally exogenous variables (numpy array) `X`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenpus of shape (t, n_x).

        Returns
        -------
        self :
            SeasonalWindowAverage fitted model.
        """
        mod = _seasonal_window_average(
            y=y,
            h=self.season_length,
            fitted=False,
            season_length=self.season_length,
            window_size=self.window_size,
        )
        self.model_ = dict(mod)
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted SeasonalWindowAverage.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val_seas(season_vals=self.model_["mean"], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

    def predict_in_sample(self):
        """Access fitted SeasonalWindowAverage insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        raise NotImplementedError

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient SeasonalWindowAverage predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n,).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.


        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _seasonal_window_average(
            y=y,
            h=h,
            fitted=fitted,
            season_length=self.season_length,
            window_size=self.window_size,
        )
        res = dict(res)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

# %% ../nbs/src/core/models.ipynb 289
def _chunk_forecast(y, aggregation_level):
    lost_remainder_data = len(y) % aggregation_level
    y_cut = y[lost_remainder_data:]
    aggregation_sums = _chunk_sums(y_cut, aggregation_level)
    sums_forecast, _ = _optimized_ses_forecast(aggregation_sums)
    return sums_forecast


@njit(nogil=NOGIL, cache=CACHE)
def _expand_fitted_demand(fitted: np.ndarray, y: np.ndarray) -> np.ndarray:
    out = np.empty_like(y)
    out[0] = np.nan
    fitted_idx = 0
    for i in range(1, y.size):
        if y[i - 1] > 0:
            fitted_idx += 1
            out[i] = fitted[fitted_idx]
        elif fitted_idx > 0:
            # if this entry is zero, the model didn't change
            out[i] = out[i - 1]
        else:
            # if we haven't seen any demand, use naive
            out[i] = y[i - 1]
    return out


@njit(nogil=NOGIL, cache=CACHE)
def _expand_fitted_intervals(fitted: np.ndarray, y: np.ndarray) -> np.ndarray:
    out = np.empty_like(y)
    out[0] = np.nan
    fitted_idx = 0
    for i in range(1, y.size):
        if y[i - 1] != 0:
            fitted_idx += 1
            if fitted[fitted_idx] == 0:
                # to avoid division by zero
                out[i] = 1
            else:
                out[i] = fitted[fitted_idx]
        elif fitted_idx > 0:
            # if this entry is zero, the model didn't change
            out[i] = out[i - 1]
        else:
            # if we haven't seen any intervals, use 1 to avoid division by zero
            out[i] = 1
    return out


def _adida(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
):
    if (y == 0).all():
        res = {"mean": np.zeros(h, dtype=np.float32)}
        if fitted:
            res["fitted"] = np.zeros(y.size, dtype=np.float32)
            res["fitted"][0] = np.nan
        return res
    y = _ensure_float(y)
    y_intervals = _intervals(y)
    mean_interval = y_intervals.mean()
    aggregation_level = round(mean_interval)
    sums_forecast = _chunk_forecast(y, aggregation_level)
    forecast = sums_forecast / aggregation_level
    res = {"mean": _repeat_val(val=forecast, h=h)}
    if fitted:
        warnings.warn("Computing fitted values for ADIDA is very expensive")
        fitted_aggregation_levels = np.round(
            y_intervals.cumsum() / np.arange(1, y_intervals.size + 1)
        )
        fitted_aggregation_levels = _expand_fitted_intervals(
            np.append(np.nan, fitted_aggregation_levels), y
        )[1:].astype(np.int32)

        sums_fitted = np.empty(y.size - 1, dtype=y.dtype)
        for i, agg_lvl in enumerate(fitted_aggregation_levels):
            sums_fitted[i] = _chunk_forecast(y[: i + 1], agg_lvl)

        res["fitted"] = np.append(np.nan, sums_fitted / fitted_aggregation_levels)
    return res

# %% ../nbs/src/core/models.ipynb 290
class ADIDA(_TS):

    def __init__(
        self,
        alias: str = "ADIDA",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """ADIDA model.

        Aggregate-Dissagregate Intermittent Demand Approach: Uses temporal aggregation to reduce the
        number of zero observations. Once the data has been agregated, it uses the optimized SES to
        generate the forecasts at the new level. It then breaks down the forecast to the original
        level using equal weights.

        ADIDA specializes on sparse or intermittent series are series with very few non-zero observations.
        They are notoriously hard to forecast, and so, different methods have been developed
        especifically for them.

        References
        ----------
        [Nikolopoulos, K., Syntetos, A. A., Boylan, J. E., Petropoulos, F., & Assimakopoulos, V. (2011). An aggregate–disaggregate intermittent demand approach (ADIDA) to forecasting: an empirical proposition and analysis. Journal of the Operational Research Society, 62(3), 544-554.](https://researchportal.bath.ac.uk/en/publications/an-aggregate-disaggregate-intermittent-demand-approach-adida-to-f).

        Parameters
        ----------
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            By default, the model will compute the native prediction
            intervals.
        """
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the ADIDA model.

        Fit an ADIDA to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).

        Returns
        -------
        self :
            ADIDA fitted model.
        """
        self.model_ = _adida(y=y, h=1, fitted=False)
        self._y = y
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted ADIDA.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val(val=self.model_["mean"][0], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception(
                "You have to instantiate the class with `prediction_intervals`"
                "to calculate them"
            )
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted ADIDA insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        fitted = _adida(y=self._y, h=1, fitted=True)["fitted"]
        res = {"fitted": fitted}
        if level is not None:
            sigma = _calculate_sigma(self._y - fitted, self._y.size)
            res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient ADIDA predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n,).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _adida(y=y, h=h, fitted=fitted)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        else:
            raise Exception(
                "You have to instantiate the class with `prediction_intervals`"
                "to calculate them"
            )
        if fitted:
            sigma = _calculate_sigma(y - res["fitted"], y.size)
            res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 302
def _croston_classic(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
):
    y = _ensure_float(y)
    # demand
    yd = _demand(y)
    if not yd.size:  # no demand
        return _naive(y=y, h=h, fitted=fitted)
    ydp, ydf = _ses_forecast(yd, 0.1)

    # intervals
    yi = _intervals(y)
    yip, yif = _ses_forecast(yi, 0.1)

    if yip != 0.0:
        mean = ydp / yip
    else:
        mean = ydp
    out = {"mean": _repeat_val(val=mean, h=h)}
    if fitted:
        ydf = _expand_fitted_demand(np.append(ydf, ydp), y)
        yif = _expand_fitted_intervals(np.append(yif, yip), y)
        out["fitted"] = ydf / yif
    return out

# %% ../nbs/src/core/models.ipynb 303
class CrostonClassic(_TS):

    def __init__(
        self,
        alias: str = "CrostonClassic",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """CrostonClassic model.

        A method to forecast time series that exhibit intermittent demand.
        It decomposes the original time series into a non-zero demand size $z_t$ and
        inter-demand intervals $p_t$. Then the forecast is given by:
        $$\hat{y}_t = \\frac{\hat{z}_t}{\hat{p}_t}$$

        where $\hat{z}_t$ and $\hat{p}_t$ are forecasted using SES. The smoothing parameter
        of both components is set equal to 0.1

        References
        ----------
        [Croston, J. D. (1972). Forecasting and stock control for intermittent demands. Journal of the Operational Research Society, 23(3), 289-303.](https://link.springer.com/article/10.1057/jors.1972.50)

        Parameters
        ----------
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            By default, the model will compute the native prediction
            intervals.
        """
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the CrostonClassic model.

        Fit an CrostonClassic to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).

        Returns
        -------
        self :
            CrostonClassic fitted model.
        """
        self.model_ = _croston_classic(y=y, h=1, fitted=True)
        self.model_["sigma"] = _calculate_sigma(y - self.model_["fitted"], y.size)
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted CrostonClassic.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val(val=self.model_["mean"][0], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception(
                "You have to instantiate the class with `prediction_intervals` to calculate them"
            )
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted CrostonClassic insample predictions.

        Parameters
        ----------
        level: List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            res = _add_fitted_pi(res=res, se=self.model_["sigma"], level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient CrostonClassic predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _croston_classic(y=y, h=h, fitted=fitted)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=dict(res), y=y, X=X, level=level)
        else:
            raise Exception(
                "You have to instantiate the class with `prediction_intervals` to calculate them"
            )
        if fitted:
            sigma = _calculate_sigma(y - res["fitted"], y.size)
            res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 314
def _croston_optimized(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
):
    y = _ensure_float(y)
    # demand
    yd = _demand(y)
    if not yd.size:
        return _naive(y=y, h=h, fitted=fitted)
    ydp, _ = _optimized_ses_forecast(yd)

    # intervals
    yi = _intervals(y)
    yip, _ = _optimized_ses_forecast(yi)

    if yip != 0.0:
        mean = ydp / yip
    else:
        mean = ydp
    out = {"mean": _repeat_val(val=mean, h=h)}
    if fitted:
        warnings.warn("Computing fitted values for CrostonOptimized is very expensive")
        ydf = np.empty(yd.size + 1, dtype=y.dtype)
        ydf[0] = np.nan
        for i in range(yd.size):
            ydf[i + 1] = _optimized_ses_forecast(yd[: i + 1])[0]

        yif = np.empty(yi.size + 1, dtype=y.dtype)
        yif[0] = np.nan
        for i in range(yi.size):
            yiff = _optimized_ses_forecast(yi[: i + 1])[0]
            if yiff == 0:
                yiff = 1.0
            yif[i + 1] = yiff

        ydf = _expand_fitted_demand(ydf, y)
        yif = _expand_fitted_intervals(yif, y)
        out["fitted"] = ydf / yif
    return out

# %% ../nbs/src/core/models.ipynb 315
class CrostonOptimized(_TS):

    def __init__(
        self,
        alias: str = "CrostonOptimized",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """CrostonOptimized model.

        A method to forecast time series that exhibit intermittent demand.
        It decomposes the original time series into a non-zero demand size $z_t$ and
        inter-demand intervals $p_t$. Then the forecast is given by:
        $$\hat{y}_t = \\frac{\hat{z}_t}{\hat{p}_t}$$

        A variation of the classic Croston's method where the smooting paramater is optimally
        selected from the range $[0.1,0.3]$. Both the non-zero demand $z_t$ and the inter-demand
        intervals $p_t$ are smoothed separately, so their smoothing parameters can be different.

        References
        ----------
        [Croston, J. D. (1972). Forecasting and stock control for intermittent demands. Journal of the Operational Research Society, 23(3), 289-303.](https://link.springer.com/article/10.1057/jors.1972.50).

        Parameters
        ----------
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            This is required for generating future prediction intervals.
        """
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the CrostonOptimized model.

        Fit an CrostonOptimized to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).

        Returns
        -------
        self :
            CrostonOptimized fitted model.
        """
        self.model_ = _croston_optimized(y=y, h=1, fitted=False)
        self._y = y
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted CrostonOptimized.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val(val=self.model_["mean"][0], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted CrostonOptimized insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        fitted = _croston_optimized(y=self._y, h=1, fitted=True)["fitted"]
        res = {"fitted": fitted}
        if level is not None:
            sigma = _calculate_sigma(self._y - fitted, self._y.size)
            res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient CrostonOptimized predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _croston_optimized(y=y, h=h, fitted=fitted)
        res = dict(res)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        if fitted:
            sigma = _calculate_sigma(y - res["fitted"], y.size)
            res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 326
def _croston_sba(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
) -> Dict[str, np.ndarray]:
    out = _croston_classic(y=y, h=h, fitted=fitted)
    out["mean"] *= 0.95
    if fitted:
        out["fitted"] *= 0.95
    return out

# %% ../nbs/src/core/models.ipynb 327
class CrostonSBA(_TS):

    def __init__(
        self,
        alias: str = "CrostonSBA",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """CrostonSBA model.

        A method to forecast time series that exhibit intermittent demand.
        It decomposes the original time series into a non-zero demand size $z_t$ and
        inter-demand intervals $p_t$. Then the forecast is given by:
        $$\hat{y}_t = \\frac{\hat{z}_t}{\hat{p}_t}$$

        A variation of the classic Croston's method that uses a debiasing factor, so that the
        forecast is given by:
        $$\hat{y}_t = 0.95  \\frac{\hat{z}_t}{\hat{p}_t}$$

        References
        ----------
        [Croston, J. D. (1972). Forecasting and stock control for intermittent demands. Journal of the Operational Research Society, 23(3), 289-303.](https://link.springer.com/article/10.1057/jors.1972.50).

        Parameters
        ----------
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            By default, the model will compute the native prediction
            intervals.
        """
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the CrostonSBA model.

        Fit an CrostonSBA to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).

        Returns
        -------
        self :
            CrostonSBA fitted model.
        """
        self.model_ = _croston_sba(y=y, h=1, fitted=True)
        self.model_["sigma"] = _calculate_sigma(y - self.model_["fitted"], y.size)
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted CrostonSBA.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val(val=self.model_["mean"][0], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception(
                "You have to instantiate the class with `prediction_intervals` to calculate them"
            )
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted CrostonSBA insample predictions.

        Parameters
        ----------
        level: List[float]
            Confidence levels (0-100) prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            res = _add_fitted_pi(res=res, se=self.model_["sigma"], level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient CrostonSBA predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _croston_sba(y=y, h=h, fitted=fitted)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=dict(res), y=y, X=X, level=level)
        else:
            raise Exception(
                "You have to instantiate the class with `prediction_intervals`"
                "to calculate them"
            )
        if fitted:
            sigma = _calculate_sigma(y - res["fitted"], y.size)
            res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 338
def _imapa(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: bool,  # fitted values
):
    if (y == 0).all():
        res = {"mean": np.zeros(h, dtype=np.float32)}
        if fitted:
            res["fitted"] = np.zeros(y.size, dtype=np.float32)
            res["fitted"][0] = np.nan
        return res
    y = _ensure_float(y)
    y_intervals = _intervals(y)
    mean_interval = y_intervals.mean().item()
    max_aggregation_level = round(mean_interval)
    forecasts = np.empty(max_aggregation_level, np.float32)
    for aggregation_level in range(1, max_aggregation_level + 1):
        lost_remainder_data = len(y) % aggregation_level
        y_cut = y[lost_remainder_data:]
        aggregation_sums = _chunk_sums(y_cut, aggregation_level)
        forecast, _ = _optimized_ses_forecast(aggregation_sums)
        forecasts[aggregation_level - 1] = forecast / aggregation_level
    forecast = forecasts.mean()
    res = {"mean": _repeat_val(val=forecast, h=h)}
    if fitted:
        warnings.warn("Computing fitted values for IMAPA is very expensive.")
        fitted_vals = np.empty_like(y)
        fitted_vals[0] = np.nan
        for i in range(y.size - 1):
            fitted_vals[i + 1] = _imapa(y[: i + 1], h=1, fitted=False)["mean"].item()
        res["fitted"] = fitted_vals
    return res

# %% ../nbs/src/core/models.ipynb 339
class IMAPA(_TS):

    def __init__(
        self,
        alias: str = "IMAPA",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """IMAPA model.

        Intermittent Multiple Aggregation Prediction Algorithm: Similar to ADIDA, but instead of
        using a single aggregation level, it considers multiple in order to capture different
        dynamics of the data. Uses the optimized SES to generate the forecasts at the new levels
        and then combines them using a simple average.

        References
        ----------
        [Syntetos, A. A., & Boylan, J. E. (2021). Intermittent demand forecasting: Context, methods and applications. John Wiley & Sons.](https://www.ifors.org/intermittent-demand-forecasting-context-methods-and-applications/).

        Parameters
        ----------
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            By default, the model will compute the native prediction
            intervals.
        """
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the IMAPA model.

        Fit an IMAPA to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).

        Returns
        -------
        self :
            IMAPA fitted model.
        """
        self.model_ = _imapa(y=y, h=1, fitted=False)
        self._y = y
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted IMAPA.

        Parameters
        ----------
        h : int
            Forecast horizon.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        """
        mean = _repeat_val(val=self.model_["mean"][0], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception(
                "You have to instantiate the class with `prediction_intervals`"
                "to calculate them"
            )
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted IMAPA insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        fitted = _imapa(y=self._y, h=1, fitted=True)["fitted"]
        res = {"fitted": fitted}
        if level is not None:
            sigma = _calculate_sigma(self._y - fitted, self._y.size)
            res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient IMAPA predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _imapa(y=y, h=h, fitted=fitted)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        else:
            raise Exception(
                "You have to instantiate the class with `prediction_intervals`"
                "to calculate them"
            )
        if fitted:
            sigma = _calculate_sigma(y - res["fitted"], y.size)
            res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 350
def _tsb(
    y: np.ndarray,  # time series
    h: int,  # forecasting horizon
    fitted: int,  # fitted values
    alpha_d: float,
    alpha_p: float,
) -> Dict[str, np.ndarray]:
    if (y == 0).all():
        res = {"mean": np.zeros(h, dtype=np.float32)}
        if fitted:
            res["fitted"] = np.zeros(y.size, dtype=np.float32)
            res["fitted"][0] = np.nan
        return res
    y = _ensure_float(y)
    yd = _demand(y)
    yp = _probability(y)
    ypf, ypft = _ses_forecast(yp, alpha_p)
    ydf, ydft = _ses_forecast(yd, alpha_d)
    res = {"mean": _repeat_val(val=ypf * ydf, h=h)}
    if fitted:
        ydft = _expand_fitted_demand(np.append(ydft, ydf), y)
        res["fitted"] = ypft * ydft
    return res

# %% ../nbs/src/core/models.ipynb 351
class TSB(_TS):

    def __init__(
        self,
        alpha_d: float,
        alpha_p: float,
        alias: str = "TSB",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        """TSB model.

        Teunter-Syntetos-Babai: A modification of Croston's method that replaces the inter-demand
        intervals with the demand probability $d_t$, which is defined as follows.

        $$
        d_t = \\begin{cases}
            1  & \\text{if demand occurs at time t} \\\
            0  & \\text{otherwise.}
        \\end{cases}
        $$

        Hence, the forecast is given by

        $$\hat{y}_t= \hat{d}_t\hat{z_t}$$

        Both $d_t$ and $z_t$ are forecasted using SES. The smooting paramaters of each may differ,
        like in the optimized Croston's method.

        References
        ----------
        [Teunter, R. H., Syntetos, A. A., & Babai, M. Z. (2011). Intermittent demand: Linking forecasting to inventory obsolescence. European Journal of Operational Research, 214(3), 606-615.](https://www.sciencedirect.com/science/article/abs/pii/S0377221711004437)

        Parameters
        ----------
        alpha_d : float
            Smoothing parameter for demand.
        alpha_p : float
            Smoothing parameter for probability.
        alias : str
            Custom name of the model.
        prediction_intervals : Optional[ConformalIntervals]
            Information to compute conformal prediction intervals.
            This is required for generating future prediction intervals.
        """
        self.alpha_d = alpha_d
        self.alpha_p = alpha_p
        self.alias = alias
        self.prediction_intervals = prediction_intervals
        self.only_conformal_intervals = True

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the TSB model.

        Fit an TSB to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).

        Returns
        -------
        self :
            TSB fitted model.
        """
        self.model_ = _tsb(
            y=y, h=1, fitted=True, alpha_d=self.alpha_d, alpha_p=self.alpha_p
        )
        self.model_["sigma"] = _calculate_sigma(y - self.model_["fitted"], y.size)
        self._store_cs(y=y, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted TSB.

        Parameters
        ----------
        h : int
            Forecast horizon.
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = _repeat_val(self.model_["mean"][0], h=h)
        res = {"mean": mean}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted TSB insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            res = _add_fitted_pi(res=res, se=self.model_["sigma"], level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient TSB predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = _tsb(y=y, h=h, fitted=fitted, alpha_d=self.alpha_d, alpha_p=self.alpha_p)
        res = dict(res)
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        if fitted:
            sigma = _calculate_sigma(y - res["fitted"], y.size)
            res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 363
def _predict_mstl_seas(mstl_ob, h, season_length):
    seasoncolumns = mstl_ob.filter(regex="seasonal*").columns
    nseasons = len(seasoncolumns)
    seascomp = np.full((h, nseasons), np.nan)
    seasonal_periods = (
        [season_length] if isinstance(season_length, int) else season_length
    )
    for i in range(nseasons):
        mp = seasonal_periods[i]
        colname = seasoncolumns[i]
        seascomp[:, i] = np.tile(
            mstl_ob[colname].values[-mp:], trunc(1 + (h - 1) / mp)
        )[:h]
    lastseas = seascomp.sum(axis=1)
    return lastseas

# %% ../nbs/src/core/models.ipynb 364
class MSTL(_TS):
    """MSTL model.

    The MSTL (Multiple Seasonal-Trend decomposition using LOESS) decomposes the time series
    in multiple seasonalities using LOESS. Then forecasts the trend using
    a custom non-seaonal model and each seasonality using a SeasonalNaive model.

    References
    ----------
    [Bandara, Kasun & Hyndman, Rob & Bergmeir, Christoph. (2021). "MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns".](https://arxiv.org/abs/2107.13462).

    Parameters
    ----------
    season_length : Union[int, List[int]
        Number of observations per unit of time. For multiple seasonalities use a list.
    trend_forecaster : model, default=AutoETS(model='ZZN')
        StatsForecast model used to forecast the trend component.
    stl_kwargs : dict
        Extra arguments to pass to [`statsmodels.tsa.seasonal.STL`](https://www.statsmodels.org/dev/generated/statsmodels.tsa.seasonal.STL.html#statsmodels.tsa.seasonal.STL).
        The `period` and `seasonal` arguments are reserved.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        season_length: Union[int, List[int]],
        trend_forecaster=AutoETS(model="ZZN"),
        stl_kwargs: Optional[Dict] = None,
        alias: str = "MSTL",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):

        # check ETS model doesnt have seasonality
        if repr(trend_forecaster) == "AutoETS":
            if trend_forecaster.model[2] != "N":
                raise Exception(
                    "Trend forecaster should not adjust " "seasonal models."
                )
        # check if trend forecaster has season_length=1
        if hasattr(trend_forecaster, "season_length"):
            if trend_forecaster.season_length != 1:
                raise Exception(
                    "Trend forecaster should not adjust "
                    "seasonal models. Please pass `season_length=1` "
                    "to your trend forecaster"
                )
        if isinstance(season_length, int):
            season_length = [season_length]
        else:
            season_length = sorted(season_length)
        self.season_length = season_length
        self.trend_forecaster = trend_forecaster
        self.prediction_intervals = prediction_intervals
        self.alias = alias

        if self.trend_forecaster.prediction_intervals is None and (
            self.prediction_intervals is not None
        ):
            self.trend_forecaster.prediction_intervals = prediction_intervals
        self.stl_kwargs = dict() if stl_kwargs is None else stl_kwargs

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the MSTL model.

        Fit MSTL to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X: array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self :
            MSTL fitted model.
        """
        self.model_ = mstl(
            x=y,
            period=self.season_length,
            stl_kwargs=self.stl_kwargs,
        )
        x_sa = self.model_[["trend", "remainder"]].sum(axis=1).values
        self.trend_forecaster = self.trend_forecaster.new().fit(y=x_sa, X=X)
        self._store_cs(y=x_sa, X=X)
        return self

    def predict(
        self,
        h: int,
        X: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
    ):
        """Predict with fitted MSTL.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        kwargs: Dict[str, Any] = {"h": h, "X": X}
        if self.trend_forecaster.prediction_intervals is None:
            kwargs["level"] = level
        res = self.trend_forecaster.predict(**kwargs)
        seas = _predict_mstl_seas(self.model_, h=h, season_length=self.season_length)
        res = {key: val + seas for key, val in res.items()}
        if level is None or self.trend_forecaster.prediction_intervals is None:
            return res
        level = sorted(level)
        if self.trend_forecaster.prediction_intervals is not None:
            res = self.trend_forecaster._add_predict_conformal_intervals(res, level)
        else:
            raise Exception(
                "You have to instantiate either the trend forecaster class or MSTL class with `prediction_intervals` to calculate them"
            )
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted MSTL insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        res = self.trend_forecaster.predict_in_sample(level=level)
        seas = self.model_.filter(regex="seasonal*").sum(axis=1).values
        res = {key: val + seas for key, val in res.items()}
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient MSTL predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        model_ = mstl(
            x=y,
            period=self.season_length,
            stl_kwargs=self.stl_kwargs,
        )
        x_sa = model_[["trend", "remainder"]].sum(axis=1).values
        kwargs = {"y": x_sa, "h": h, "X": X, "X_future": X_future, "fitted": fitted}
        if fitted or self.trend_forecaster.prediction_intervals is None:
            kwargs["level"] = level
        res = self.trend_forecaster.forecast(**kwargs)
        if level is not None:
            level = sorted(level)
            if self.trend_forecaster.prediction_intervals is not None:
                res = self.trend_forecaster._add_conformal_intervals(
                    fcst=res, y=x_sa, X=X, level=level
                )
            elif f"lo-{level[0]}" not in res:
                raise Exception(
                    "You have to instantiate either the trend forecaster class or MSTL class with `prediction_intervals` to calculate them"
                )
        # reseasonalize results
        seas_h = _predict_mstl_seas(model_, h=h, season_length=self.season_length)
        seas_insample = model_.filter(regex="seasonal*").sum(axis=1).values
        res = {
            key: val + (seas_insample if "fitted" in key else seas_h)
            for key, val in res.items()
        }
        return res

    def forward(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Apply fitted MSTL model to a new time series.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        if not hasattr(self.trend_forecaster, "model_"):
            raise Exception("You have to use the `fit` method first")
        model_ = mstl(
            x=y,
            period=self.season_length,
            stl_kwargs=self.stl_kwargs,
        )
        x_sa = model_[["trend", "remainder"]].sum(axis=1).values
        kwargs = {"y": x_sa, "h": h, "X": X, "X_future": X_future, "fitted": fitted}
        if fitted or self.trend_forecaster.prediction_intervals is None:
            kwargs["level"] = level
        res = self.trend_forecaster.forward(**kwargs)
        if level is not None:
            level = sorted(level)
            if self.trend_forecaster.prediction_intervals is not None:
                res = self.trend_forecaster._add_conformal_intervals(
                    fcst=res, y=x_sa, X=X, level=level
                )
        # reseasonalize results
        seas_h = _predict_mstl_seas(model_, h=h, season_length=self.season_length)
        seas_insample = model_.filter(regex="seasonal*").sum(axis=1).values
        res = {
            key: val + (seas_insample if "fitted" in key else seas_h)
            for key, val in res.items()
        }
        return res

# %% ../nbs/src/core/models.ipynb 380
class TBATS(_TS):
    """Trigonometric Box-Cox transform, ARMA errors, Trend and Seasonal components (TBATS) model.

    TBATS is an innovations state space model framework used for forecasting time series with multiple seasonalities. It uses a Box-Cox tranformation, ARMA errors, and a trigonometric representation of the seasonal patterns based on Fourier series.

    The name TBATS is an acronym for the key features of the model: Trigonometric, Box-Cox transform, ARMA errors, Trend, and Seasonal components.

    References
    ----------
    - [De Livera, A. M., Hyndman, R. J., & Snyder, R. D. (2011). Forecasting time series with complex seasonal patterns using exponential smoothing. Journal of the American statistical association, 106(496), 1513-1527.](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=f3de25596ab60ef0e886366826bf58a02b35a44f)

    - [De Livera, Alysha M (2017). Modeling time series with complex seasonal patterns using exponential smoothing. Monash University. Thesis.](https://doi.org/10.4225/03/589299681de3d)

    Parameters
    ----------
    season_length : int or list of int.
        Number of observations per unit of time. Ex: 24 Hourly data.
    use_boxcox : bool (default=True)
        Whether or not to use a Box-Cox transformation.
    bc_lower_bound : float (default=0.0)
        Lower bound for the Box-Cox transformation.
    bc_upper_bound : float (default=1.0)
        Upper bound for the Box-Cox transformation.
    use_trend : bool (default=True)
        Whether or not to use a trend component.
    use_damped_trend : bool (default=False)
        Whether or not to dampen the trend component.
    use_arma_errors : bool (default=False)
        Whether or not to use a ARMA errors.
    alias : str
        Custom name of the model.
    """

    @_old_kw_to_pos(["seasonal_periods"], [1])
    def __init__(
        self,
        season_length: Union[int, List[int]],
        use_boxcox: Optional[bool] = True,
        bc_lower_bound: float = 0.0,
        bc_upper_bound: float = 1.0,
        use_trend: Optional[bool] = True,
        use_damped_trend: Optional[bool] = False,
        use_arma_errors: bool = False,
        alias: str = "TBATS",
        *,
        seasonal_periods=None,  # noqa: ARG002
    ):
        if isinstance(season_length, int):
            season_length = [season_length]
        self.season_length = list(season_length)
        self.use_boxcox = use_boxcox
        self.bc_lower_bound = bc_lower_bound
        self.bc_upper_bound = bc_upper_bound
        self.use_trend = use_trend
        self.use_damped_trend = use_damped_trend
        self.use_arma_errors = use_arma_errors
        self.alias = alias

    def __repr__(self):
        return self.alias

    def fit(self, y: np.ndarray, X: Optional[np.ndarray] = None):
        """Fit TBATS model.

        Fit TBATS model to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : numpy.array, optional (default=None)
            Ignored

        Returns
        -------
        self :
            TBATS model.
        """
        self.model_ = tbats_selection(
            y=y,
            seasonal_periods=self.season_length,
            use_boxcox=self.use_boxcox,
            bc_lower_bound=self.bc_lower_bound,
            bc_upper_bound=self.bc_upper_bound,
            use_trend=self.use_trend,
            use_damped_trend=self.use_damped_trend,
            use_arma_errors=self.use_arma_errors,
        )
        return self

    def predict(
        self, h: int, X: Optional[np.ndarray] = None, level: Optional[List[int]] = None
    ):
        """Predict with fitted TBATS model.

        Parameters
        ----------
        h : int
            Forecast horizon.
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        fcst = tbats_forecast(self.model_, h)
        res = {"mean": fcst["mean"]}
        if level is not None:
            level = sorted(level)
            sigmah = _compute_sigmah(self.model_, h)
            pred_int = _calculate_intervals(res, level, h, sigmah)
            res = {**res, **pred_int}
        if self.model_["BoxCox_lambda"] is not None:
            res_trans = {
                k: inv_boxcox(v, self.model_["BoxCox_lambda"]) for k, v in res.items()
            }
            for k, v in res_trans.items():
                res_trans[k] = np.where(np.isnan(v), res[k], v)
        else:
            res_trans = res
        return res_trans

    def predict_in_sample(self, level: Optional[Tuple[int]] = None):
        """Access fitted TBATS model predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"].ravel()}
        if level is not None:
            se = _calculate_sigma(self.model_["errors"], self.model_["errors"].shape[1])
            fitted_pred_int = _add_fitted_pi(res, se, level)
            res = {**res, **fitted_pred_int}
        if self.model_["BoxCox_lambda"] is not None:
            res_trans = {
                k: inv_boxcox(v, self.model_["BoxCox_lambda"]) for k, v in res.items()
            }
            for k, v in res_trans.items():
                res_trans[k] = np.where(np.isnan(v), res[k], v)
        else:
            res_trans = res
        return res_trans

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient TBATS model.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mod = tbats_selection(
            y=y,
            seasonal_periods=self.season_length,
            use_boxcox=self.use_boxcox,
            bc_lower_bound=self.bc_lower_bound,
            bc_upper_bound=self.bc_upper_bound,
            use_trend=self.use_trend,
            use_damped_trend=self.use_damped_trend,
            use_arma_errors=self.use_arma_errors,
        )
        fcst = tbats_forecast(mod, h)
        res = {"mean": fcst["mean"]}
        if fitted:
            res["fitted"] = mod["fitted"].ravel()
        if level is not None:
            level = sorted(level)
            sigmah = _compute_sigmah(mod, h)
            pred_int = _calculate_intervals(res, level, h, sigmah)
            res = {**res, **pred_int}
            if fitted:
                se = _calculate_sigma(mod["errors"], mod["errors"].shape[1])
                fitted_pred_int = _add_fitted_pi(res, se, level)
                res = {**res, **fitted_pred_int}
        if mod["BoxCox_lambda"] is not None:
            res_trans = {k: inv_boxcox(v, mod["BoxCox_lambda"]) for k, v in res.items()}
            for k, v in res_trans.items():
                res_trans[k] = np.where(np.isnan(v), res[k], v)
        else:
            res_trans = res
        return res_trans

# %% ../nbs/src/core/models.ipynb 388
class AutoTBATS(TBATS):
    """AutoTBATS model.

    Automatically selects the best TBATS model from all feasible combinations of the parameters use_boxcox, use_trend, use_damped_trend, and use_arma_errors.
    Selection is made using the AIC.
    Default value for use_arma_errors is True since this enables the evaluation of models with and without ARMA errors.

    References
    ----------
    - [De Livera, A. M., Hyndman, R. J., & Snyder, R. D. (2011). Forecasting time series with complex seasonal patterns using exponential smoothing. Journal of the American statistical association, 106(496), 1513-1527.](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=f3de25596ab60ef0e886366826bf58a02b35a44f)

    - [De Livera, Alysha M (2017). Modeling time series with complex seasonal patterns using exponential smoothing. Monash University. Thesis.](https://doi.org/10.4225/03/589299681de3d)

    Parameters
    ----------
    seasonal_periods : int or list of int.
        Number of observations per unit of time. Ex: 24 Hourly data.
    use_boxcox : bool (default=None)
        Whether or not to use a Box-Cox transformation. By default tries both.
    bc_lower_bound : float (default=0.0)
        Lower bound for the Box-Cox transformation.
    bc_upper_bound : float (default=1.0)
        Upper bound for the Box-Cox transformation.
    use_trend : bool (default=None)
        Whether or not to use a trend component. By default tries both.
    use_damped_trend : bool (default=None)
        Whether or not to dampen the trend component. By default tries both.
    use_arma_errors : bool (default=True)
        Whether or not to use a ARMA errors. Default is True and this evaluates both models.
    alias : str
        Custom name of the model.
    """

    @_old_kw_to_pos(["seasonal_periods"], [1])
    def __init__(
        self,
        season_length: Union[int, List[int]],
        use_boxcox: Optional[bool] = None,
        bc_lower_bound: float = 0.0,
        bc_upper_bound: float = 1.0,
        use_trend: Optional[bool] = None,
        use_damped_trend: Optional[bool] = None,
        use_arma_errors: bool = True,
        alias: str = "AutoTBATS",
        *,
        seasonal_periods=None  # noqa: ARG002
    ):
        super().__init__(
            season_length=season_length,
            use_boxcox=use_boxcox,
            bc_lower_bound=bc_lower_bound,
            bc_upper_bound=bc_upper_bound,
            use_trend=use_trend,
            use_damped_trend=use_damped_trend,
            use_arma_errors=use_arma_errors,
            alias=alias,
        )

# %% ../nbs/src/core/models.ipynb 398
class Theta(AutoTheta):
    """Standard Theta Method.

    References
    ----------
    [Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016). "Models for optimising the theta method and their relationship to state space models". International Journal of Forecasting](https://www.sciencedirect.com/science/article/pii/S0169207016300243)

    Parameters
    ----------
    season_length : int
        Number of observations per unit of time. Ex: 24 Hourly data.
    decomposition_type : str
        Sesonal decomposition type, 'multiplicative' (default) or 'additive'.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        season_length: int = 1,
        decomposition_type: str = "multiplicative",
        alias: str = "Theta",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        super().__init__(
            season_length=season_length,
            model="STM",
            decomposition_type=decomposition_type,
            alias=alias,
            prediction_intervals=prediction_intervals,
        )

# %% ../nbs/src/core/models.ipynb 412
class OptimizedTheta(AutoTheta):
    """Optimized Theta Method.

    References
    ----------
    [Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016). "Models for optimising the theta method and their relationship to state space models". International Journal of Forecasting](https://www.sciencedirect.com/science/article/pii/S0169207016300243)

    Parameters
    ----------
    season_length : int
        Number of observations per unit of time. Ex: 24 Hourly data.
    decomposition_type : str
        Sesonal decomposition type, 'multiplicative' (default) or 'additive'.
    alias : str
        Custom name of the model. Default `OptimizedTheta`.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        season_length: int = 1,
        decomposition_type: str = "multiplicative",
        alias: str = "OptimizedTheta",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        super().__init__(
            season_length=season_length,
            model="OTM",
            decomposition_type=decomposition_type,
            alias=alias,
            prediction_intervals=prediction_intervals,
        )

# %% ../nbs/src/core/models.ipynb 426
class DynamicTheta(AutoTheta):
    """Dynamic Standard Theta Method.

    References
    ----------
    [Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016). "Models for optimising the theta method and their relationship to state space models". International Journal of Forecasting](https://www.sciencedirect.com/science/article/pii/S0169207016300243)

    Parameters
    ----------
    season_length : int
        Number of observations per unit of time. Ex: 24 Hourly data.
    decomposition_type : str
        Sesonal decomposition type, 'multiplicative' (default) or 'additive'.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        season_length: int = 1,
        decomposition_type: str = "multiplicative",
        alias: str = "DynamicTheta",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        super().__init__(
            season_length=season_length,
            model="DSTM",
            decomposition_type=decomposition_type,
            alias=alias,
            prediction_intervals=prediction_intervals,
        )

# %% ../nbs/src/core/models.ipynb 440
class DynamicOptimizedTheta(AutoTheta):
    """Dynamic Optimized Theta Method.

    References
    ----------
    [Jose A. Fiorucci, Tiago R. Pellegrini, Francisco Louzada, Fotios Petropoulos, Anne B. Koehler (2016). "Models for optimising the theta method and their relationship to state space models". International Journal of Forecasting](https://www.sciencedirect.com/science/article/pii/S0169207016300243)

    Parameters
    ----------
    season_length : int
        Number of observations per unit of time. Ex: 24 Hourly data.
    decomposition_type : str
        Sesonal decomposition type, 'multiplicative' (default) or 'additive'.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        season_length: int = 1,
        decomposition_type: str = "multiplicative",
        alias: str = "DynamicOptimizedTheta",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        super().__init__(
            season_length=season_length,
            model="DOTM",
            decomposition_type=decomposition_type,
            alias=alias,
            prediction_intervals=prediction_intervals,
        )

# %% ../nbs/src/core/models.ipynb 455
class GARCH(_TS):
    """Generalized Autoregressive Conditional Heteroskedasticity (GARCH) model.

    A method for modeling time series that exhibit non-constant volatility over time.
    The GARCH model assumes that at time $t$, $y_t$ is given by:

    $$y_t = v_t \sigma_t$$

    with

    $$\sigma_t^2 = w + \sum_{i=1}^p a_i y_{t-i}^2 + \sum_{j=1}^q b_j \sigma_{t-j}^2$$.

    Here $v_t$ is a sequence of iid random variables with zero mean and unit variance.
    The coefficients $w$, $a_i$, $i=1,...,p$, and $b_j$, $j=1,...,q$ must satisfy the following conditions:

    1. $w > 0$ and $a_i, b_j \geq 0$ for all $i$ and $j$.
    2. $\sum_{k=1}^{max(p,q)} a_k + b_k < 1$. Here it is assumed that $a_i=0$ for $i>p$ and $b_j=0$ for $j>q$.

    The ARCH model is a particular case of the GARCH model when $q=0$.

    References
    ----------
    [Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica: Journal of the econometric society, 987-1007.](http://www.econ.uiuc.edu/~econ508/Papers/engle82.pdf)

    [Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of econometrics, 31(3), 307-327.](https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=7da8bfa5295375c1141d797e80065a599153c19d)

    [James D. Hamilton. Time Series Analysis Princeton University Press, Princeton, New Jersey, 1st Edition, 1994.](https://press.princeton.edu/books/hardcover/9780691042893/time-series-analysis)

    Parameters
    ----------
    p : int
        Number of lagged versions of the series.
    q: int
        Number of lagged versions of the volatility.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        p: int = 1,
        q: int = 1,
        alias: str = "GARCH",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.p = p
        self.q = q
        if q != 0:
            self.alias = alias + "(" + str(p) + "," + str(q) + ")"
        else:
            self.alias = alias + "(" + str(p) + ")"
        self.prediction_intervals = prediction_intervals

    def __repr__(self):
        return self.alias

    def fit(self, y: np.ndarray, X: Optional[np.ndarray] = None):
        """Fit GARCH model.

        Fit GARCH model to a time series (numpy array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).

        Returns
        -------
        self :
            GARCH model.
        """
        self.model_ = garch_model(y, p=self.p, q=self.q)
        self.model_["actual_residuals"] = y - self.model_["fitted"]
        self._store_cs(y, X)
        return self

    def predict(
        self, h: int, X: Optional[np.ndarray] = None, level: Optional[List[int]] = None
    ):
        """Predict with fitted GARCH model.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        fcst = garch_forecast(self.model_, h)
        res = {"mean": fcst["mean"], "sigma2": fcst["sigma2"]}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            quantiles = _quantiles(level)
            lo = res["mean"].reshape(-1, 1) - quantiles * res["sigma2"].reshape(-1, 1)
            hi = res["mean"].reshape(-1, 1) + quantiles * res["sigma2"].reshape(-1, 1)
            lo = lo[:, ::-1]
            lo = {f"lo-{l}": lo[:, i] for i, l in enumerate(reversed(level))}
            hi = {f"hi-{l}": hi[:, i] for i, l in enumerate(level)}
            res = {**res, **lo, **hi}
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted GARCH model predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            residuals = self.model_["actual_residuals"]
            se = _calculate_sigma(residuals, len(residuals) - 1)
            res = _add_fitted_pi(res=res, se=se, level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient GARCH model.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mod = garch_model(y, p=self.p, q=self.q)
        fcst = garch_forecast(mod, h)
        keys = ["mean", "sigma2"]
        if fitted:
            keys.append("fitted")
        res = {key: fcst[key] for key in keys}
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_predict_conformal_intervals(res, level)
            else:
                quantiles = _quantiles(level)
                lo = res["mean"].reshape(-1, 1) - quantiles * res["sigma2"].reshape(
                    -1, 1
                )
                hi = res["mean"].reshape(-1, 1) + quantiles * res["sigma2"].reshape(
                    -1, 1
                )
                lo = lo[:, ::-1]
                lo = {f"lo-{l}": lo[:, i] for i, l in enumerate(reversed(level))}
                hi = {f"hi-{l}": hi[:, i] for i, l in enumerate(level)}
                res = {**res, **lo, **hi}
            if fitted:
                se = _calculate_sigma(y - mod["fitted"], len(y) - 1)
                res = _add_fitted_pi(res=res, se=se, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 468
class ARCH(GARCH):
    """Autoregressive Conditional Heteroskedasticity (ARCH) model.

    A particular case of the GARCH(p,q) model where $q=0$.
    It assumes that at time $t$, $y_t$ is given by:

    $$y_t = \epsilon_t \sigma_t$$

    with

    $$\sigma_t^2 = w0 + \sum_{i=1}^p a_i y_{t-i}^2$$.

    Here $\epsilon_t$ is a sequence of iid random variables with zero mean and unit variance.
    The coefficients $w$ and $a_i$, $i=1,...,p$ must be nonnegative and $\sum_{k=1}^p a_k < 1$.

    References
    ----------
    [Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica: Journal of the econometric society, 987-1007.](http://www.econ.uiuc.edu/~econ508/Papers/engle82.pdf)

    [James D. Hamilton. Time Series Analysis Princeton University Press, Princeton, New Jersey, 1st Edition, 1994.](https://press.princeton.edu/books/hardcover/9780691042893/time-series-analysis)

     Parameters
    ----------
    p : int
        Number of lagged versions of the series.
    alias : str
        Custom name of the model.
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        By default, the model will compute the native prediction
        intervals.
    """

    def __init__(
        self,
        p: int = 1,
        alias: str = "ARCH",
        prediction_intervals: Optional[ConformalIntervals] = None,
    ):
        self.p = p
        self.alias = alias
        super().__init__(p, q=0, alias=alias)

    def __repr__(self):
        return self.alias

# %% ../nbs/src/core/models.ipynb 479
class SklearnModel(_TS):
    """scikit-learn model wrapper

    Parameters
    ----------
    model : sklearn.base.BaseEstimator
        scikit-learn estimator
    prediction_intervals : Optional[ConformalIntervals]
        Information to compute conformal prediction intervals.
        This is required for generating future prediction intervals.
    alias : str, optional (default=None)
        Custom name of the model. If `None` will use the model's class.
    """

    def __init__(
        self,
        model,
        prediction_intervals: Optional[ConformalIntervals] = None,
        alias: Optional[str] = None,
    ):
        self.model = model
        self.prediction_intervals = prediction_intervals
        self.alias = alias

    def __repr__(self):
        if self.alias is not None:
            return self.alias
        return self.model.__class__.__name__

    def fit(
        self,
        y: np.ndarray,
        X: np.ndarray,
    ) -> "SklearnModel":
        """Fit the model.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Exogenous of shape (t, n_x).

        Returns
        -------
        self : SklearnModel
            Fitted SklearnModel object.
        """
        from sklearn.base import clone

        self.model_ = {"model": clone(self.model)}
        self.model_["model"].fit(X, y)
        self._store_cs(y=y, X=X)
        self.model_["fitted"] = self.model_["model"].predict(X)
        residuals = y - self.model_["fitted"]
        self.model_["sigma"] = _calculate_sigma(residuals, y.size)
        return self

    def predict(
        self,
        h: int,
        X: np.ndarray,
        level: Optional[List[int]] = None,
    ) -> Dict[str, Any]:
        """Predict with fitted SklearnModel.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Exogenous of shape (h, n_x).
        level: List[int]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"mean": self.model_["model"].predict(X)}
        if level is None:
            return res
        level = sorted(level)
        if self.prediction_intervals is not None:
            res = self._add_predict_conformal_intervals(res, level)
        else:
            raise Exception("You must pass `prediction_intervals` to compute them.")
        return res

    def predict_in_sample(self, level: Optional[List[int]] = None) -> Dict[str, Any]:
        """Access fitted SklearnModel insample predictions.

        Parameters
        ----------
        level : List[int]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        res = {"fitted": self.model_["fitted"]}
        if level is not None:
            level = sorted(level)
            res = _add_fitted_pi(res=res, se=self.model_["sigma"], level=level)
        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: np.ndarray,
        X_future: np.ndarray,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ) -> Dict[str, Any]:
        """Memory Efficient SklearnModel predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        h : int
            Forecast horizon.
        X : array-like
            Insample exogenous of shape (t, n_x).
        X_future : array-like
            Exogenous of shape (h, n_x).
        level : List[int]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        from sklearn.base import clone

        model = clone(self.model)
        model.fit(X, y)
        res = {"mean": model.predict(X_future)}
        if fitted:
            res["fitted"] = model.predict(X)
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                raise Exception("You must pass `prediction_intervals` to compute them.")
            if fitted:
                residuals = y - res["fitted"]
                sigma = _calculate_sigma(residuals, y.size)
                res = _add_fitted_pi(res=res, se=sigma, level=level)
        return res

    def forward(
        self,
        y: np.ndarray,
        h: int,
        X: np.ndarray,
        X_future: np.ndarray,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Apply fitted SklearnModel to a new/updated time series.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        h : int
            Forecast horizon.
        X : array-like
            Insample exogenous of shape (t, n_x).
        X_future : array-like
            Exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `constant` for point predictions and `level_*` for probabilistic predictions.
        """
        if not hasattr(self, "model_"):
            raise Exception("You have to use the `fit` method first")
        res = {"mean": self.model_["model"].predict(X_future)}
        if fitted:
            res["fitted"] = self.model_["model"].predict(X)
        if level is not None:
            level = sorted(level)
            if self.prediction_intervals is not None:
                res = self._add_conformal_intervals(fcst=res, y=y, X=X, level=level)
            else:
                raise Exception("You must pass `prediction_intervals` to compute them.")
            if fitted:
                se = _calculate_sigma(y - res["fitted"], y.size)
                res = _add_fitted_pi(res=res, se=se, level=level)
        return res

# %% ../nbs/src/core/models.ipynb 490
class ConstantModel(_TS):

    def __init__(self, constant: float, alias: str = "ConstantModel"):
        """Constant Model.

        Returns Constant values.

        Parameters
        ----------
        constant: float
            Custom value to return as forecast.
        alias: str
            Custom name of the model.
        """
        self.constant = constant
        self.alias = alias

    def __repr__(self):
        return self.alias

    def fit(
        self,
        y: np.ndarray,
        X: Optional[np.ndarray] = None,
    ):
        """Fit the Constant model.

        Fit an Constant Model to a time series (numpy.array) `y`.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (t, ).
        X : array-like
            Optional exogenous of shape (t, n_x).

        Returns
        -------
        self:
            Constant fitted model.
        """
        self.n_y = len(y)
        return self

    def predict(
        self,
        h: int,  # forecasting horizon
        X: Optional[np.ndarray] = None,  # exogenous regressors
        level: Optional[List[int]] = None,  # confidence level
    ):
        """Predict with fitted ConstantModel.

        Parameters
        ----------
        h : int
            Forecast horizon.
        X : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = np.full(h, self.constant, dtype=np.float32)
        res = {"mean": mean}

        if level is not None:
            for lv in sorted(level):
                res[f"lo-{lv}"] = mean
                res[f"hi-{lv}"] = mean

        return res

    def predict_in_sample(self, level: Optional[List[int]] = None):
        """Access fitted Constant Model insample predictions.

        Parameters
        ----------
        level : List[float]
            Confidence levels (0-100) for prediction intervals.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `fitted` for point predictions and `level_*` for probabilistic predictions.
        """
        fitted = np.full(self.n_y, self.constant, dtype=np.float32)
        res = {"fitted": fitted}
        if level is not None:
            for lv in sorted(level):
                res[f"fitted-lo-{lv}"] = fitted
                res[f"fitted-hi-{lv}"] = fitted

        return res

    def forecast(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Memory Efficient Constant Model predictions.

        This method avoids memory burden due from object storage.
        It is analogous to `fit_predict` without storing information.
        It assumes you know the forecast horizon in advance.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n,).
        h: int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels (0-100) for prediction intervals.
        fitted : bool
            Whether or not to return insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `mean` for point predictions and `level_*` for probabilistic predictions.
        """
        mean = np.full(h, self.constant, dtype=np.float32)
        res = {"mean": mean}

        if fitted:
            fitted_vals = np.full(self.n_y, self.constant, dtype=np.float32)
            res["fitted"] = fitted_vals

        if level is not None:
            for lv in sorted(level):
                res[f"lo-{lv}"] = mean
                res[f"hi-{lv}"] = mean
                if fitted:
                    res[f"fitted-lo-{lv}"] = fitted_vals
                    res[f"fitted-hi-{lv}"] = fitted_vals
        return res

    def forward(
        self,
        y: np.ndarray,
        h: int,
        X: Optional[np.ndarray] = None,
        X_future: Optional[np.ndarray] = None,
        level: Optional[List[int]] = None,
        fitted: bool = False,
    ):
        """Apply Constant model predictions to a new/updated time series.

        Parameters
        ----------
        y : numpy.array
            Clean time series of shape (n, ).
        h : int
            Forecast horizon.
        X : array-like
            Optional insample exogenous of shape (t, n_x).
        X_future : array-like
            Optional exogenous of shape (h, n_x).
        level : List[float]
            Confidence levels for prediction intervals.
        fitted : bool
            Whether or not returns insample predictions.

        Returns
        -------
        forecasts : dict
            Dictionary with entries `constant` for point predictions and `level_*` for probabilistic predictions.
        """
        res = self.forecast(
            y=y, h=h, X=X, X_future=X_future, level=level, fitted=fitted
        )
        return res

# %% ../nbs/src/core/models.ipynb 504
class ZeroModel(ConstantModel):

    def __init__(self, alias: str = "ZeroModel"):
        """Returns Zero forecasts.

        Returns Zero values.

        Parameters
        ----------
        alias: str
            Custom name of the model.
        """
        super().__init__(constant=0, alias=alias)

# %% ../nbs/src/core/models.ipynb 518
class NaNModel(ConstantModel):

    def __init__(self, alias: str = "NaNModel"):
        """NaN Model.

        Returns NaN values.

        Parameters
        ----------
        alias: str
            Custom name of the model.
        """
        super().__init__(constant=np.nan, alias=alias)