Multivariate Time Series Forecasting

Multivariate time series forecasting works similarly to univariate time series forecasting (covered here and here). The main difference is that you must specify the index of a target univariate to forecast, e.g. for a 5-variable time series you may want to forecast the value of the 3rd variable (we specify this by indicating target_seq_index = 2). To begin, we will load the multivariate SeattleTrail dataset for time series forecasting.

from merlion.utils import TimeSeries
from ts_datasets.forecast import SeattleTrail

time_series, metadata = SeattleTrail()[0]
train_data = TimeSeries.from_pd(time_series[metadata["trainval"]])
test_data = TimeSeries.from_pd(time_series[~metadata["trainval"]])

print(f"Time series is {train_data.dim}-dimensional")

Time series is 5-dimensional

Model Initialization and Training

For the purposes of this tutorial, we will be using 3 models:

1. DefaultForeacster (which automatically detects whether the input time series is univariate or multivariate);

2. ARIMA (a classic univariate algorithm) trained to forecast a specific univariate; and

3. A ForecasterEnsemble which selects the better of the two models.

All models are trained with a maximum allowed forecasting horizon of 100 steps. Note that all multivariate forecasting models can be used for univariate time series, and by specifying target_seq_index appropriately, univariate models can be used for multivariate time series as well. Moreover, the API is identical in all cases.

from merlion.evaluate.forecast import ForecastMetric
from merlion.models.factory import ModelFactory
from merlion.models.ensemble.combine import ModelSelector

# Time series is sampled hourly, so max_forecast_steps = 24 means we can predict up to 1 day in the future
target_seq_index = 2
max_forecast_steps = 24
kwargs = dict(target_seq_index=target_seq_index, max_forecast_steps=max_forecast_steps)

model1 = ModelFactory.create("DefaultForecaster", **kwargs)
model2 = ModelFactory.create("Arima", **kwargs)

# This ModelSelector combiner picks the best model based on sMAPE
model3 = ModelFactory.create("ForecasterEnsemble", models=[model1, model2],
                             combiner=ModelSelector(metric=ForecastMetric.sMAPE))

for model in [model1, model2, model3]:
    print(f"Training {type(model).__name__}...")
    train_pred, train_stderr = model.train(train_data)

Training DefaultForecaster...
Training Arima...
Training ForecasterEnsemble...

Validation: 100%|██████████| 31550400/31550400 [00:20<00:00, 1558063.14it/s]

Model Inference and Quantitative Evaluation

Like univariate models, we may call model.forecast() to get a forecast and potentially a standard error for the model. We can use these to evaluate the model’s performance. Note that the model selector successfully picks the better of the two models.

from merlion.evaluate.forecast import ForecastMetric

target_univariate = test_data.univariates[test_data.names[target_seq_index]]
target = target_univariate[:max_forecast_steps].to_ts()

for model in [model1, model2, model3]:
    forecast, stderr = model.forecast(target.time_stamps)
    rmse = ForecastMetric.RMSE.value(ground_truth=target, predict=forecast)
    smape = ForecastMetric.sMAPE.value(ground_truth=target, predict=forecast)
    print(f"{type(model).__name__}")
    print(f"RMSE:  {rmse:.4f}")
    print(f"sMAPE: {smape:.4f}")
    print()

DefaultForecaster
RMSE:  6.6216
sMAPE: 121.1709

Arima
RMSE:  10.2208
sMAPE: 140.2772

ForecasterEnsemble
RMSE:  6.6216
sMAPE: 121.1709

We can also use a ForecastEvaluator to evaluate a model in a manner that simulates live deployment. Here, we train an initial model on the training data, and we obtain its predictions on the training data using a sliding window of 1 day (horizon="1d" means that we want the model to predict 1 day in the future at each time step, and cadence="1d" means that we obtain a prediction from the model once per day). Note that we never actually re-train the model (retrain_freq=None).

from merlion.evaluate.forecast import ForecastEvaluator, ForecastEvaluatorConfig

for model in [model1, model2]:
    print(f"{type(model).__name__} Sliding Window Evaluation")
    evaluator = ForecastEvaluator(model=model, config=ForecastEvaluatorConfig(
        horizon="1d", cadence="1d", retrain_freq=None))
    train_result, test_pred = evaluator.get_predict(train_vals=train_data, test_vals=test_data)
    rmse = evaluator.evaluate(ground_truth=test_data, predict=test_pred, metric=ForecastMetric.RMSE)
    smape = evaluator.evaluate(ground_truth=test_data, predict=test_pred, metric=ForecastMetric.sMAPE)
    print(f"RMSE:  {rmse:.4f}")
    print(f"sMAPE: {smape:.4f}")
    print()

DefaultForecaster Sliding Window Evaluation

ForecastEvaluator: 100%|██████████| 31528800/31528800 [01:22<00:00, 383184.32it/s]

RMSE:  12.0339
sMAPE: 99.4165

Arima Sliding Window Evaluation

ForecastEvaluator: 100%|██████████| 31528800/31528800 [03:45<00:00, 139586.71it/s]

RMSE:  13.1032
sMAPE: 112.2607