多変量確率的ボラティリティモデルで相関係数の時変性をとらえる

はじめに

この記事はマケデコ Advent Calendar 2025の9日目の記事です。

資産間の価格変動の相関係数を求める方法としては、過去一定期間のリターンからローリング相関を計算するものが広く用いられています。しかし、この方法には、ローリングのウィンドウの期間の長さによって値が変わったり、市場の急変に追随するのが遅かったりするという課題が存在します。

これを解決する方法として、ボラティリティと相関を同時にモデリングする、多変量の確率的ボラティリティモデル（Multivariate Stochastic Volatility (MSV) モデル）を紹介します。

MSVモデルをPythonとStanで実装し、2008年～2025年のTOPIXと東証REIT指数の相関係数を推定してみました。相関係数はおおむね0.4程度ですが、2016年～2020年や2025年は0.2程度まで低下していました。

金融工学の論文実装でアドカレを書いている人があんまりいないのでこういうのもありでしょう。マケデコのアドカレは気付けば3年連続で書いています。よかったら過去記事も読んでみてください。

毎年ボラティリティネタで書いてますね。これはわたしの趣味です。

ちなみにこの記事は、2023年の記事のStochastic Volatilityモデルを単一資産から複数資産に拡張したものですので、この過去の記事も参考になると思います。

相関係数のモデリングの重要性

資産間の相関係数はアセットアロケーションなどに使われるため、個々の資産のボラティリティとならんで重要なのはいうまでもありません。ではどうやって相関係数を求めればよいのでしょうか？

よく知られているのは、各資産の過去n日のリターンどうしの相関係数を求め、1日ずつスライドさせていく、いわゆるローリング相関係数です。

2008/5〜2025/12のTOPIXと東証REIT指数（両方配当なし）のデータを用いて、それぞれの対数リターンの過去n営業日ローリング相関係数をn=250でプロットしてみると以下のようになります。これは過去1年のローリングに相当します。なお、TOPIXと東証REIT指数のデータはJ-Quantsから取りました。

これは広く用いられている方法ですが、問題点があります。

特に2点目と3点目ですが、ローリング相関は、ウィンドウ期間中の相関係数が一定であることを暗黙の前提とします。その場合、ウィンドウ期間中の標本相関係数は真の相関の推定値として解釈できます。

しかし、市場は急に変動することがあるためその仮定は現実的ではありません。極端なリターンがあると相関係数が急変動する一方、相関係数の構造変化には追随が遅れてしまいます。極端なリターンの例としては、2011年3月の東日本大震災のときの市場の急変動により、プロットのとおり相関係数が急に上がっていることが挙げられます。日次リターンではなく月次リターンなどより長い期間のリターンから求めることで極端なリターンへの対応ができますが、構造変化への追随はより遅れます。

これらの問題に対する一つの解決策は、それぞれの資産の日次のリターンから、各資産のボラティリティとその相関係数をモデリングすることです。

様々な方法が提案されていますが1、その一つがMSVモデルです。これは、観測されるリターンから、潜在変数である各資産のボラティリティと資産間の相関係数を同時にモデリングするものです。

複雑な状態空間モデルですが、Stanなどで実装すればボラティリティや相関係数の推定値を得ることができます。またモデルの拡張も容易です。

MSVモデルをアセットアロケーションに適用した研究としては、Aguilar and West (2000) や、それを発展させてスパース性を持たせたZhou, et al. (2014) が有名です。特に前者では、MSVモデルを用いたアセットアロケーションは、1992年のイギリスのERM脱退時のポンド相場のような構造変化時によいパフォーマンスを示すことが報告されています。

モデル

いま、資産1と資産2の$t$日における日次終値をそれぞれ$S_{t,1}, S_{t,2}$とします。

このとき、$t$日におけるリターンを、対数リターンの100倍として、それぞれ$r_{t,1} = 100 (\log(S_{t,1}) - \log(S_{t-1,1})), r_{t,2} = 100 (\log(S_{t,2}) - \log(S_{t-1,2}))$とします。100倍するのは、パーセント表記にするということです。

以下のように定式化したMSVモデルを考えます。

$$ \begin{aligned} y_{t,i} &= \exp(h_{t,i}/2) \epsilon_{t,i}, \quad i = 1, 2 \\ (\epsilon_{t,1}, \epsilon_{t,2})’ &\sim N\left(\begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 1 & \rho_t \\ \rho_t & 1 \end{pmatrix}\right) \\ h_{t+1,i} &= \mu_i + \phi_i (h_{t,i} - \mu_i) + \eta_{t,i}, \quad \eta_{t,i} \sim N(0, \sigma_{\eta,i}^2) \\ g_{t+1} &= g_t + \zeta_t, \quad \zeta_t \sim N(0, \sigma_\zeta^2) \\ \rho_t &= \tanh(g_t) \end{aligned} $$

ただし、$y_{t,i} = r_{t, i}$とします。これは日次リターンが有意に正でも負でもないということです。多くの実証研究ではそれが示されていますが、仮にそうではない場合は、第1式に定数項か、その定数を時変にしてその状態方程式を加えます。

このとき、資産$i (i = 1, 2)$の$t$日におけるボラティリティは$\sigma_{t,i} = \exp(h_{t,i}/2)$、相関係数は$\rho_t$となります。

分散の対数がAR(1)過程に従い、相関係数はランダムウォークする$g_t$を$\tanh({g_{t}})$で相関係数が取るべき-1から1までの間に押し込めたものとして同時にモデリングしています。

なお、MSVモデルにはさまざまなバリエーションが存在します2。多変量の観測系列があり、各系列のボラティリティが確率的で時変する潜在変数であればMSVモデルといえます。

第1式と第2式はおおむね共通していますが、例えば相関係数をつくる$g$は上の式ではなく平均回帰性を持つようにAR(1)過程にするなど、細かい差異がいろいろ存在します。大森 (2019) に詳しく書かれています。

実装

上述のモデルをStanで書いて、Pythonからキックすることでパラメータを推定します。この実装のセクションはStanの知識が必要ですので、Stanに慣れていない方はこの章は読み飛ばしても大丈夫ですが、自分でやってみたい場合は参考にしてください。最近はStan以外にもベイズ推定のフレームワークがいろいろあるので好きなものを使ってください。

いま、以下のようにTOPIXと東証REIT指数の終値と対前日対数リターン（100倍したもの）をpolars.DataFrameで持っているとします。

shape: (4_303, 5)
┌────────────┬────────────┬───────────┬───────────┬───────────┐
│ Date       ┆ CloseTopix ┆ RetTopix  ┆ CloseReit ┆ RetReit   │
│ ---        ┆ ---        ┆ ---       ┆ ---       ┆ ---       │
│ date       ┆ f64        ┆ f64       ┆ f64       ┆ f64       │
╞════════════╪════════════╪═══════════╪═══════════╪═══════════╡
│ 2008-05-08 ┆ 1372.95    ┆ -1.469897 ┆ 1585.12   ┆ 0.608743  │
│ 2008-05-09 ┆ 1341.76    ┆ -2.297952 ┆ 1550.03   ┆ -2.238583 │
│ 2008-05-12 ┆ 1342.79    ┆ 0.076735  ┆ 1535.05   ┆ -0.971133 │
│ 2008-05-13 ┆ 1360.05    ┆ 1.277192  ┆ 1547.59   ┆ 0.813593  │
│ 2008-05-14 ┆ 1373.04    ┆ 0.95058   ┆ 1556.18   ┆ 0.553522  │
│ …          ┆ …          ┆ …         ┆ …         ┆ …         │
│ 2025-12-01 ┆ 3338.33    ┆ -1.194338 ┆ 1995.68   ┆ -1.399209 │
│ 2025-12-02 ┆ 3341.06    ┆ 0.081744  ┆ 1998.79   ┆ 0.155715  │
│ 2025-12-03 ┆ 3334.32    ┆ -0.201936 ┆ 1986.95   ┆ -0.59412  │
│ 2025-12-04 ┆ 3398.21    ┆ 1.898006  ┆ 1973.22   ┆ -0.693407 │
│ 2025-12-05 ┆ 3362.56    ┆ -1.054623 ┆ 1962.18   ┆ -0.561063 │
└────────────┴────────────┴───────────┴───────────┴───────────┘

上のMSVモデルを以下の通り書き下し3、“model.stan”として保存します。事前分布はKim, Shephard and Chib (1998), 大森 (2019) にならっています。

Stanコード（クリックすると折りたたみが開きます）
// Stanの収束をよくするためのテクニックを以下の通り入れている
// 1. phiのlower, upperを-0.999から0.999に縛っている
// -1 - 1だと両端でサンプリングが不安定になって収束が悪いため
//（sigma系のパラメータのacfやESSが微妙になる）
// 2. h_raw, g_rawの非中心パラメータ化（「再パラメータ化」の一種）
// sigma_etaが小さいとき、hとsigma_etaの事後分布が強く相関するので、h_rawとsigma_etaを分離する

data {
  int<lower=0> n; // 時点数
  int<lower=1> p; // 次元（p=2 であることを前提としている）
  matrix[p, n] y; // リターン（2×n）
}

parameters {
  matrix[p, n] h_raw;
  vector[n] g_raw;

  vector[p] mu;
  vector<lower=0.0005, upper=0.9995>[p] phi_raw;
  vector<lower=0>[p] sigma_eta_sq;
  real<lower=0> sigma_zeta;
}

transformed parameters {
  vector<lower=-1, upper=1>[p] phi = 2 * phi_raw - 1;
  vector<lower=0>[p] sigma_eta = sqrt(sigma_eta_sq);

  matrix[p, n] h;
  vector[n] g;
  vector<lower=-1, upper=1>[n] rho;

  // 非中心パラメータ化
  for (i in 1:p) {
    h[i, 1] = mu[i] + (sigma_eta[i] / sqrt(1 - phi[i]^2)) * h_raw[i, 1];
  }
  g[1] = 10 * g_raw[1];
  for (t in 2:n) {
    for (i in 1:p) {
      h[i, t] = mu[i] + phi[i] * (h[i, t-1] - mu[i]) + sigma_eta[i] * h_raw[i, t];
    }
    g[t] = g[t-1] + sigma_zeta * g_raw[t];
  }

  for (t in 1:n) {
    rho[t] = tanh(g[t]);
  }
}

model {
  // 論文のとおり事前分布を設定する
  // sigma_zetaは論文とはモデルが違うのでflat priorにしている
  mu ~ normal(0, 1);
  phi_raw ~ beta(20, 1.5);
  sigma_eta_sq ~ inv_gamma(2.5, 0.025);

  to_vector(h_raw) ~ std_normal();
  g_raw ~ std_normal();

  // 観測方程式
  for (t in 1:n) {
    real sigma1 = exp(h[1, t] / 2.0);
    real sigma2 = exp(h[2, t] / 2.0);

    matrix[p, p] Sigma;
    Sigma[1, 1] = sigma1^2;
    Sigma[2, 2] = sigma2^2;
    Sigma[1, 2] = rho[t] * sigma1 * sigma2;
    Sigma[2, 1] = Sigma[1, 2];

    y[, t] ~ multi_normal(rep_vector(0.0, p), Sigma);
  }
}

generated quantities {
  matrix[p, n] volatility;
  vector[n] log_lik;

  for (t in 1:n) {
    real sigma1 = exp(h[1, t] / 2.0);
    real sigma2 = exp(h[2, t] / 2.0);

    volatility[1, t] = sigma1;
    volatility[2, t] = sigma2;

    matrix[p, p] Sigma;
    Sigma[1, 1] = sigma1^2;
    Sigma[2, 2] = sigma2^2;
    Sigma[1, 2] = rho[t] * sigma1 * sigma2;
    Sigma[2, 1] = Sigma[1, 2];

    log_lik[t] = multi_normal_lpdf(y[, t] | rep_vector(0.0, p), Sigma);
  }
}

MSVモデルは素朴に実装するとMCMCの収束が悪くなりがちなので、非中心パラメータ化のようなテクニックを使っています（Stanコードのコメントを参照）。

iter_warmup=1000, iter_sampling=1000, thin=1, chain=4で、私の環境では15分くらいで推定できました。なお、使用した環境はPython=3.13.7, CmdStan=2.36, cmdstanpy=1.3.0です。

MCMCサンプリングが収束したこと（パラメータ推定に問題がないこと）を確認しましたが、割愛します4。

結果

結果の相関係数です。線の上下の帯は95%ベイズ信用区間です。

おおむね0.4程度を取っていますが、2016-2020年や2025年は0.2程度まで低下していますね。

2023年時点での過去10年の月次相関係数を計算すると、国内株と国内REITは0.4程度とのことです。（参考: 資産配分を再考、最適な投信の組み合わせは？ | 東証マネ部！）

これと整合的な結果が得られましたが、時期によってはそれより小さいことも大きいこともあることが分かり、市場の変化に追随できることを示せました。また、各時点での相関係数の信用区間を示すこともでき、不確実性の程度が分かるのもよい点です。

基本的に国内株も国内REITも金利と逆相関の動きをしますが、株は景気や業績の影響も強く受けるため、2025年のように金利上昇局面でも買われることがありますね。このような場合に相関関係が崩れるのだと思われます。

補足

冒頭のローリング相関のプロットと異なり、滑らかな曲線を描いています。これは、各時点$t$の相関係数を全期間のリターンデータを用いて推定していますが、ローリング相関のプロットは$t$までのデータ（さらに直近n期間のみ）から求めていることの違いによるものです。

状態空間モデルの文脈で言うと、前者は「平滑化」と呼ばれる推定量に近いものであり、後者は「フィルタ化」と呼ばれる推定量に近いです（ただし、ローリング相関は状態空間モデルではないので、あくまで$t$の相関係数を$t$までのデータから求めているという意味でフィルタ化に近いということです）。

MCMCで求めたパラメータの推定値は平滑化推定量に近いものです。将来の時点のデータを使って過去の時点のパラメータを求めているためバックテストには使えませんが、後から振り返ってデータを解釈する目的では適していますね。

バックテストに使うためには、MCMCで毎日回すか、あるいは粒子フィルタのようなフィルタ系の手法で推定してフィルタ化推定量を求めるという解決策があります。後者ですが、推定にかかる時間もフィルタ系の方が早いので、毎日や毎週のオンライン予測には適しています。

なお、資産価格の系列数が多くなると共分散行列の要素数が多くなるためパラメータの推定が困難になります。その場合は少ない因子に分解して推定するなどの方法を取ります。これも大森 (2019) を参照してください。

おわりに

金融工学の論文は実装すればすぐ収益に結びつくようなものではありませんが、論文の読解と実装を地道に積み重ねていくことで知見が得られると思っています。やはり、巨人の肩の上に立つことには意味があるんだと思います。

参考文献

  1. 他に、分足データを用いて相関係数の推定値を求めるRealized Correlationという手法もあります。2024年のアドカレで書いたRealized Volatilityの記事の多変量版の拡張になります。 ↩︎

  2. ここで書いたモデルは、大森 (2019) の「動学的均一相関MSVモデル」を2変量に限定して、$g$をAR(1)ではなくランダムウォークとし、かつ相関係数が正だけではなく負の値も取れるように拡張したものにおおむね相当します。 ↩︎

  3. なお、TOPIXのリターンは平均0.020, 標準偏差1.345, 東証REIT指数のリターンは平均0.005, 標準偏差1.336であり、有意に正でも負でもないため、前掲のモデルの第1式に定数項を入れないことは妥当な定式化と言えます。 ↩︎

  4. 内容はこちらを参考にしてください: [R] [stan] bayesplot を使ったモンテカルロ法の実践ガイド - ill-identified diary, [R][Stan]マルコフ連鎖モンテカルロ法の実践ガイド2: ランクプロット他 - ill-identified diary ↩︎