Image by Author | Piktochart

Multivariate time series analysis techniques are used to analyze different aspects of multiple time series, for instance, modeling interdependencies between them, identifying possible causal relationships such that one time series can “predict” another, and so on.

There are several techniques to analyze multiple time series, each one specialized in certain aspects. This article describes the practical application of two of them: Granger Causality to analyze causal or predictive relationships, and Vector Autoregression (VAR) to model interdependencies based on linear relationships.

Granger Causality Analysis

While we will use again the Chicago rides dataset as in previous parts of this series, in this occasion, we will extract a bivariate time series: in other words, data describing two time series along the same timeline: one for daily bus rides (‘bus‘), and one for daily train rides (‘rail_boardings‘).

We first load the dataset, preprocess it, and visualize both time series:

import numpy as np
import pandas as pd

url = "https://raw.githubusercontent.com/gakudo-ai/open-datasets/refs/heads/main/CTA_-_Ridership_-_Daily_Boarding_Totals.csv"
df = pd.read_csv(url)

df['service_date'] = pd.to_datetime(df['service_date'])
df.set_index('service_date', inplace=True)
df.sort_index(inplace=True)

data = df[['bus', 'rail_boardings']].dropna()

data.plot(figsize=(14, 7), title="Bus and Rail Ridership Over Time")

Granger causality analysis in Python is easily conducted using the grangercausalitytests() method, available after importing statsmodels.tsa.stattools.grangercausalitytests.

Before doing so, it is important to understand a key parameter in this method: the number of lags. Lags represent the number of previous time steps in one series used to try predicting the current value for another time series: for instance, to analyze whether recent bus ride figures up to five days ago can help predict today’s train rides or vice versa (five lags). This parameter is set with the maxlag argument, whose value is set as smaller to test short-term relationships between the time series, and larger for long-term relationships (the latter may introduce noise if set as too large).

Let’s apply the test with maxlag=5 and interpret the results:

from statsmodels.tsa.stattools import grangercausalitytests

maxlag = 5
print("Granger Causality Test (Do Bus Boardings Predict Rail Boardings?):")
grangercausalitytests(data[['rail_boardings', 'bus']], maxlag=maxlag, verbose=True)

print("\nGranger Causality Test (Do Rail Boardings Predict Bus Boardings?):")
grangercausalitytests(data[['bus', 'rail_boardings']], maxlag=maxlag, verbose=True)

Here’s an excerpt of the result (you’ll get a test result for each possible number of lags up to maxlag, and notice we are performing two sets of tests, for analyzing causality in both directions):

Granger Causality
number of lags (no zero) 1
ssr based F test:         F=162.1714, p=0.0000  , df_denom=8032, df_num=1
ssr based chi2 test:   chi2=162.2320, p=0.0000  , df=1
likelihood ratio test: chi2=160.6159, p=0.0000  , df=1
parameter F test:         F=162.1714, p=0.0000  , df_denom=8032, df_num=1

In essence, both large F-values and p-values close to zero (under 0.05) are indicators of a strong causal relationship. You’ll notice that in this example there is a very strong causal relationship between both time series: past bus rides are very useful in predicting future train rides, and vice versa.

Vector Autoregression

Vector autoregression (VAR) builds a machine learning model that captures interdependencies between our two time series, and can be used to forecast future rides based on such interdependencies. This code fits a VAR model, after which we forecast bus and train rides for the next 15 days.

from statsmodels.tsa.api import VAR
import matplotlib.pyplot as plt

model = VAR(data)
results = model.fit(maxlags=5)

forecast_steps = 15
forecast = results.forecast(data.values[-results.k_ar:], steps=forecast_steps)
forecast_index = pd.date_range(data.index[-1], periods=forecast_steps + 1, freq='D')[1:]
forecast_df = pd.DataFrame(forecast, index=forecast_index, columns=['bus_forecast', 'rail_forecast'])

plt.figure(figsize=(10, 6))
data_2020 = data.loc['2022-01-01':]
plt.plot(data_2020['bus'], label="Bus (Actual)", alpha=0.7)
plt.plot(data_2020['rail_boardings'], label="Rail (Actual)", alpha=0.7)
plt.plot(forecast_df['bus_forecast'], label="Bus (Forecast)", linestyle="--", color='green')
plt.plot(forecast_df['rail_forecast'], label="Rail (Forecast)", linestyle="--", color='red')

plt.title("VAR Forecast for Bus and Rail Ridership (2022)")
plt.xlabel("Date")
plt.ylabel("Ridership")
plt.legend()
plt.grid()
plt.show()