SCALECAST 0.11.0

The practitioner's time series forecasting library

🌄 Scalecast: The practitioner's time series forecasting library

About

Scalecast is a light-weight modeling procedure, wrapper, and results container meant for those who are looking for the fastest way possible to apply, tune, and validate many different model classes for forecasting applications. In the Data Science industry, it is often asked of practitioners to deliver predictions and ranges of predictions for several lines of businesses or data slices, 100s or even 1000s. In such situations, it is common to see a simple linear regression or some other quick procedure applied to all lines due to the complexity of the task. This works well enough for people who need to deliver something, but more can be achieved.

The scalecast package was designed to address this situation and offer advanced machine learning models that can be applied, optimized, and validated quickly. Unlike many libraries, the predictions produced by scalecast are always dynamic by default, not averages of one-step forecasts, so you don't run into the situation where the estimator looks great on the test-set but can't generalize to real data. What you see is what you get, with no attempt to oversell results. If you download a library that looks like it's able to predict the COVID pandemic in your test-set, you probably have a one-step forecast happening under-the-hood. You can't predict the unpredictable, and you won't see such things with scalecast.

The library provides the Forecaster (for one series) and MVForecaster (for multiple series) wrappers around the following estimators:

• Any regression model from Sklearn, including Sklearn APIs (like Xgboost and LightGBM)
• Recurrent neural nets from Keras TensorFlow
• Classic econometric models from statsmodels: Holt-Winters Exponential Smoothing and ARIMA
• Facebook Prophet
• LinkedIn Silverkite
• Average, weighted average, and spliced models

Installation

• Only the base package is needed to get started:
  pip install scalecast
• Optional add-ons:
  pip install fbprophet (prophet model--see here to resolve a common installation issue if using Anaconda)
  pip install greykite (silverkite model)
  pip install tqdm (progress bar with notebook)
  pip install ipython (widgets with notebook)
  pip install ipywidgets (widgets with notebook)
  jupyter nbextension enable --py widgetsnbextension (widgets with notebook)
  jupyter labextension install @jupyter-widgets/jupyterlab-manager (widgets with Lab)

Links

Links
📚 Read the Docs	Official scalecast docs
📋 Examples	Official scalecast notebooks
📓 TDS Article 1	Univariate Forecasting
📓 TDS Article 2	Multivariate Forecasting
📓 TDS Article 3	Feature Reduction
🛠️ Change Log	See what's changed

Example

Let's say we wanted to forecast each of the 1-year, 5-year, 10-year, 20-year, and 30-year corporate bond rates through the next 12 months. There are two ways we could do this with scalecast:

1. Forecast each series individually (univariate): they will only have their own histories and any exogenous regressors we add to make forecasts. One series is predicted forward dynamically at a time.
2. Forecast all series together (multivariate): they will have their own histories, exogenous regressors, and each other's histories to make forecasts. All series will be predicted forward dynamically at the same time.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scalecast.Forecaster import Forecaster
from scalecast.MVForecaster import MVForecaster
from scalecast import GridGenerator
from scalecast.multiseries import export_model_summaries
import pandas_datareader as pdr

sns.set(rc={'figure.figsize':(14,7)})

df = pdr.get_data_fred(
    ['HQMCB1YR','HQMCB5YR','HQMCB10YR','HQMCB20YR','HQMCB30YR'],
    start='2000-01-01',
    end='2022-03-01'
)

f_dict = {c: Forecaster(y=df[c],current_dates=df.index) for c in df}

Option 1 - Univariate

Select Models

models = (
    'arima',       # linear time series model
    'elasticnet',  # linear model with regularization
    'knn',         # nearest neighbor model
    'xgboost',     # boosted tree model
)

Create Grids

These grids will be used to tune each model. To get example grids, you can use:

GridGenerator.get_example_grids()

This saves a Grids.py file to your working directory by default, which scalecast knows how to read. In this example, we create our own grids:

arima_grid = dict(
    order = [(1,1,1),(0,1,1),(0,1,0)],
    seasonal_order = [(1,1,1,12),(0,1,1,12),(0,1,0,12)],
    Xvars = [None,'all']
)

elasticnet_grid = dict(
    l1_ratio = [0.25,0.5,0.75,1],
    alpha = np.linspace(0,1,100),
)

knn_grid = dict(
    n_neighbors = np.arange(2,100,2)
)

xgboost_grid = dict(
     n_estimators=[150,200,250],
     scale_pos_weight=[5,10],
     learning_rate=[0.1,0.2],
     gamma=[0,3,5],
     subsample=[0.8,0.9],
)

grid_list = [arima_grid,elasticnet_grid,knn_grid,xgboost_grid]
grids = dict(zip(models,grid_list))

Select test length, validation length, and forecast horizon

def prepare_fcst(f):
    f.set_test_length(0.2)
    f.set_validation_length(12)
    f.generate_future_dates(12)

Add seasonal regressors

These are regressors like month, quarter, dayofweek, dayofyear, minute, hour, etc. Raw integer values, dummy variables, or fourier transformed variables. They are determined by the series' own histories.

def add_seasonal_regressors(f):
    f.add_seasonal_regressors('month',raw=False,sincos=True)
    f.add_seasonal_regressors('year')
    f.add_seasonal_regressors('quarter',raw=False,dummy=True,drop_first=True)

Choose Autoregressive Terms

A better way to do this would be to examine each series individually for autocorrelation. This example uses three lags for each series and one seasonal seasonal lag (assuming 12-month seasonality).

def add_ar_terms(f):
    f.add_ar_terms(3)       # lags
    f.add_AR_terms((1,12))  # seasonal lags

Write the forecast procedure

def tune_test_forecast(k,f,models):
    for m in models:
        print(f'forecasting {m} for {k}')
        f.set_estimator(m)
        f.ingest_grid(grids[m])
        f.tune()
        f.auto_forecast()

Run a forecast loop

for k, f in f_dict.items():
    prepare_fcst(f)
    add_seasonal_regressors(f)
    add_ar_terms(f)
    f.integrate(critical_pval=0.01) # takes differences in series until they are stationary using the adf test
    tune_test_forecast(k,f,models)

forecasting arima for HQMCB1YR
forecasting elasticnet for HQMCB1YR
forecasting knn for HQMCB1YR
forecasting xgboost for HQMCB1YR
forecasting arima for HQMCB5YR
forecasting elasticnet for HQMCB5YR
forecasting knn for HQMCB5YR
forecasting xgboost for HQMCB5YR
forecasting arima for HQMCB10YR
forecasting elasticnet for HQMCB10YR
forecasting knn for HQMCB10YR
forecasting xgboost for HQMCB10YR
forecasting arima for HQMCB20YR
forecasting elasticnet for HQMCB20YR
forecasting knn for HQMCB20YR
forecasting xgboost for HQMCB20YR
forecasting arima for HQMCB30YR
forecasting elasticnet for HQMCB30YR
forecasting knn for HQMCB30YR
forecasting xgboost for HQMCB30YR

Visualize results

Since there are 5 series to visualize, it might be undesirable to write a plot function for each one. Instead, scalecast lets you leverage Jupyter widgets by using this function:

from scalecast.notebook import results_vis
results_vis(f_dict)

Because we aren't able to show widgets through markdown, this readme shows a visualization for the 30-year rate only:

Integrated Results

f.plot_test_set(ci=True,order_by='LevelTestSetMAPE')
plt.title(f'{k} test-set results',size=16)
plt.show()

Level Results

f.plot_test_set(level=True,order_by='LevelTestSetMAPE')
plt.title(f'{k} test-set results',size=16)
plt.show()

Comparing level test-set MAPE values, the K-nearest Neighbor model performed best, although, as we can see, predicting bond rates accurately is difficult if not impossible. To make the forecasts look better, we can set dynamic_testing=False in the manual_forecast() or auto_forecast() methods when calling forecasts. This will make test-set predictions an average on one-step forecasts. By default, everything is dynamic with scalecast to give a more realistic sense of how the models perform. To see our future predictions:

f.plot(level=True,models='knn')
plt.title(f'{k} forecast',size=16)
plt.show()

View Results

We can print a dataframe that shows how each model performed on each series.

results = export_model_summaries(f_dict,determine_best_by='LevelTestSetMAPE')
results.columns

Index(['ModelNickname', 'Estimator', 'Xvars', 'HyperParams', 'Scaler',
       'Observations', 'Tuned', 'DynamicallyTested', 'Integration',
       'TestSetLength', 'TestSetRMSE', 'TestSetMAPE', 'TestSetMAE',
       'TestSetR2', 'LastTestSetPrediction', 'LastTestSetActual', 'CILevel',
       'CIPlusMinus', 'InSampleRMSE', 'InSampleMAPE', 'InSampleMAE',
       'InSampleR2', 'ValidationSetLength', 'ValidationMetric',
       'ValidationMetricValue', 'models', 'weights', 'LevelTestSetRMSE',
       'LevelTestSetMAPE', 'LevelTestSetMAE', 'LevelTestSetR2', 'best_model',
       'Series'],
      dtype='object')

results[['ModelNickname','Series','LevelTestSetMAPE','LevelTestSetR2','HyperParams']]

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

</style>

	ModelNickname	Series	LevelTestSetMAPE	LevelTestSetR2	HyperParams
0	elasticnet	HQMCB1YR	3.178337	-1.156049	{'l1_ratio': 0.25, 'alpha': 0.0}
1	knn	HQMCB1YR	3.526606	-1.609881	{'n_neighbors': 6}
2	arima	HQMCB1YR	5.052915	-4.458788	{'order': (1, 1, 1), 'seasonal_order': (0, 1, ...
3	xgboost	HQMCB1YR	5.881190	-6.700988	{'n_estimators': 150, 'scale_pos_weight': 5, '...
4	xgboost	HQMCB5YR	0.372265	0.328466	{'n_estimators': 250, 'scale_pos_weight': 5, '...
5	knn	HQMCB5YR	0.521959	0.119923	{'n_neighbors': 26}
6	elasticnet	HQMCB5YR	0.664711	-0.263214	{'l1_ratio': 0.25, 'alpha': 0.0}
7	arima	HQMCB5YR	1.693632	-7.335685	{'order': (1, 1, 1), 'seasonal_order': (0, 1, ...
8	elasticnet	HQMCB10YR	0.145834	0.390825	{'l1_ratio': 0.25, 'alpha': 0.010101010101010102}
9	knn	HQMCB10YR	0.175513	0.341443	{'n_neighbors': 26}
10	xgboost	HQMCB10YR	0.465610	-2.875923	{'n_estimators': 150, 'scale_pos_weight': 5, '...
11	arima	HQMCB10YR	0.569411	-4.968624	{'order': (1, 1, 1), 'seasonal_order': (0, 1, ...
12	elasticnet	HQMCB20YR	0.096262	0.475912	{'l1_ratio': 0.25, 'alpha': 0.030303030303030304}
13	knn	HQMCB20YR	0.103565	0.486593	{'n_neighbors': 26}
14	xgboost	HQMCB20YR	0.105995	0.441079	{'n_estimators': 200, 'scale_pos_weight': 5, '...
15	arima	HQMCB20YR	0.118033	0.378309	{'order': (1, 1, 1), 'seasonal_order': (0, 1, ...
16	knn	HQMCB30YR	0.075318	0.560975	{'n_neighbors': 22}
17	elasticnet	HQMCB30YR	0.089038	0.538360	{'l1_ratio': 0.25, 'alpha': 0.030303030303030304}
18	xgboost	HQMCB30YR	0.098148	0.495318	{'n_estimators': 200, 'scale_pos_weight': 5, '...
19	arima	HQMCB30YR	0.099816	0.498638	{'order': (1, 1, 1), 'seasonal_order': (0, 1, ...

Option 2: Multivariate

Select Models

Only sklearn models are available with multivariate forecasting, so we can replace ARIMA with mlr.

mv_models = (
    'mlr',
    'elasticnet',
    'knn',
    'xgboost',
)

Create Grids

We can use three of the same grids as we did in univariate forecasting and create a new MLR grid, with a modification to also search the optimal lag numbers. The lags argument can be an int, list, or dict type and all series will use the other series' lags (as well as their own lags) in each model that is called, unless we tell the models to not use a certain series' lags for whatever reason. Again, for mv forecasting, we can save default grids:

GridGenerator.get_mv_grids()

This creates the MVGrids.py file in our working directory by default, which scalecast knows how to read.

mlr_grid = dict(lags = np.arange(1,13,1))
elasticnet_grid['lags'] = np.arange(1,13,1)
knn_grid['lags'] = np.arange(1,13,1)
xgboost_grid['lags'] = np.arange(1,13,1)

mv_grid_list = [mlr_grid,elasticnet_grid,knn_grid,xgboost_grid]
mv_grids = dict(zip(mv_models,mv_grid_list))

Create multivariate forecasting object

• Need to change test and validation length
• Regressors are already carried forward from the underlying Forecaster objects
• Integrated levels are also carried forward from the underlying Forecaster objects

mvf = MVForecaster(
    *f_dict.values(),
    names = f_dict.keys(),
)

mvf.set_test_length(.2)
mvf.set_validation_length(12)
mvf

MVForecaster(
    DateStartActuals=2000-02-01T00:00:00.000000000
    DateEndActuals=2022-03-01T00:00:00.000000000
    Freq=MS
    N_actuals=266
    N_series=5
    SeriesNames=['HQMCB1YR', 'HQMCB5YR', 'HQMCB10YR', 'HQMCB20YR', 'HQMCB30YR']
    ForecastLength=12
    Xvars=['monthsin', 'monthcos', 'year', 'quarter_2', 'quarter_3', 'quarter_4']
    TestLength=53
    ValidationLength=12
    ValidationMetric=rmse
    ForecastsEvaluated=[]
    CILevel=0.95
    BootstrapSamples=100
    CurrentEstimator=mlr
    OptimizeOn=mean
)

Choose how to optimize the models when tuning hyperparameters

Default behavior is use the mean performance of each model on all series. We don't have to run the line below to keep this behavior, but we also have the option to use this code to optimize on the min/max error across all series, a weighted average of the error across the series, or to only consider one series' error over all others. See the docs for more info.

mvf.set_optimize_on('mean')

Write Forecasting Procedure

• Instead of grid search, we will use randomized grid search to speed up evaluation times

for m in mv_models:
    print(f'forecasting {m}')
    mvf.set_estimator(m)
    mvf.ingest_grid(mv_grids[m])
    mvf.limit_grid_size(100,random_seed=20) # do this because now grids are larger and this speeds it up
    mvf.tune()
    mvf.auto_forecast()

forecasting mlr
forecasting elasticnet
forecasting knn
forecasting xgboost

Set best model

mvf.set_best_model(determine_best_by='LevelTestSetMAPE')
mvf.best_model

'knn'

The elasticnet model was chosen based on its average test-set MAPE performance on all series.

Visualize results

Multivariate forecasting allows us to view all series and all models together. This could get jumbled, so let's just see the mlr and elasticnet results, knowing we can see the others if we want later.

Integrated Results

mvf.plot_test_set(ci=True,models=['mlr','elasticnet'])
plt.title(f'test-set results',size=16)
plt.show()

Level Results

mvf.plot_test_set(level=True,models=['mlr','elasticnet'])
plt.title(f'test-set results',size=16)
plt.show()

Once again, in this object, we can also set dynamic_testing=False in the manual_forecast() or auto_forecast() methods when calling forecasts. This would make our plots and metrics look better, but not be realistic in the sense it wouldn't show how well our forecasts could predict more than one period in the future. Let's see model forecasts into the future using the elasticent model only:

mvf.plot(level=True,models='elasticnet')
plt.title(f'forecasts',size=16)
plt.show()

View Results

We can print a dataframe that shows how each model did on each series.

mvresults = mvf.export_model_summaries()
mvresults.columns

Index(['Series', 'ModelNickname', 'Estimator', 'Xvars', 'HyperParams', 'Lags',
       'Scaler', 'Observations', 'Tuned', 'DynamicallyTested', 'Integration',
       'TestSetLength', 'TestSetRMSE', 'TestSetMAPE', 'TestSetMAE',
       'TestSetR2', 'LastTestSetPrediction', 'LastTestSetActual', 'CILevel',
       'CIPlusMinus', 'InSampleRMSE', 'InSampleMAPE', 'InSampleMAE',
       'InSampleR2', 'ValidationSetLength', 'ValidationMetric',
       'ValidationMetricValue', 'LevelTestSetRMSE', 'LevelTestSetMAPE',
       'LevelTestSetMAE', 'LevelTestSetR2', 'OptimizedOn', 'MetricOptimized',
       'best_model'],
      dtype='object')

mvresults[['ModelNickname','Series','LevelTestSetMAPE','LevelTestSetR2','HyperParams','Lags']]

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

</style>

	ModelNickname	Series	LevelTestSetMAPE	LevelTestSetR2	HyperParams	Lags
0	knn	HQMCB1YR	2.324660	-0.246684	{'n_neighbors': 38}	1
1	mlr	HQMCB1YR	3.793842	-2.039209	{}	1
2	elasticnet	HQMCB1YR	3.793793	-2.039144	{'l1_ratio': 1.0, 'alpha': 0.0}	1
3	xgboost	HQMCB1YR	5.647091	-6.032171	{'n_estimators': 200, 'scale_pos_weight': 10, ...	6
4	knn	HQMCB5YR	0.576072	0.015072	{'n_neighbors': 38}	1
5	mlr	HQMCB5YR	0.824692	-0.856614	{}	1
6	elasticnet	HQMCB5YR	0.824665	-0.856494	{'l1_ratio': 1.0, 'alpha': 0.0}	1
7	xgboost	HQMCB5YR	0.940324	-1.357810	{'n_estimators': 200, 'scale_pos_weight': 10, ...	6
8	knn	HQMCB10YR	0.210301	0.166048	{'n_neighbors': 38}	1
9	mlr	HQMCB10YR	0.243091	-0.063537	{}	1
10	elasticnet	HQMCB10YR	0.243075	-0.063448	{'l1_ratio': 1.0, 'alpha': 0.0}	1
11	xgboost	HQMCB10YR	0.229634	-0.576659	{'n_estimators': 200, 'scale_pos_weight': 10, ...	6
12	knn	HQMCB20YR	0.110740	0.436453	{'n_neighbors': 38}	1
13	mlr	HQMCB20YR	0.107826	0.438056	{}	1
14	elasticnet	HQMCB20YR	0.107812	0.438077	{'l1_ratio': 1.0, 'alpha': 0.0}	1
15	xgboost	HQMCB20YR	0.288137	-3.211295	{'n_estimators': 200, 'scale_pos_weight': 10, ...	6
16	knn	HQMCB30YR	0.087734	0.568717	{'n_neighbors': 38}	1
17	mlr	HQMCB30YR	0.084045	0.558262	{}	1
18	elasticnet	HQMCB30YR	0.084037	0.558233	{'l1_ratio': 1.0, 'alpha': 0.0}	1
19	xgboost	HQMCB30YR	0.300581	-4.218186	{'n_estimators': 200, 'scale_pos_weight': 10, ...	6

Backtest results

To test how well, on average, our models would have done across the last-10 12-month forecast horizons, we can use the backtest() method. It works for both the Forecaster and MVForecaster objects.

mvf.backtest('elasticnet')
mvf.export_backtest_metrics('elasticnet')

<style scoped> .dataframe tbody tr th:only-of-type { vertical-align: middle; }

.dataframe tbody tr th {
    vertical-align: top;
}

.dataframe thead th {
    text-align: right;
}

</style>

		iter1	iter2	iter3	iter4	iter5	iter6	iter7	iter8	iter9	iter10	mean
series	metric
HQMCB1YR	RMSE	0.487971	0.295814	0.147322	0.164049	0.091141	0.060353	0.119651	0.028003	0.046719	0.154311	0.159533
MAE	0.28571	0.171958	0.097406	0.136998	0.070702	0.056345	0.117123	0.02111	0.041883	0.149434	0.114867
R2	-0.074922	-0.032213	0.237629	-1.07959	-0.758305	-2.192486	-35.041028	-0.680347	-2.329534	-13.099094	-5.504989
MAPE	0.408617	0.31018	0.237814	0.426764	0.239001	0.226237	0.475545	0.082124	0.166775	0.593249	0.316631
HQMCB5YR	RMSE	0.570817	0.65871	0.604627	0.644545	0.537364	0.410313	0.401821	0.458499	0.408862	0.273161	0.496872
MAE	0.402543	0.525583	0.522105	0.585369	0.466181	0.344805	0.326988	0.392991	0.327104	0.203996	0.409767
R2	-0.205337	-2.315526	-4.08979	-7.256139	-5.021421	-3.150129	-3.676984	-4.558652	-3.030571	-0.79983	-3.410438
MAPE	0.213215	0.316168	0.353514	0.43293	0.355474	0.274756	0.269099	0.342425	0.28483	0.178973	0.302138
HQMCB10YR	RMSE	0.374449	0.474401	0.507182	0.635606	0.534922	0.4349	0.483086	0.586451	0.589008	0.352687	0.497269
MAE	0.310979	0.392149	0.464092	0.601646	0.481887	0.376152	0.403262	0.512031	0.498309	0.277614	0.431812
R2	-0.109954	-3.014724	-7.328628	-12.560674	-5.475112	-2.715311	-3.472925	-4.833621	-4.089653	-0.683906	-4.428451
MAPE	0.111612	0.140797	0.175868	0.234974	0.18859	0.148315	0.160139	0.207382	0.201964	0.112856	0.168250
HQMCB20YR	RMSE	0.370336	0.255213	0.318524	0.493896	0.413632	0.283031	0.379169	0.523698	0.605537	0.303488	0.394652
MAE	0.33989	0.18065	0.280991	0.471872	0.375782	0.245164	0.316363	0.465334	0.539688	0.251894	0.346763
R2	-0.556491	-0.335237	-2.450294	-7.18817	-3.561754	-0.947127	-2.450003	-5.43214	-6.657672	-0.731047	-3.030993
MAPE	0.105202	0.053226	0.086066	0.147459	0.116833	0.076122	0.097674	0.144917	0.16848	0.079162	0.107514
HQMCB30YR	RMSE	0.37143	0.217251	0.267355	0.453609	0.397683	0.244522	0.370656	0.506118	0.607924	0.297473	0.373402
MAE	0.34647	0.176251	0.223837	0.431921	0.360314	0.206274	0.31087	0.450336	0.546683	0.252607	0.330556
R2	-0.810715	0.070295	-0.935815	-4.488289	-2.754592	-0.381914	-2.206821	-5.173769	-7.18541	-0.80564	-2.467267
MAPE	0.106906	0.053282	0.066859	0.131948	0.109309	0.062304	0.093255	0.135931	0.165438	0.076923	0.100215

Correlation Matrices

• If you want to see how correlated the series are in your MVForecaster object, you can use these correlation matrices

All Series, no lags

heatmap_kwargs = dict(
    disp='heatmap',
    vmin=-1,
    vmax=1,
    annot=True,
    cmap = "Spectral",
)
mvf.corr(**heatmap_kwargs)
plt.show()

Two series, with lags

mvf.corr_lags(y='HQMCB1YR',x='HQMCB30YR',lags=12,**heatmap_kwargs)
plt.show()

There is much more than can be done with this package! Be sure to read the docs and see the examples!