PaperlayerView exact revision on arXiv

Source-linked AI summary

High-Dimensional Multivariate Forecasting with Low-Rank Gaussian Copula Processes

David Salinas, Michael Bohlke-Schneider, Laurent Callot, Roberto Medico, Jan Gasthaus

arXiv:1910.03002v2cs.LGstat.ML

TL;DR

High-dimensional multivariate forecasting needs to capture dependencies without the computational burden of full covariance estimation. The paper combines an RNN with a low-rank Gaussian copula process, reporting improved accuracy while scaling to much larger dimensions.

  • Problem

    Existing high-dimensional forecasting methods often assume independence because estimating large, time-varying covariance matrices is computationally and numerically difficult.

  • Method

    The model combines an LSTM-RNN with non-parametric marginal CDF transformations and a low-rank Gaussian copula covariance structure.

  • Results

    The method improves accuracy on most datasets, with CRPS and CRPS-Sum improving by average amounts of 10% and 40%, respectively, over the second-best models.

  • Takeaways & Limitations

    The approach supports probabilistic forecasting with time-varying covariance matrices for high-dimensional collections of time series.

  • Takeaways & Limitations

    The empirical marginal distributions assume stationarity; trends may require preprocessing such as detrending or differencing.

Abstract

from arXiv · show

Predicting the dependencies between observations from multiple time series is critical for applications such as anomaly detection, financial risk management, causal analysis, or demand forecasting. However, the computational and numerical difficulties of estimating time-varying and high-dimensional covariance matrices often limits existing methods to handling at most a few hundred dimensions or requires making strong assumptions on the dependence between series. We propose to combine an RNN-based time series model with a Gaussian copula process output model with a low-rank covariance structure to reduce the computational complexity and handle non-Gaussian marginal distributions. This permits to drastically reduce the number of parameters and consequently allows the modeling of time-varying correlations of thousands of time series. We show on several real-world datasets that our method provides significant accuracy improvements over state-of-the-art baselines and perform an ablation study analyzing the contributions of the different components of our model.

1 Introduction

High-dimensional probabilistic forecasting must represent dependencies among many heterogeneous time series, yet covariance estimation makes independence assumptions common and limits existing models. The paper proposes a low-rank Gaussian-copula approach that jointly models temporal dynamics and time-varying dependence while reducing parameter growth.

  • High-dimensional forecasting often assumes conditional independence because estimating large covariance matrices is difficult.
  • Independence assumptions are inadequate when correlations matter, including portfolio risk, retail competition, and multivariate anomaly detection.
  • Existing multivariate models remain unsuited to applications with more than a few hundred dimensions because their parameter counts grow rapidly with dimension.
  • The proposed model combines an LSTM-RNN, Gaussian copula processes, and low-rank-plus-diagonal covariance to jointly learn temporal dynamics and time-varying covariance with fewer parameters.
  • The method separately estimates each marginal distribution non-parametrically and uses the resulting CDF transformations in a Gaussian copula.
  • The paper reports scaling to dimensions up to an order of magnitude larger than previously reported and improved accuracy over state-of-the-art methods on six real-world datasets.

2 Autoregressive RNN Model for Probabilistic Multivariate Forecasting

The forecasting model uses tied, per-series LSTM state updates and a joint Gaussian-copula emission distribution. Its state-dependent mean and covariance produce autoregressive multivariate forecasts, with parameters learned by maximum likelihood.

  • Each series has an independently evolved LSTM state, while the transition parameters are shared across all time series.
  • The joint emission transforms every series through an invertible marginal mapping before modeling the transformed observations with a Gaussian distribution.
  • The functions µ(ht) and Σ(ht) map the collected state to the mean and covariance of the transformed observations.
  • The joint predictive distribution is factorized into sequential conditionals, allowing forecasts to be generated by sampling each time step and updating the states.
  • Model parameters are learned from observed data by maximum likelihood optimization of the forecasting loss.
  • During training, random fixed-size slices are sampled from long series, while prediction initializes states using only the final context window.

3 Gaussian Copula

The Gaussian copula separates marginal behavior from cross-series dependence by transforming each marginal to a uniform or Gaussian scale. Empirical CDF estimates make this construction applicable when marginal distributions are unknown.

  • A copula represents the joint distribution of variables with uniform marginal distributions.
  • Sklar’s theorem expresses a joint cumulative distribution as a copula applied to its univariate marginal CDFs.
  • For continuous marginals, the copula is uniquely determined by probability integral transforms, which are uniformly distributed.
  • The Gaussian copula applies the inverse standard-normal CDF to marginal CDF values and models the resulting variables with a multivariate normal distribution.
  • The method estimates unknown marginals using empirical CDFs and linearly interpolates them so the transformations remain differentiable.
  • The estimated marginal distributions use a fixed history of m = 100 past observations in the experiments.

4 Low-rank Gaussian Process Parametrization

The model uses a low-rank-plus-diagonal covariance and a Gaussian-process view of the emission distribution to make high-dimensional forecasting computationally manageable.

  • Emission distribution: The emission distribution is a Gaussian over marginally transformed observations, with mean and covariance determined from the LSTM state.Each transformation maps a series using its empirical CDF so that its transformed values are marginally standard normal.
  • Gaussian-process view: Viewing the emission distribution as a time-varying low-rank Gaussian process allows training on subsets of series in each mini-batch.This further alleviates memory constraints for very high-dimensional collections.
  • Covariance parametrization: The structured covariance Σ = D + V V^T uses O(N × r) parameters and evaluates the likelihood in O(Nr^2+r^3) operations.D is diagonal, V has dimensions N × r, and the experiments use r = 10.
  • Shared parametrization: Shared functions parameterize the mean, diagonal, and low-rank covariance components using each series’ state and feature vector.Features can be known a priori, learned embeddings, or a combination of both.

5 Experiments

Experiments evaluate covariance recovery and forecasting accuracy on synthetic and real-world multivariate datasets, including comparisons across architectures, distributions, and transformations.

  • Synthetic experiment: The synthetic experiment tests whether the method recovers complex time-varying low-rank covariance patterns.The data are generated from a normal distribution with time-varying mean and covariance.
  • Synthetic experiment: Figure 3 compares true and predicted covariance entries for N = 4 and N = 8 time series.The plots contain 10 lower-triangle entries including the diagonal for N = 4, and 28 for N = 8.
  • Evaluation design: The study compares correlated-sample baselines and ablations varying RNN architecture, output distribution, and data transformation.The transformations include copula, mean scaling, and no transformation; distributions include independent, low-rank, and full-rank normal outputs.
  • Results: 10% and 40% average improvements are reported for CRPS and CRPS-Sum, respectively, versus the second-best models.The comparison is against Vec-LSTM-lowrank-Copula for CRPS and Vec-LSTM-ind-scaling for CRPS-Sum.
  • Results: On traffic, GP models use 44K parameters versus 1.1M for Vec-LSTM with a low-rank distribution and 38M with a full-rank distribution.The paper also reports that GP-Copula significantly improves accuracy on most datasets.
  • Qualitative covariance analysis: Predicted covariance matrices reconstruct the evolving topology of spatial relationships across hours in the Taxi dataset.The covariance graph is obtained by retaining edges whose correlation exceeds a constant threshold.

6 Conclusion

The paper presents probabilistic forecasts for high-dimensional multivariate time series using low-rank covariance structure and a low-rank Gaussian copula process. It argues that estimating high-dimensional, time-varying covariance matrices may support applications including anomaly detection, imputation, and graph analysis.

  • Low-rank approximation avoids the potentially very large parameter count of full covariance matrices in high-dimensional forecasting.
  • A low-rank Gaussian copula process supports stable direct optimization of autoregressive-model parameters.
  • The proposed techniques may enable estimating high-dimensional time-varying covariance matrices for anomaly detection, imputation, and graph analysis.

Supplementary material

The supplementary material details the low-rank Gaussian likelihood, marginal CDF treatment, baseline models, training settings, and rank analysis. The low-rank covariance reduces parameter and likelihood-computation costs, while intermediate ranks can provide better test performance than higher ranks.

  • A Multivariate Likelihood: The covariance parametrization Σ = D + V V T uses diagonal D and low-rank V, with rank r much smaller than dimension d.
  • A Multivariate Likelihood: The restricted covariance has O(dr) parameters, and its Gaussian likelihood can be computed in O(dr2 + r3) time.
  • A Multivariate Likelihood: The likelihood calculation uses C = I_r + V T D^-1 V and the Woodbury identity to evaluate determinant and Mahalanobis terms efficiently.
  • Gaussian Copula: The truncated empirical CDF clips estimates to [δ_m, 1 − δ_m], with δ_m chosen to balance bias and variance; jitter is added for discrete data.
  • C Effect of rank on low-rank approximation: On electricity, training loss decreases with rank, but test performance is best at ranks 32/64; ranks 128 and 256 show a larger train–test gap consistent with over-fitting.

F Dataset details

The supplementary material describes dataset summaries, preprocessing, evaluation metrics, sampling, and parameter-comparison tables. Forecast accuracy is assessed both marginally and on sums across series using sample-based CRPS measures.

  • F Dataset details: Table 3 summarizes forecast horizon, domain, observation frequency, series dimension N, and available training time steps T for the datasets.
  • F Dataset details: Taxi preprocessing removes outliers, retains the most active spatial region, bins observations every 30 minutes, and aggregates locations on a spatial grid.
  • F Dataset details: Taxi uses January 2015 for training and January 2016 for testing.
  • Evaluation metrics: CRPS measures predictive-distribution compatibility and is a proper scoring rule minimized when the predictive distribution matches the data-generating distribution.
  • Evaluation metrics: Models generate 400 samples to estimate quantiles, and marginal CRPS is averaged across dimensions and forecast steps.
  • Evaluation metrics: CRPS-Sum evaluates the predictive distribution of the sum across dimensions, using summed samples before sorting into quantiles.

G.2 Mean Squared Error (MSE)

The MSE evaluates squared deviations between targets and predicted distribution means across all series and the complete prediction range. MSE-sum is defined analogously to CRPS-sum.

  • MSE averages squared errors over all N time series and the full prediction range from t = T − t0 + 1 through T.
  • MSE uses the target z and predicted distribution mean ẑ, while MSE-sum is defined analogously to CRPS-sum.

H Comparison with forecasting methods with diagonal covariance

The approach is evaluated against deep-learning forecasting baselines, while learned correlations are examined across datasets with differing dependence patterns.

  • Baselines: The evaluation compares the approach with DeepAR, MQCNN, and DeepGLO on selected datasets.DeepAR and MQCNN represent deep-learning forecasting baselines; DeepGLO is included on two author-provided datasets.
  • Baselines: None of the competing approaches models correlations across time series in their forecasts, while DeepGLO provides only point forecasts.
  • Learned correlations: The learned correlation matrices show nearly homoscedastic behavior for exchange-rate data, with most correlations close to 1.The currencies are measured relative to the US dollar, producing highly correlated series.
  • Learned correlations: The remaining datasets exhibit heteroscedastic correlations, including day-night patterns in solar, electricity, traffic, and taxi data.The model’s learned correlations are illustrated across all datasets in Figure 2.

J Effect of the number of evaluation samples on CRPS and inference runtime

Increasing the number of evaluation samples improves CRPS only modestly beyond 100 samples, while inference runtime increases linearly.

  • CRPS: More than 100 evaluation samples have only a small effect on CRPS.
  • Runtime: Inference runtime increases linearly as the number of evaluation samples grows.

K Additional experiments details

Additional experiments specify time features, lag inputs, implementation settings, and reproducibility details for the model evaluation.

  • Input features: Time is represented with calendar features matched to the dataset frequency and concatenated to the LSTM input at each time-step.Hourly, daily, and minute-granularity datasets use different calendar-feature combinations.
  • Evaluation: The GluonTS implementation is evaluated using CRPS-sum and CRPS accuracy metrics, with means and standard errors reported over three runs.
  • Input features: Lag inputs are selected according to time frequency, including [1, 24, 168] for hourly and for daily data.Thirty-minute data use [1, 2, 4, 12, 24, 48].
  • Implementation: Experiments run on an AWS c5.4xlarge instance with 16 cores and 32GB RAM, and all RNN models train and evaluate in under five hours.
  • Implementation: Training failures in Table 4 arise from out-of-memory errors or NaNs caused by unstable models.
  • Reproducibility: The model was re-implemented in GluonTS, with a static code version released to support reproducibility over time.The released version is separate from current GluonTS releases, whose reproducibility is not guaranteed.
Loading 1910.03002v2…