End to End Walkthrough with Polars

​

Introducing Polars: A High-Performance DataFrame Library

1. Performance Efficiency: Constructed using Rust, Polars exhibits an exemplary ability to manage substantial datasets with remarkable speed and minimal memory usage.
2. Lazy Evaluation: Polars operates on the principle of ‘lazy evaluation’, creating an optimized logical plan of operations for efficient execution, a feature that mirrors the functionality of Apache Spark.
3. Parallel Execution: Demonstrating the capability to exploit multi-core CPUs, Polars facilitates parallel execution of operations, substantially accelerating data processing tasks.

Prerequesites This Guide assumes basic familiarity with StatsForecast. For a minimal example visit the Quick Start

1. Install packages.
2. Read the data.
3. Explore the data.
4. Train many models for every unique combination of time series.
5. Evaluate the model’s performance using cross-validation.
6. Select the best model for every unique time series.

Not Covered in this guide

• Forecasting at scale using clusters on the cloud.
  • Forecast the M5 Dataset in 5min using Ray clusters.
  • Forecast the M5 Dataset in 5min using Spark clusters.
  • Learn how to predict 1M series in less than 30min.
• Training models on Multiple Seasonalities.
  • Learn to use multiple seasonality in this Electricity Load forecasting tutorial.
• Using external regressors or exogenous variables
  • Follow this tutorial to include exogenous variables like weather or holidays or static variables like category or family.
• Comparing StatsForecast with other popular libraries.
  • You can reproduce our benchmarks here.

​

Install libraries

​

Read the data

• The unique_id (string, int or category) represents an identifier for the series.
• The ds (datestamp or int) column should be either an integer indexing time or a datestampe ideally like YYYY-MM-DD for a date or YYYY-MM-DD HH:MM:SS for a timestamp.
• The y (numeric) represents the measurement we wish to forecast.

import polars as pl

Y_df = pl.read_parquet('https://datasets-nixtla.s3.amazonaws.com/m4-hourly.parquet')
Y_df.head()

Note Processing time is dependent on the available computing resources. Running this example with the complete dataset takes around 10 minutes in a c5d.24xlarge (96 cores) instance from AWS.

uids = Y_df['unique_id'].unique(maintain_order=True)[:10] # Select 10 ids to make the example faster
Y_df = Y_df.filter(pl.col('unique_id').is_in(uids))
Y_df = Y_df.group_by('unique_id').tail(7 * 24) #Select last 7 days of data to make example faster

​

Explore Data with the plot method

Note The StatsForecast.plot method uses matplotlib as a default engine. You can change to plotly by setting engine="plotly".

from statsforecast import StatsForecast

StatsForecast.plot(Y_df)

​

Train multiple models for many series

• Auto Forecast: Automatic forecasting tools search for the best parameters and select the best possible model for a series of time series. These tools are useful for large collections of univariate time series. Includes automatic versions of: Arima, ETS, Theta, CES.
• Exponential Smoothing: 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. Examples: SES, Holt’s Winters, SSO.
• Benchmark models: classical models for establishing baselines. Examples: Mean, Naive, Random Walk
• Intermittent or Sparse models: suited for series with very few non-zero observations. Examples: CROSTON, ADIDA, IMAPA
• Multiple Seasonalities: suited for signals with more than one clear seasonality. Useful for low-frequency data like electricity and logs. Examples: MSTL.
• Theta Models: fit two theta lines to a deseasonalized time series, using different techniques to obtain and combine the two theta lines to produce the final forecasts. Examples: Theta, DynamicTheta

• AutoARIMA: Automatically selects the best ARIMA (AutoRegressive Integrated Moving Average) model using an information criterion. Ref: AutoARIMA.
• HoltWinters: triple exponential smoothing, Holt-Winters’ method is an extension of exponential smoothing for series that contain both trend and seasonality. Ref: HoltWinters
• SeasonalNaive: Memory Efficient Seasonal Naive predictions. Ref: SeasonalNaive
• HistoricAverage: arthimetic mean. Ref: HistoricAverage.
• DynamicOptimizedTheta: The theta family of models has been shown to perform well in various datasets such as M3. Models the deseasonalized time series. Ref: DynamicOptimizedTheta.

from statsforecast.models import (
    HoltWinters,
    CrostonClassic as Croston, 
    HistoricAverage,
    DynamicOptimizedTheta as DOT,
    SeasonalNaive
)

# Create a list of models and instantiation parameters
models = [
    HoltWinters(),
    Croston(),
    SeasonalNaive(season_length=24),
    HistoricAverage(),
    DOT(season_length=24)
]

• models: a list of models. Select the models you want from models and import them.
• freq: a string indicating the frequency of the data. (See panda’s available frequencies.) This is also available with Polars.
• n_jobs: n_jobs: int, number of jobs used in the parallel processing, use -1 for all cores.
• fallback_model: a model to be used if a model fails.

# Instantiate StatsForecast class as sf
sf = StatsForecast( 
    models=models,
    freq=1, 
    n_jobs=-1,
    fallback_model=SeasonalNaive(season_length=7),
    verbose=True,
)

Note StatsForecast achieves its blazing speed using JIT compiling through Numba. The first time you call the statsforecast class, the fit method should take around 5 seconds. The second time -once Numba compiled your settings- it should take less than 0.2s.

• h (int): represents the forecast h steps into the future. In this case, 12 months ahead.
• level (list of floats): this optional parameter is used for probabilistic forecasting. Set the level (or confidence percentile) of your prediction interval. For example, level=[90] means that the model expects the real value to be inside that interval 90% of the times.

Note The forecast method is compatible with distributed clusters, so it does not store any model parameters. If you want to store parameters for every model you can use the fit and predict methods. However, those methods are not defined for distrubed engines like Spark, Ray or Dask.

forecasts_df = sf.forecast(df=Y_df, h=48, level=[90])
forecasts_df.head()

Forecast:   0%|          | 0/10 [Elapsed: 00:00]

sf.plot(Y_df,forecasts_df)

# Plot to unique_ids and some selected models
sf.plot(Y_df, forecasts_df, models=["HoltWinters","DynamicOptimizedTheta"], unique_ids=["H10", "H105"], level=[90])

# Explore other models 
sf.plot(Y_df, forecasts_df, models=["SeasonalNaive"], unique_ids=["H10", "H105"], level=[90])

​

Evaluate the model’s performance

Tip Setting n_windows=1 mirrors a traditional train-test split with our historical data serving as the training set and the last 48 hours serving as the testing set.

• df: training data frame
• h (int): represents h steps into the future that are being forecasted. In this case, 24 hours ahead.
• step_size (int): step size between each window. In other words: how often do you want to run the forecasting processes.
• n_windows(int): number of windows used for cross validation. In other words: what number of forecasting processes in the past do you want to evaluate.

cv_df = sf.cross_validation(
    df=Y_df,
    h=24,
    step_size=24,
    n_windows=2
)

• unique_id: series identifier
• ds: datestamp or temporal index
• cutoff: the last datestamp or temporal index for the n_windows. If n_windows=1, then one unique cuttoff value, if n_windows=2 then two unique cutoff values.
• y: true value
• "model": columns with the model’s name and fitted value.

cv_df.head()

from utilsforecast.losses import mse

def evaluate_cv(df, metric):
    models = [c for c in df.columns if c not in ('unique_id', 'ds', 'cutoff', 'y')]
    evals = mse(cv_df, models=models)
    pos2model = dict(enumerate(models))
    return evals.with_columns(
        best_model=pl.concat_list(models).list.arg_min().replace_strict(pos2model)
    )

Warning You can also use Mean Average Percentage Error (MAPE), however for granular forecasts, MAPE values are extremely hard to judge and not useful to assess forecasting quality.

evaluation_df = evaluate_cv(cv_df, mse)
evaluation_df.head()

evaluation_df['best_model'].value_counts()

seasonal_ids = evaluation_df.filter(pl.col('best_model') == 'SeasonalNaive')['unique_id']
sf.plot(Y_df,forecasts_df, unique_ids=seasonal_ids, models=["SeasonalNaive","DynamicOptimizedTheta"])

​

Select the best model for every unique series

def get_best_model_forecast(forecasts_df, evaluation_df):
    models = {
        c.replace('-lo-90', '').replace('-hi-90', '')
        for c in forecasts_df.columns
        if c not in ('unique_id', 'ds')
    }
    model2pos = {m: i for i, m in enumerate(models)}
    with_best = forecasts_df.join(evaluation_df[['unique_id', 'best_model']], on='unique_id')
    return with_best.select(
        'unique_id',
        'ds',
        *[
            (
                pl.concat_list([f'{m}{suffix}' for m in models])
                .list.get(pl.col('best_model').replace_strict(model2pos))
                .alias(f'best_model{suffix}')
            )
            for suffix in ('', '-lo-90', '-hi-90')
        ]
    )

prod_forecasts_df = get_best_model_forecast(forecasts_df, evaluation_df)
prod_forecasts_df.head()

sf.plot(Y_df, prod_forecasts_df, level=[90])