Multivariate Time Series Forecasting using Python

In today's data-driven world, predicting the future behavior of complex systems often requires understanding how multiple variables interact over time. Multivariate time series forecasting allows us to simultaneously predict several interdependent variables, offering richer insights than traditional univariate methods.

Introduction

Time series analysis is a cornerstone of predictive analytics, enabling us to forecast future trends based on historical data. While univariate models focus on a single variable, many real-world phenomena are influenced by a web of interconnected factors. Consider economic indicators, stock markets, or even the performance metrics in a complex AI system like my own GrowthAI project; predicting one without accounting for others often leads to suboptimal results. This is where multivariate time series forecasting shines, leveraging the relationships between variables to produce more accurate and holistic predictions. Instead of treating each variable in isolation, it recognizes that their movements are often correlated, with changes in one potentially signaling changes in others.

This article will delve into the realm of multivariate time series forecasting using Python. We'll explore core concepts, practical implementation steps, and various modeling techniques, from classical statistical methods like Vector Autoregression (VAR) to modern deep learning approaches such as Long Short-Term Memory (LSTM) networks. My goal is to provide a comprehensive guide, backed by concrete examples, to equip you with the knowledge to tackle your own interdependent forecasting challenges.

Why It Matters

The ability to predict multiple variables simultaneously holds immense value across countless domains. In finance, forecasting the closing prices of several correlated stocks, like AAPL, MSFT, NFLX, and GOOG, provides a more complete market picture than isolated predictions. For supply chain management, anticipating demand for various product lines while also considering raw material prices can optimize inventory and production schedules. Similarly, in an AI system like GrowthAI, where we track multiple key performance indicators (KPIs) like user engagement, conversion rates, and retention, multivariate forecasting helps us understand the collective health and trajectory of the system, identifying cascading effects that univariate models would miss.

The primary advantage of multivariate forecasting is its capacity to capture the intricate dependencies and interactions between variables. When variables exhibit strong correlations or Granger causality, ignoring these relationships can lead to less accurate predictions. By modeling these interdependencies, multivariate approaches can often achieve superior prediction accuracy compared to a collection of independent univariate forecasts. This holistic view is particularly crucial in dynamic environments where a single event can ripple across an entire system, allowing us to anticipate not just what will happen to one variable, but how that event will affect related outcomes. It allows Kishna to make more informed decisions by understanding the systemic impact rather than isolated trends.

Core Concepts

Before diving into implementation, understanding a few core concepts is essential. At its heart, multivariate time series forecasting involves predicting multiple time-dependent variables concurrently using their historical values and the historical values of other related series. Unlike univariate models that forecast a single variable, this approach explicitly models the relationships between all variables in the system.

One critical concept is stationarity. A stationary time series has statistical properties (mean, variance, autocorrelation) that do not change over time. Many statistical forecasting models, including VAR, assume stationarity. If a series is non-stationary, as is often the case with financial data like stock prices, it typically needs to be transformed, most commonly through differencing. Differencing involves calculating the difference between consecutive observations (e.g., $Y_t - Y_{t-1}$). As we'll see, the Augmented Dickey-Fuller (ADF) test is a standard tool for assessing stationarity.

Another key aspect is lag selection. In models like Vector Autoregression (VAR), the number of past observations (lags) used to predict future values is crucial. Criteria like the Akaike Information Criterion (AIC) or Bayesian Information Criterion (BIC) help determine the optimal lag order by balancing model fit with complexity. For evaluating forecasts, common metrics include Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and Mean Absolute Percentage Error (MAPE), which quantify the difference between predicted and actual values on a hold-out test set.

Architecture and How It Works

The typical workflow for multivariate time series forecasting is a structured process designed to ensure robust and accurate predictions. As KishnaKushwaha would approach any data-intensive project, it begins with meticulous data handling and preparation:

Step-by-Step Implementation

Let's walk through a practical implementation of multivariate time series forecasting using a Vector Autoregression (VAR) model in Python. We will use a dataset of OHLCV (Open, High, Low, Close, Volume) data for four major stocks: AAPL, MSFT, NFLX, and GOOG. The dataset covers approximately three months, from February 7, 2023, to May 5, 2023, with 62 daily observations per stock and no missing values. You can download this ideal dataset from here (Coming Soon)" target="_blank">this link and save it as stocks_ohlcv.csv in your project directory.

Data Loading and Initial Exploration

First, we'll load the data and pivot it to get the closing prices for each stock in a wide format, which is suitable for multivariate analysis. We'll use pandas for data manipulation and statsmodels for our VAR model.

import pandas as pd
import numpy as np
from statsmodels.tsa.api import VAR
from statsmodels.tsa.stattools import adfuller
import matplotlib.pyplot as plt

# Load data
df = pd.read_csv('stocks_ohlcv.csv', parse_dates=['Date'])
# Keep closing prices for the four stocks
close = df.pivot(index='Date', columns='Ticker', values='Close')

print("First 5 rows of closing prices:")
print(close.head())
print(f"Dataset shape: {close.shape}")
        

Checking for Stationarity and Differencing

Next, we perform the Augmented Dickey-Fuller (ADF) test on each stock's closing price series to check for stationarity. As previously noted, the raw closing prices are non-stationary, so we'll apply first-order differencing to achieve stationarity. This transformation is crucial for the validity of the VAR model's assumptions.

# Check stationarity
print("\nADF Test p-values for raw closing prices:")
for col in close.columns:
    adf = adfuller(close[col].dropna())
    print(f'{col}: ADF p-value = {adf[1]:.3f}')

# Difference to achieve stationarity
close_diff = close.diff().dropna()

print("\nFirst 5 rows of differenced closing prices:")
print(close_diff.head())

# Re-check stationarity on differenced data
print("\nADF Test p-values for differenced closing prices:")
for col in close_diff.columns:
    adf = adfuller(close_diff[col].dropna())
    print(f'{col}: ADF p-value = {adf[1]:.3f}')
        

After differencing, the ADF test p-values for all series will be significantly less than 0.05, confirming stationarity. This transformation stabilizes the mean of the time series, making it suitable for modeling with VAR.

Model Selection and Training (VAR)

With stationary data, we can now fit the VAR model. A crucial step here is selecting the optimal lag order. The select_order method from statsmodels helps us determine this by evaluating various information criteria, such as the Akaike Information Criterion (AIC). Once the optimal lag is identified, we fit the VAR model to our differenced data.

# Select lag order using AIC
model = VAR(close_diff)
lag_selection = model.select_order(maxlags=15) # Test up to 15 lags
optimal_lag = lag_selection.aic # Using AIC as per common practice and KB
print(f'\nOptimal lag (AIC): {optimal_lag}')

# Fit VAR model with optimal lag
var_fit = model.fit(optimal_lag)
print("\nVAR Model Summary:")
print(var_fit.summary())
        

Forecasting and Inverse Transformation

After fitting, we generate our 5-day ahead forecasts. Since our model was trained on differenced data, the forecasts will also be in terms of differences. To get the actual price levels, we need to inverse-transform these forecasts by cumulatively adding them to the last known closing prices. This step restores the predictions to their original scale, making them interpretable.

# Forecast next 5 days
forecast_steps = 5
# The .y attribute holds the fitted values. We pass the *last* known data point for forecasting.
forecast = var_fit.forecast(y=var_fit.y, steps=forecast_steps)

# Create a DataFrame for the differenced forecasts
last_date = close_diff.index[-1]
forecast_index = pd.date_range(start=last_date + pd.Timedelta(days=1), periods=forecast_steps, freq='B') # 'B' for business day frequency
forecast_df = pd.DataFrame(forecast, index=forecast_index, columns=close_diff.columns)

# Undo differencing to get price level forecasts
last_known_prices = close.iloc[-1] # Last known actual closing prices
forecast_level = last_known_prices + forecast_df.cumsum()

print("\nForecasted Closing Prices (USD) for next 5 business days:")
print(forecast_level)
        

The model generated 5-day ahead forecasts for May 6-10, 2023. For AAPL, the forecasted closing prices were approximately 174.37, 169.68, 168.19, 160.84, 167.65 USD. Similar forecasts were produced for GOOG, MSFT, and NFLX, providing a forward-looking view for each intertwined stock.

Visualizing the Forecast

Finally, we visualize the historical closing prices alongside our 5-day forecast. This provides an intuitive check on the model's performance and allows for visual validation, ensuring the forecasts make sense in the context of past trends.

# Plot historical and forecasted prices
plt.figure(figsize=(14, 7))
ax = close.plot(label='Historical', ax=plt.gca(), linewidth=1.5)
forecast_level.plot(ax=ax, label='Forecast', linestyle='--', linewidth=2)
plt.title('VAR Forecast of Stock Closing Prices (AAPL, MSFT, NFLX, GOOG)', fontsize=16)
plt.xlabel('Date', fontsize=12)
plt.ylabel('Closing Price (USD)', fontsize=12)
plt.legend(title='Ticker', bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(True, linestyle=':', alpha=0.7)
plt.tight_layout()
plt.show()
        

This visualization allows us to see how the predicted prices align with the recent historical trajectory of the stocks. Such plots are essential for sanity-checking any forecasting model, especially when dealing with critical financial data.

After executing the script, you'll observe how the forecasted lines extend from the historical data, providing a visual representation of the predicted future movements. This graphical output, as Kishna has found invaluable in projects like Intervu for understanding user engagement trends, offers immediate insights into the model's projections.

Practical Examples

Example 1: Stock Price Forecasting with VAR

The complete step-by-step implementation above serves as a robust example of forecasting daily closing prices of AAPL, MSFT, NFLX, and GOOG jointly using a VAR model on differenced series. This scenario is highly relevant for financial analysts, portfolio managers, or anyone needing to understand the collective movements of a set of interdependent assets. The VAR model's strength lies in its ability to explicitly model the linear relationships between these stocks, recognizing that a significant movement in one stock might predict subsequent movements in others.

For instance, if AAPL has a strong positive day, the VAR model might capture a historical tendency for MSFT to also perform well shortly after. By differencing the data, we ensure stationarity, making the statistical properties of the series stable over time, which is a key assumption for VAR. The AIC-based lag selection optimizes the model's ability to capture these past dependencies without overfitting. The resulting 5-day forecasts provide actionable insights, such as AAPL forecasted to be around $174.37 on the first day of the forecast horizon, followed by GOOG at $113.62, MSFT at $320.18, and NFLX at $327.60. These joint predictions are far more insightful than individual forecasts.

Example 2: Multivariate LSTM for Complex Patterns

While VAR models are excellent for linear relationships, many real-world time series exhibit complex, non-linear patterns. This is where deep learning models like Long Short-Term Memory (LSTM) networks come into play. LSTMs are a type of recurrent neural network (RNN) particularly suited for sequence prediction tasks because they can learn long-term dependencies. In my work on Intervu, forecasting user behavior metrics based on complex interaction sequences, LSTMs proved invaluable.

Consider a scenario where we want to predict the next day's OHLCV values for multiple stocks simultaneously, incorporating not just closing prices but also open, high, low, and volume from the past 30 days. This would be a multivariate, multi-step forecasting problem. The LSTM can process multiple input features (all OHLCV columns for all stocks) over a sequence (30 days) and output predictions for multiple features for the next day. The provided snippet demonstrates the basic structure for a multivariate LSTM model using Keras.

# Multivariate LSTM example (Keras)
import numpy as np
import pandas as pd
from sklearn.preprocessing import MinMaxScaler
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense, Reshape
import tensorflow as tf # Added for Reshape layer

# Prepare multivariate sequence
def create_multivariate_sequences(data, n_steps_in, n_steps_out):
    X, y = [], []
    for i in range(len(data)):
        end_ix = i + n_steps_in
        out_end_ix = end_ix + n_steps_out
        if out_end_ix > len(data):
            break
        X.append(data[i:end_ix, :])
        y.append(data[end_ix:out_end_ix, :])
    return np.array(X), np.array(y)

# Load and scale closing prices (using simplified data for demonstration)
# In a real scenario, you'd use all OHLCV features and concatenate them.
close_data_for_lstm = pd.read_csv('stocks_ohlcv.csv', parse_dates=['Date']).pivot(index='Date', columns='Ticker', values='Close')

scaler = MinMaxScaler()
scaled = scaler.fit_transform(close_data_for_lstm.values)

n_in, n_out = 3, 1 # Example: use past 3 days to predict next 1 day
X, y = create_multivariate_sequences(scaled, n_in, n_out)

# Train/test split
split = int(0.8 * len(X))
X_train, X_test = X[:split], X[split:]
y_train, y_test = y[:split], y[split:]

# Build LSTM
model = Sequential()
model.add(LSTM(64, activation='relu', input_shape=(n_in, X.shape[2])))
model.add(Dense(n_out * close_data_for_lstm.shape[1]))
model.add(Reshape((n_out, close_data_for_lstm.shape[1])))
model.compile(optimizer='adam', loss='mse')

# Fit the model (epochs reduced for quick execution)
model.fit(X_train, y_train, epochs=10, batch_size=1, validation_data=(X_test, y_test), verbose=0)

# Forecast
pred_scaled = model.predict(X_test[-1:])  # shape (1, n_out, n_features)
pred = scaler.inverse_transform(pred_scaled[0])
print('\nMultivariate LSTM - Next-day forecast for all stocks (example):', pred)
        

This code illustrates how to structure a multivariate input for an LSTM and predict multiple outputs. LSTMs are particularly powerful for scenarios involving long sequences and where the relationships between variables are intricate and non-linear. The Reshape layer is crucial for ensuring the output matches the desired (n_out, n_features) format for inverse scaling. If you're exploring this area, check out AI interview prep tool — Try Intervu free.

Example 3: State-Space Models for Dynamic Systems

State-space models, also known as Dynamic Linear Models (DLMs) or Kalman filters, offer a flexible framework for modeling dynamic systems with unobserved components. They are particularly useful when dealing with measurement noise, missing data, and systems where underlying 'states' evolve over time. These models represent the observed series as functions of hidden states, which themselves follow a dynamic process. This approach is powerful for capturing the joint evolution of variables like electricity demand and temperature, or in my Voice Agent project, understanding the latent intent and sentiment evolving during a conversation.

The statsmodels library provides implementations of state-space models, including VARMAX, which can be configured to act as a state-space model. For instance, we could apply a Dynamic Linear Model to capture the joint evolution of electricity demand and temperature across hourly intervals, allowing for real-time updates as new data arrives. The following snippet shows a basic setup using VARMAX as a proxy for a state-space model:

# State-Space (Dynamic Linear Model) example using statsmodels
import pandas as pd
import numpy as np
import statsmodels.api as sm

# Assume df has columns: demand, temperature
# For this example, we'll create synthetic data. 
# In a real scenario, you'd load actual energy.csv or similar.
np.random.seed(42)
num_obs = 100
dates = pd.date_range(start='2023-01-01', periods=num_obs, freq='H')
df_energy = pd.DataFrame({
    'demand': np.random.normal(loc=100, scale=10, size=num_obs) + np.sin(np.arange(num_obs)/10) * 20,
    'temperature': np.random.normal(loc=20, scale=5, size=num_obs) + np.cos(np.arange(num_obs)/15) * 8
}, index=dates)

# Build a VARMAX model (which includes state-space representation)
# order=(1, 0) specifies a VAR(1) model without moving average terms.
# You can extend this with exogenous variables (X) or seasonal components.
model_ss = sm.tsa.VARMAX(df_energy, order=(1, 0)) 
results_ss = model_ss.fit(disp=False)

# Forecast next 12 periods
forecast_ss = results_ss.get_forecast(steps=12)
pred_mean_ss = forecast_ss.predicted_mean
pred_ci_ss = forecast_ss.conf_int()

print('\nState-Space Model (VARMAX) - Forecasted Mean for next 12 periods:')
print(pred_mean_ss)
print('\nState-Space Model (VARMAX) - 95% Confidence Intervals:')
print(pred_ci_ss)
        

This example demonstrates how to set up a basic state-space model for joint forecasting. State-space models are particularly powerful because they can incorporate different types of components (trend, seasonality, regression, noise) into a single, coherent framework, providing both point forecasts and uncertainty estimates (confidence intervals).

Conclusion

Multivariate time series forecasting is an indispensable tool in the modern data scientist's toolkit, particularly when dealing with systems where multiple variables interact and influence each other. As Kishna has observed in the development of projects like GrowthAI and Voice Agent, recognizing and modeling these interdependencies is key to unlocking deeper insights and making more robust predictions. From understanding stock market dynamics to optimizing resource allocation, the ability to forecast multiple series concurrently provides a significant advantage over univariate approaches.

We've explored the fundamental concepts, walked through a step-by-step implementation of a VAR model in Python for stock price forecasting, and briefly touched upon advanced techniques like LSTMs and State-Space models. By applying techniques such as stationarity testing, differencing, and careful model selection, you can build powerful predictive models that capture the complex dance between your data points. The journey into multivariate forecasting is rich and rewarding, offering endless opportunities to uncover hidden patterns and drive more intelligent decision-making.

Explore more technical guides and tutorials on our articles page.