Multivariate Time Series Forecasting

• Multivariate time series forecasting is the process of predicting future values in a collection of interdependent variables by modeling both temporal and cross-sectional dependencies.
• Classical models like VARMA offer parameter interpretability, while machine learning approaches such as random forests and LSTM capture nonlinear patterns in complex data.
• Effective forecasting relies on robust feature engineering, including sliding windows, lagged variables, and careful handling of missing data.

Multivariate time series forecasting is the task of predicting future values in a collection of temporally indexed, interdependent variables, where each time point consists of a vector in ℝ^d rather than a scalar observation. This formulation generalizes univariate forecasting and is fundamental in domains ranging from econometrics to cyber-physical systems, open-source software health monitoring, and security risk management. The multivariate nature introduces additional challenges and opportunities in modeling cross-variable dependencies, handling missing data, and designing robust evaluation strategies.

1. Problem Definition and Formalization

A multivariate time series is a sequence $X = \{x_t\}_{t=1}^T$, where each $x_t \in \mathbb{R}^d$ represents the measurements of $d$ variables at timestamp $t$. Given a historical window $X_{1:W} = \{x_{t-W+1}, ..., x_t \}$, the forecasting task is to predict future observation(s) $Y = \{x_{t+h_1}, ..., x_{t+h_k}\}$ for one or more forecast horizons $H = \{h_1,...,h_k\}$. The forecast can target:

• The full future vector sequence (multi-step, multi-output)
• A subset of variables (targeted selection)
• Aggregated or transformed representations (e.g., trends, buckets, or derived indices)

Multivariate forecasting encompasses scenarios where dependencies exist both across time (temporal autocorrelation) and across variables (cross-correlation), demanding models that effectively capture and exploit both.

2. Classical Statistical Models

Classical approaches leverage parametric and semiparametric models that explicitly encode the temporal and cross-sectional dependencies.

• Vector Autoregressive Moving Average (VARMA): The prototypical linear model for stationary multivariate time series. VARMA(p, q) represents each variable as a linear function of past lags of all variables (autoregressive part) plus moving-average innovations:

$$y_t = \sum_{i=1}^p A_i y_{t-i} + \sum_{j=1}^q M_j \epsilon_{t-j} + \epsilon_t,$$

where $y_t \in \mathbb{R}^d$ and $A_i, M_j$ are $d \times d$ coefficient matrices (Tsakpinis et al., 26 Jan 2026).

• Feature Engineering for Statistical Methods: In practical VARMA application to software-maintenance forecasting, lagged aggregates and rolling statistics over windows are used to mimic the AR and MA components. Derived features such as cross-products (e.g., number of commits × number of issue comments) can further encode cross-dependencies (Tsakpinis et al., 26 Jan 2026).

These models are effective for stationary, low-noise signals and offer substantial interpretability through parameters and impulse-response analysis. However, their modeling flexibility is limited for highly nonlinear relationships or non-stationary dependencies.

3. Machine Learning and Deep Learning Models

To overcome the limitations of parametric models, machine learning and neural architectures have been applied to multivariate time series forecasting.

• Random Forests (RF): Non-parametric tree ensembles that ingest the flattened windowed feature matrix (using all variables and lags as features). They are robust to nonlinearity, variable interactions, and can perform both regression and classification tasks (e.g., predicting raw values, trend classes, or buckets) (Tsakpinis et al., 26 Jan 2026).
• LSTM and Sequence Models: Recurrent neural networks, particularly Long Short-Term Memory (LSTM) architectures, natively consume multivariate temporal sequences as tensor input (shape (W,F), with W = window length, F = number of features). Stacked LSTM layers with optional dense heads are typical. Training employs loss functions appropriate to the target (MSE for regression, cross-entropy for classification). Regularization (dropout) and early stopping are used to avoid overfitting (Tsakpinis et al., 26 Jan 2026).

An important empirical result is that, in certain contexts (e.g., monthly aggregated open-source maintenance indicators), simple RF and VARMA baselines can match or surpass deep learning models, indicating that model selection must consider the autocorrelation structure and variable interactions present in the data (Tsakpinis et al., 26 Jan 2026).

4. Forecast Target Representations and Evaluation Metrics

Multivariate forecasting targets admit diverse representations beyond direct numerical prediction, enhancing both accuracy and interpretability:

• Raw Score Prediction: Directly forecast the value vector $x_{t+h}$ (MSE loss).
• Bucketed/Discretized Levels: Predict discrete categories (e.g., [low, moderate, high]; cross-entropy loss), enabling coarse-grained risk assessment.
• Trend Slope and Type: Predict temporal derivatives, either numerical or discretized (e.g., upward, stable, downward trends; cross-entropy loss).

Aggregated and trend-based targets are often more forecastable (accuracy > 0.95 for bucketed classes in maintenance activity) than raw values due to the reduced impact of high-frequency noise (Tsakpinis et al., 26 Jan 2026).

Evaluation commonly uses accuracy (after discretization where needed), complemented by MAE, RMSE, precision, recall, and F1 scores. Sliding-window cross-validation with varying training-window lengths and forecast horizons is standard for empirical benchmarking.

Forecast Target	Typical Accuracy	Remarks
Raw Score	~0.78–0.80	Regression; hardest
Bucketed	>0.95	3-class
Trend Slope	~0.65–0.69	Regression (differences)
Trend Type	>0.80	3-class

Confusion matrices for discrete targets usually show most errors are between adjacent categories, indicating effective ordinal capture (Tsakpinis et al., 26 Jan 2026).

5. Feature Engineering and Data Preprocessing

Success in multivariate time series forecasting relies critically on well-chosen input features and data preprocessing protocols:

• Temporal Aggregation: Downsampling or aggregating (e.g., daily to monthly) can attenuate high-frequency noise and enhance regularity, especially when domain activity (such as open-source maintenance) is naturally bursty with low mean inter-event intervals (Tsakpinis et al., 26 Jan 2026).
• Sliding Windows and Lagged Features: Feature matrices are typically constructed by concatenating all variables over a fixed window of lags, possibly engineering rolling means, variances, or cross-feature products to augment the representation for models that expect tabular input (VARMA, RF) (Tsakpinis et al., 26 Jan 2026).
• Handling Missing Data and Activity Gates: In longitudinal studies, care is taken to define valid activity intervals (e.g., only compute scores for periods after the initial warmup window and before archival), with $g(t)$ gates for active/inactive status (Tsakpinis et al., 26 Jan 2026).
• Numerical Stability and Smoothing: Raw indicator variables such as commit counts or issue events may be highly variable; rolling aggregation (e.g., 90-day sliding window) is used for stability.

6. Empirical Findings and Best Practices

Case studies in OSS maintenance forecasting demonstrate that:

• Bucketed and categorical trend-type forecasts are substantially easier and more stable than predicting raw numerical scores, enabling high-confidence “traffic light” style operationalization in proactive risk dashboards (Tsakpinis et al., 26 Jan 2026).
• Simpler statistical/ML models (VARMA, RF) perform on par with more complex deep learning models (LSTM), with the latter showing little added value unless strong nonlinearity or intricate temporal patterns are present (Tsakpinis et al., 26 Jan 2026).
• Key sources of error include unpredictable exogenous shocks, idiosyncratic project behaviors, and incomplete activity indicators. Future work is needed for personalized models and integration of richer signals.

7. Applications and Limitations

Multivariate time series forecasting enables forward-looking risk and resource allocation in domains where entity state is inherently multidimensional and temporally dependent.

Notable applications:

• Proactive assessment of software maintenance risk in open-source dependency networks, where maintenance activity scores (%-activity, trend types) are derived from commits, issues, and comments over sliding windows (Tsakpinis et al., 26 Jan 2026).
• Infrastructure monitoring, economic and financial forecasting, sensor fusion in cyber-physical systems.

Limitations:

• Predictive accuracy is limited by the input feature set: omission of key signals (e.g., community communication outside core channels) can cap attainable performance.
• Structural breaks and regime changes (e.g., sudden abandonment or takeovers) are not easily forecastable from regular time series patterns (Tsakpinis et al., 26 Jan 2026).

A strategic approach is to tailor representation (target selection) and model class to the rate and regularity at which critical events occur in the system of study.

---

References

• "Forecasting the Maintained Score from the OpenSSF Scorecard for GitHub Repositories linked to PyPI libraries" (Tsakpinis et al., 26 Jan 2026)