Multiple Seasonal-Trend decomposition with LOESS (MSTL)

Fast, effective forecasting of time series exhibiting multiple seasonality and trend.

Introduction

The MSTL algorithm, introduced in this paper, provides a way of applying Seasonal-Trend decomposition to multiple seasonalities. This allows effective modelling of time series with multiple complex components.

As well as the MSTL algorithm this crate also provides the MSTLModel struct, which is capable of running the MSTL algorithm over some time series data, then modelling the final trend component using a given trend forecaster. It can then recombine the forecasted trend with the decomposed seasonalities to generate in-sample and out-of-sample forecasts.

The latter use case is the main entrypoint of this crate.

Usage

use augurs::{mstl::MSTLModel, prelude::*};

# fn main() -> Result<(), Box<dyn std::error::Error>> {
// Input data must be a `&[f64]` for the MSTL algorithm.
let y = &[1.5, 3.0, 2.5, 4.2, 2.7, 1.9, 1.0, 1.2, 0.8];

// Define the number of seasonal observations per period.
// In this example we have hourly data and both daily and
// weekly seasonality.
let periods = vec![3, 4];
// Create an MSTL model using a naive trend forecaster.
// Note: in real life you may want to use a different
// trend forecaster - see below.
let mstl = MSTLModel::naive(periods);
// Fit the model. Note this consumes `mstl` and returns
// a fitted version.
let fit = mstl.fit(y)?;

// Obtain in-sample and out-of-sample predictions, along
// with prediction intervals.
let in_sample = fit.predict_in_sample(0.95)?;
let out_of_sample = fit.predict(10, 0.95)?;
# Ok(())
# }

Using alternative trend models

The MSTLModel is generic over the trend model used. As long as the model passed implements the TrendModel trait from this crate, it can be used to model the trend after decomposition. For example, the AutoETS struct from the ets crate can be used instead. First, add the augurs_ets crate to your Cargo.toml with the mstl feature enabled:

[dependencies]
augurs_ets = { version = "*", features = ["mstl"] }

use augurs_ets::AutoETS;
use augurs_mstl::MSTLModel;

let y = vec![1.5, 3.0, 2.5, 4.2, 2.7, 1.9, 1.0, 1.2, 0.8];

let periods = vec![24, 24 * 7];
let trend_model = AutoETS::non_seasonal();
let mstl = MSTLModel::new(periods, trend_model);
let fit = mstl.fit(y)?;

let in_sample = fit.predict_in_sample(0.95)?;
let out_of_sample = fit.predict(10, 0.95)?;

(Note that the above example doesn't compile for this crate due to a circular dependency, but would work in a separate crate!)

Implementing a trend model

To use your own trend model, you'll need a struct that implements the TrendModel trait. See below for an example of a model that predicts a constant value for all time points in the horizon.

use std::borrow::Cow;

use augurs_core::{Forecast, ForecastIntervals};
use augurs_mstl::{FittedTrendModel, TrendModel};

// The unfitted model. Sometimes this will need state!
#[derive(Debug)]
struct ConstantTrendModel {
    // The constant value to predict.
    constant: f64,
}

// The fitted model. This will invariable need state.
#[derive(Debug)]
struct FittedConstantTrendModel {
    // The constant value to predict.
    constant: f64,
    // The number of values in the training data.
    y_len: usize,
}

impl TrendModel for ConstantTrendModel {
    fn name(&self) -> Cow<'_, str> {
        Cow::Borrowed("Constant")
    }

    fn fit(
        &self,
        y: &[f64],
    ) -> Result<
        Box<dyn FittedTrendModel + Sync + Send>,
        Box<dyn std::error::Error + Send + Sync + 'static>,
    > {
        // Your algorithm should do whatever it needs to do to fit the model.
        // You need to return a boxed implementation of `FittedTrendModel`.
        Ok(Box::new(FittedConstantTrendModel {
            constant: self.constant,
            y_len: y.len(),
        }))
    }
}

impl FittedTrendModel for FittedConstantTrendModel {
    fn predict_inplace(
        &self,
        horizon: usize,
        level: Option<f64>,
        forecast: &mut Forecast,
    ) -> Result<(), Box<dyn std::error::Error + Send + Sync + 'static>> {
        forecast.point = vec![self.constant; horizon];
        if let Some(level) = level {
            let mut intervals = forecast
                .intervals
                .get_or_insert_with(|| ForecastIntervals::with_capacity(level, horizon));
            intervals.lower = vec![self.constant; horizon];
            intervals.upper = vec![self.constant; horizon];
        }
        Ok(())
    }

    fn predict_in_sample_inplace(
        &self,
        level: Option<f64>,
        forecast: &mut Forecast,
    ) -> Result<(), Box<dyn std::error::Error + Send + Sync + 'static>> {
        forecast.point = vec![self.constant; self.y_len];
        if let Some(level) = level {
            let mut intervals = forecast
                .intervals
                .get_or_insert_with(|| ForecastIntervals::with_capacity(level, self.y_len));
            intervals.lower = vec![self.constant; self.y_len];
            intervals.upper = vec![self.constant; self.y_len];
        }
        Ok(())
    }

    fn training_data_size(&self) -> Option<usize> {
        Some(self.y_len)
    }
}

Credits

This implementation is based heavily on both the R implementation and the statsforecast implementation. It also makes heavy use of the [stlrs][] crate.

References

Bandara, Kasun & Hyndman, Rob & Bergmeir, Christoph. (2021). “MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns”.

License

Dual-licensed to be compatible with the Rust project. Licensed under the Apache License, Version 2.0 <http://www.apache.org/licenses/LICENSE-2.0> or the MIT license <http://opensource.org/licenses/MIT>, at your option.