論文サマリー：対尤度を用いた多変量順序時系列のモデリング

著者: Anna Nalpantidi, Dimitris Karlis
所属: アテネ経済大学統計学部（ギリシャ）
日付: 2026 年 2 月（arXiv 投稿日）

1. 研究の背景と問題設定

順序データ（Ordinal Data）の時系列分析は、環境問題、健康研究、経済学など多岐にわたる分野で応用されています。しかし、**多変量順序時系列（Multiple Ordinal Time Series）**の同時分布をモデル化する手法は、以下の理由から限られており、発展が困難でした。

1. 分布の構築難易度: 順序変数の多変量分布を容易に構築・適用できる手法が存在しない。
2. 複雑な依存構造: 各時系列内の自己相関（Autocorrelation）、系列間の同時相関（Cross-correlation）、およびラグを伴う相互相関をすべて考慮する必要がある。
3. 計算コスト: 高次元の結合分布（特にコピュラを用いた場合）を直接評価すると、計算量が爆発的に増大し、実用的ではない。

既存の研究では、高次元コピュラを直接仮定する方法や、潜在連続変数を介したモデルが提案されていますが、次元が増えると確率質量関数の計算が困難になるという課題が残っていました。

2. 提案手法（Methodology）

本研究は、**対尤度（Pairwise Likelihood）**アプローチを採用し、高次元モデルの計算的課題を回避しつつ、多変量順序時系列を効率的にモデル化する枠組みを提案しています。

2.1 モデルの構造

• 周辺モデル（Marginal Models）:
  各時系列 $Z_{kt}$ は、順序回帰ロジットモデル（Ordinal Autoregressive Logit Model）に従うと仮定します。
  $$\text{logit}(\gamma^{(k)}_{sjt}) = \alpha_{k0j} + \sum_{m=1}^{K} \sum_{g=1}^{p} \alpha'_{kmg} \tilde{Z}_{m,t-g}$$
  ここで、$\gamma$ は累積確率、$\tilde{Z}$ は二値表現されたラグ変数です。重要な点は、各系列のモデルにシステム内のすべての系列のラグ値を含めることで、系列間の相互依存性を捕捉していることです。
• 依存構造の捕捉（コピュラ）:
  同時刻（ラグ 0）における系列間の依存構造を捉えるために、**コピュラ（Copula）**を使用します。
  • 従来のアプローチでは $K$ 変数のコピュラを直接指定して計算する必要がありましたが、本研究ではすべての可能なペア（2 変量）の周辺分布のみを指定します。
  • 各ペア $(Z_{rt}, Z_{st})$ に対して、コピュラ $C(\cdot, \cdot; \phi)$ を用いた結合累積分布関数を定義し、そこから同時確率質量関数を導出します。

2.2 推定手法：対尤度と重み付き平均

高次元の対尤度を直接最大化するのではなく、以下の 2 段階のプロセスで推定を行います（Fieuws and Verbeke, 2006 の枠組みを拡張）。

1. 個別最大化:
   $K(K-1)/2$ 個のすべてのペアについて、条件付き対尤度を個別に最大化し、パラメータ推定値 $\tilde{\theta}_{rs}$ を得ます。これにより、高次元コピュラの計算を回避し、異なるペアに対して異なるコピュラ族を選択する柔軟性も得られます。
2. 推定値の統合（重み付き平均）:
   個別に得られた推定値を統合して、最終的なパラメータ $\hat{\theta}$ を算出します。
   • 単純な平均ではなく、ヘッセ行列（Hessian Matrix）の逆行列に基づく重みを用いた重み付き平均を提案しています。
   • 各ペアの尤度関数を 2 次テイラー展開し、重み付き最小二乗基準を最大化することで、推定値の分散を最小化する効率的な統合を行います。
   • 式 (8) に示されるように、ヘッセ行列 $H(\tilde{\theta}_{rs})$ を重みとして使用します。

2.3 標準誤差の算出

対尤度は合成尤度（Composite Likelihood）の一種であるため、Fieuws and Verbeke (2006) の漸近理論に基づき、情報行列（J）とバリアンス行列（K）を用いて標準誤差を算出します。

3. シミュレーション研究の結果

3 変量ケース（状態数 3）において、サンプルサイズ $T=100, 500, 1000$ でシミュレーションを行い、以下の結果が得られました。

• 不偏性と一致性: 全てのケースで、推定値の中央値は真値に非常に近く、不偏推定量であることが確認されました。サンプルサイズが増加するにつれて分散が減少し、一致性が確認されました。
• 重み付き平均の有効性: パラメータ推定値全体では単純平均と重み付き平均の差は僅かでしたが、依存パラメータ（コピュラパラメータ $\phi$）において、重み付き平均が単純平均よりも明確に優れた性能（より小さな分散と MSE）を示しました。
• 計算効率: 高次元コピュラの直接計算を回避できるため、計算コストが大幅に削減されます。

4. 実データ分析：EU 各国の失業率

1998 年から 2023 年までの 6 カ国（ポルトガル、スロバキア、フィンランド、ギリシャ、スペイン、イタリア）の四半期別失業率データを対象に分析を行いました。

• データ変換: 失業率を四分位点に基づいて 4 つの順序カテゴリ（1: 低〜4: 高）に変換しました。
• モデル適合: 各系列のラグ 1 までの自己相関と相互相関を捕捉するため、1 階の順序回帰ロジットモデルと Frank コピュラ（正負の相関を両方扱える）を使用しました。
• 結果:
  • 各国の失業状態の自己相関は強く、特にスペインとギリシャ、スペインとイタリアの間で強い相関が確認されました。
  • 推定されたパラメータは統計的に有意であり、モデルがデータ構造を適切に捉えていることが示されました。

5. 予測（Forecasting）の課題とアプローチ

対尤度アプローチでは、完全な多変量分布が存在しないため、予測には特有の課題があります。各系列に対して複数の二変量モデルから予測が得られるため、これらを統合する必要があります。

本研究では 2 つの予測統合手法を提案・比較しました：

1. Method A（状態の多数決）: 各二変量モデルからランダムに状態をサンプリングし、そのモード（最頻値）を最終予測とする。
2. Method B（確率の平均）: 各二変量モデルから得られる状態ごとの条件付き確率を平均し、その平均確率分布からモードを予測する。

結果: 両手法とも類似した結果を示しましたが、安定した状態（同じ状態が長く続くデータ）を持つ国（ギリシャ、スペインなど）では予測精度が低下する傾向がありました。これは、順序データにおける「長期の連続性」を捉えるための特別なモデル（DAR モデル等）の必要性を示唆しています。

6. 結論と意義

主要な貢献

1. 計算効率化: 高次元コピュラの直接計算を回避し、対尤度（ペアワイズ）アプローチを用いることで、多変量順序時系列の推定を計算的に実行可能にしました。
2. 柔軟なモデリング: 各ペアに対して異なるコピュラ族を選択できる柔軟性を持ち、自己相関と相互相関を同時にモデル化できます。
3. 効率的な推定戦略: 個別推定値をヘッセ行列に基づく重み付き平均で統合する手法を提案し、特に依存パラメータの推定精度を向上させました。
4. 予測フレームワーク: 不完全な情報（二変量情報のみ）から多変量時系列を予測するための実用的な統合手法を提示しました。

意義と今後の展望

この手法は、経済指標、気象データ、医療データなど、順序尺度で観測される多変量時系列の分析において強力なツールとなります。
今後の課題としては、異なる系列に対して異なる周辺モデルを適用する一般化、異なるペアに対して異なるコピュラ族を許容するモデル選択問題、および高次ラグモデルへの拡張が挙げられます。また、予測精度向上のため、状態の持続性を考慮したモデルの組み合わせも検討の余地があります。