Add other models from the paper?

[amyheather]

https://github.com/TomMonks/swast-benchmarking

Naive benchmark: Seasonal Naive method - using the most recent observation from the same day of the week and carrying it forward (e.g., if forecasting next Tuesday, use the observation from the most recent Tuesday as the predicted value).

Candidate models:

1. Holt-Winters Exponential Smoothing
2. Automatic Autoregressive Integrated Moving Average (autoARIMA)
3. Lagged regression (autoregression) with holidays and seasonal indexes
4. Lagged regularised regression (elastic-net) with seasonal indexes
5. Regression with holidays and ARIMA errors
6. Regression with holidays, seasonal indexes and ARIMA errors
7. Harmonic Regression (fourier terms) with holidays
8. Singular Spectrum Analysis
9. Facebook Prophet
10. Trigonometric seasonality, Box-Cox transformation (TBATS)
11. Comb: Simple Exponential Smoothing, Linear Trend, damped trend
12. Ensemble of:
    • Holt-Winters Exponental Smoothing (1).
    • Regression with holiday and ARIMA errors (5).
13. Ensemble of:
    • Regression with holiday and ARIMA errors (5).
    • Facebook Prophet (9).
14. Ensemble of:
    • Holt-Winters Exponential Smoothing (1).
    • Regression with holidays and ARIMA errors (5).
    • Facebook Prophet (9)

There appears to be some errors in the table in the paper. The above list of candidates was corrected based on the descriptions from the notebooks. In the paper, it says model 12 is actually 1+2 (not 1+5), and that model 13 is 5+7 (not 5+9).

[amyheather]

I add seasonal naive, which is good to have I think.

I also started adding exponential smoothing - but then, wasn't sure if it's actually the best use of time implementing them all? If know what came out best last time? Although... that was all pre-COVID.... so maybe still worth revisiting?

"""Holt-Winters Exponential Smoothing with additive errors and seasonality."""

from dataclasses import dataclass

import pandas as pd
from statsmodels.tsa.statespace.exponential_smoothing import (
    ExponentialSmoothing,
)

from .structures import CustomRepr


@dataclass(kw_only=True, repr=False)
class ExponentialSmoothingParams(CustomRepr):
    """Parameters for the exponential smoothing model.

    Parameters
    ----------
    trend
    damped_trend
    seasonal
    interval_width : float
        Width of the prediction intervals - for example, 0.95 will produce
        95% prediction intervals.

    """

    trend: bool = True
    damped_trend: bool = True
    seasonal: int = 7
    interval_width: float = 0.95


def exponential_smoothing(train, params, test, horizon):
    """Fit exponential smoothing model and generate forecast.

    Parameters
    ----------
    train : pd.DataFrame
        Historic training data.
    params : ExponentialSmoothingParams
        Parameters controlling the model.
    test : pd.DataFrame | None
        Data containing the dates to forecast in a `ds` column. Typically
        a held-out test set or cross-validation fold.
    horizon : int | None
        Number of days to forecast, after the final date in `train`. To be
        used if there is no test set (e.g., if actually predicting into future
        with no data to compare against).

    Returns
    -------
    pd.DataFrame
        Forecast dataframe.

    """
    if (test is None) == (horizon is None):
        raise ValueError("Provide exactly one of 'test' or 'horizon'.")

    # Convert training data into required format
    train = train.sort_values("ds").set_index("ds")["y"]

    # Fit model
    model = ExponentialSmoothing(
        endog=train,
        trend=params.trend,
        damped_trend=params.damped_trend,
        seasonal=params.seasonal,
    )
    fitted = model.fit()

    # Get dates for forecast
    if test is not None:
        forecast_dates = test["ds"].reset_index(drop=True)
    else:
        forecast_dates = pd.date_range(
            start=train.index.max() + pd.Timedelta(days=1),
            periods=horizon,
            freq="D",
        )

    pred = fitted.get_forecast(len(forecast_dates))
    pred.summary_frame(alpha=1 - params.interval_width)
    print(pred)

[amyheather]

In discussion with Tom, he felt this is likely not worth pursueing.