Modelling multivariate ordinal time series using pairwise likelihood

1. Problem Statement

The paper addresses the challenge of jointly modeling multiple ordinal time series. While univariate ordinal time series modeling is well-established (using Markov models, regression, latent continuous models, etc.), extending these to a multivariate setting is difficult due to two primary factors:

1. Complexity of Joint Distribution: Constructing a multivariate distribution for ordinal variables that is both flexible and computationally tractable is inherently difficult.
2. Dependence Structure: A valid multivariate framework must simultaneously account for:
   • Autocorrelation: Serial dependence within each individual time series.
   • Cross-correlation: Dependence between different time series at the same time point (contemporaneous) and across lags.

Traditional approaches often rely on high-dimensional copulas to link marginal distributions. However, calculating the joint probability mass function for a $K$-dimensional ordinal vector requires evaluating $2^K$ terms, leading to severe computational bottlenecks as the number of series ($K$) increases.

2. Methodology

The authors propose a Pairwise Likelihood approach combined with Copulas to circumvent the computational intractability of full multivariate likelihoods.

A. Model Specification

• Marginal Models: Each ordinal time series $Z_{kt}$ is modeled marginally using an Ordinal Autoregressive Logit Model of order $p$. Crucially, the marginal model for series $k$ includes lagged values from all $K$ series in the system to capture cross-dependence in the mean structure.
  $$\text{logit}(\gamma_{sjt}^{(k)}) = \alpha_{k0j} + \sum_{m=1}^K \sum_{g=1}^p \alpha'_{kmg} \tilde{Z}_{m, t-g}$$
  where $\gamma_{sjt}^{(k)}$ is the cumulative probability $P(Z_{kt} \leq s_j | \mathcal{F}_{t-1})$.
• Copula Linkage: To capture the contemporaneous (lag-0) cross-correlation, the authors use copulas. For any pair of series $(Z_{rt}, Z_{st})$, the joint cumulative distribution is defined as:
  $$F_{Z_{rt}, Z_{st}}(z_{rt}, z_{st} | \mathcal{F}_{t-1}) = C(F_{Z_{rt}}(z_{rt}), F_{Z_{st}}(z_{st}) | \mathcal{F}_{t-1}, \phi)$$
  where $C$ is a bivariate copula and $\phi$ is the dependence parameter.

B. Estimation via Pairwise Likelihood

Instead of maximizing the full $K$-dimensional likelihood, the authors maximize the Conditional Pairwise Log-Likelihood:
$$\ell_{prw}(\theta; z) = \sum_{t=2}^T \sum_{r=1}^{K-1} \sum_{s=r+1}^K \log f_{Z_{rt}, Z_{st}}(z_{rt}, z_{st}; \theta_{rs} | \mathcal{F}_{t-1})$$
This reduces the problem to estimating $K(K-1)/2$ bivariate models.

C. Two-Stage Estimation Procedure

To avoid the difficulty of optimizing the high-dimensional composite likelihood directly, the authors adopt a two-stage approach (following Fieuws and Verbeke, 2006):

1. Independent Maximization: Each of the $K(K-1)/2$ bivariate models is fitted independently to obtain individual parameter estimates $\tilde{\theta}_{rs}$.
2. Weighted Aggregation: The individual estimates are combined to form a unique set of parameters for the full system.
   • Instead of a simple mean, the authors propose a Weighted Mean based on the Hessian matrix of the log-likelihoods.
   • Using a second-order Taylor expansion, the weights are derived from the negative Hessian matrix $H(\tilde{\theta}_{rs})$.
   • The final estimator $\tilde{\theta}_{wm}$ is calculated as:
     $$\tilde{\theta}_{wm} = \left( \sum_{r,s} H(\tilde{\theta}_{rs}^f) \right)^{-1} \left( \sum_{r,s} H(\tilde{\theta}_{rs}^f) \tilde{\theta}_{rs}^f \right)$$
     where $\tilde{\theta}_{rs}^f$ is the vector padded with zeros to match the dimension of the full parameter vector $\theta$.

D. Forecasting

Forecasting is challenging because the full joint distribution is not explicitly available. The authors propose two methods to synthesize forecasts from the bivariate models:

• Method A (Mode Aggregation): Generate forecasts from each bivariate model individually, then select the final state as the mode (most frequent state) of these individual forecasts.
• Method B (Probability Aggregation): Average the conditional state probabilities ($\pi_{jt}$) from all bivariate models to create a combined probability distribution, then select the state with the highest probability (or the median) from this aggregated distribution.

3. Key Contributions

1. Computational Efficiency: The pairwise likelihood approach transforms an intractable $O(2^K)$ problem into a manageable set of $O(K^2)$ bivariate problems, making high-dimensional ordinal time series modeling feasible.
2. Flexible Dependence Modeling: The framework allows for different copula families to be selected for different pairs of time series, accommodating heterogeneous dependence structures (e.g., positive vs. negative correlations).
3. Improved Estimation Strategy: The introduction of a Hessian-weighted mean for synthesizing estimates provides a statistically rigorous way to combine independent bivariate results, outperforming simple averaging in simulation studies.
4. Forecasting Framework: The paper addresses the specific ambiguity of forecasting in composite likelihood settings by proposing and comparing two distinct aggregation strategies.

4. Results

Simulation Study

• Setup: A trivariate case ($K=3$) with 3 states was simulated using a Gumbel copula ($\phi=2$) and sample sizes $T = 100, 500, 1000$.
• Findings:
  • Both the simple mean and weighted mean estimators were unbiased and consistent (variance decreased as $T$ increased).
  • The weighted mean consistently outperformed the simple mean, particularly for the dependence parameter $\phi$, showing lower Mean Squared Error (MSE).
  • As sample size increased, MSE approached zero for both methods.

Real Data Application

• Dataset: Quarterly unemployment rates for six EU countries (Portugal, Slovakia, Finland, Greece, Spain, Italy) from 1998–2023.
• Preprocessing: Continuous rates were discretized into 4 ordinal states based on quantiles.
• Model: Ordinal autoregressive logit models with lag-1 effects from all 6 countries, linked by Frank copulas (chosen for their ability to handle both positive and negative correlations).
• Results:
  • The model successfully captured strong autocorrelation and cross-country correlations (e.g., strong links between Spain-Greece and Spain-Italy).
  • Forecasting (2021–2023): Both Method A and Method B provided similar results.
    • High Accuracy: Portugal, Slovakia, and Finland were forecasted with high probability of the correct state.
    • Lower Accuracy: Greece, Spain, and Italy showed poorer performance. The authors attribute this to the data exhibiting long "runs" of the same state (persistence), which the standard autoregressive logit model may not capture as well as specialized models (e.g., DAR models).

5. Significance and Conclusion

This paper provides a robust and scalable framework for analyzing complex multivariate ordinal data, which is prevalent in economics, health, and environmental sciences. By decoupling the estimation of marginal dynamics from the dependence structure via pairwise likelihood, the authors overcome the "curse of dimensionality" associated with full multivariate copulas.

Key Takeaways:

• The pairwise likelihood approach is a viable alternative to full likelihood for high-dimensional ordinal time series.
• Weighted averaging of estimates is superior to simple averaging, leveraging the curvature of the likelihood surface (Hessian) to improve precision.
• The methodology is flexible: it allows for different marginal models and different copula families for different pairs, though the authors note that model selection becomes a non-trivial task in such a flexible setting.
• Future work could explore incorporating covariates, higher-order lags, and specialized models for long-run persistence in ordinal data.