第97回 分散減少法 ― 対照変量、制御変量、準乱数

位置: 第10部 第97回 / 前提: [96] / 数学編: [M32] / この回でできるようになること: 対照変量法や制御変量法を使って、同じ計算量でモンテカルロ法の誤差を小さくすることができる

0. 物語 ― 相関の低い制御変量は効かない

あるチームがモンテカルロ法でオプション価格を計算していました。標本を増やせば精度は上がりますが、計算に何時間もかかるので、もっと賢い方法を試そうとしたのです。教科書に「期待値がわかっている変数で補正すれば、分散が減る」と書いてあります。制御変量法と呼ばれる手法です。

チームはこう考えました。「何でもいいから、期待値がわかる変数を持ってくればいい」。そこで、オプションのペイオフとは何の関係もない乱数をひとつ用意して、補正に使いました。期待値がゼロだとわかっているから、条件は満たしています。

結果は散々でした。10万本のシミュレーションを回しても、素朴に計算したときと標準誤差がまったく同じです。分散の減少率はゼロ。何も起きませんでした。

失敗の原因は「関係のない変数を持ってきた」ことです。制御変量法が効くかどうかは、補正に使う変数と、本来計算したい量との相関で決まります。相関がほぼゼロなら、いくら期待値がわかっていても分散は減りません。

一方、同じ問題で、ペイオフと強くつながっている変数を使ったチームがいました。こちらは相関が0.89もあり、分散を約79パーセント減らすことに成功しています。同じ計算量で標準誤差がおよそ半分以下になったのです。

分散を減らすには「正しい相手」を選ぶ必要があります。この回では、どんな変数を選べばどれだけ分散が減るのか、その仕組みを数式で明らかにします。

1. 動機 ― 標本を増やさずに精度を上げたい

この回を読み終えると、対照変量法や制御変量法を使って、同じ計算量でモンテカルロ法の誤差を小さくすることができます。

前回（第96回）で素朴なモンテカルロ法を動かし、標準誤差が $1/\sqrt{N}$ でしか縮まないことを確認しました。精度を10倍にしたければ標本を100倍に増やす必要があり、計算コストが急激に膨らみます。

しかし「同じ標本数のまま、推定量の分散だけ下げる」方法があります。鍵は乱数の間にある相関を意図的に作り出すことです。今回は3つの手法を取り上げます。

• 対照変量法（antithetic variates）: 乱数とその符号反転をペアで使う。
• 制御変量法（control variates）: 期待値が既知の変数で推定量を補正する。
• 準乱数（quasi-random numbers）: 空間を均一に埋める点列で擬似乱数を置き換える。

いずれも「同じ計算量で、より小さい標準誤差を得る」ことが目的です。

この記事で証明すること

• THM-VR-1 正規乱数とその符号反転をペアにすると、ペイオフが単調な場合に推定量の分散が小さくなります。
• THM-VR-2 期待値が既知の変数で推定量を補正するとき、分散を最小にする係数は共分散を分散で割った値であり、そのとき分散は相関係数の二乗の割合だけ減ります。
• THM-VR-3 低食い違い列を使うと、次元が低い問題では誤差の縮み方が標本数の逆数に近づきます。

2. 設定と仮定 ― 同じ計算量で比べる

第96回と同じ設定を使います。

• 原資産価格 $S_0=100$、行使価格 $K=105$、無リスク金利 $r=0.05$、ボラティリティ $\sigma =0.2$、満期 $T=1$。
• ヨーロピアンコールのペイオフは $\max (S_T-K,0)$。
• BS 公式による理論価格は $8.0214$（小数4桁に丸め）。
• 計算量はペイオフの評価回数で測ります。対照変量法は $N/2$ 個の乱数から $N/2$ 組のペア（ペイオフの評価は $N$ 回）を作るので、比較は評価回数を揃えて行います。
• 乱数のシードは $42$ に固定します。

制御変量には原資産の満期価格 $S_T$ を使います。$S_T$ はリスク中立測度のもとで

$${\mathbb{E}}^{\mathbb{Q}}[S_T]=S_0{e}^{rT}$$

であり、この期待値は解析的に求まります。

3. 主張 ― 相関を使うと分散が減る

定理 THM-VR-1（対照変量法）

標準正規乱数 $Z$ から満期株価を

$$S_T=S_0\exp ((r-\frac{1}{2}{\sigma }^2)T+\sigma \sqrt{T}Z)$$

で生成します。$Z$ を $-Z$ に置き換えた株価を $S_T^{\prime }$ とし、対応するペイオフをそれぞれ $f(Z)$, $f(-Z)$ と書きます。

定理 THM-VR-1（対照変量法） $f$ が $Z$ の単調関数で $\mathrm{Var}[f(Z)]<\infty$ であるとき、推定量

$${\hat{\mu }}_{\mathrm{AV}}=\frac{1}{2}(f(Z)+f(-Z))$$

の分散は、独立な2つの標本 $Z_1,Z_2$ による素朴な推定量 $\frac{1}{2}(f(Z_1)+f(Z_2))$ の分散以下である。すなわち

$$\mathrm{Var}\left[{\hat{\mu }}_{\mathrm{AV}}\right]\le \frac{1}{2}\mathrm{Var}[f(Z)].$$

定理 THM-VR-2（制御変量法の最適係数）

推定したい量を $X=e^{-rT}f(Z)$、期待値が既知の変数を $Y$（ここでは $S_T$）とし、補正した推定量を

$$X^{(c)}=X-c(Y-\mathbb{E}[Y])$$

と定義します。

定理 THM-VR-2（制御変量法の最適係数） $\mathrm{Var}[X^{(c)}]$ を最小にする $c$ は

$$c^{\ast }=\frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}(Y)}$$

であり、そのときの分散は

$$\mathrm{Var}[X^{(c^{\ast })}]=(1-{\rho }_{XY}^2)\mathrm{Var}[X]$$

である。ここで ${\rho }_{XY}$ は $X$ と $Y$ の相関係数である。

定理 THM-VR-3（準乱数）

定理 THM-VR-3（準乱数の誤差評価） $d$ 次元の低食い違い列（Sobol 列など）を $N$ 点使った積分近似の誤差は

$$\left|\frac{1}{N}\sum _{i=1}^{N}g({\mathbf{x}}_i)-{\int }_{[0,1{]}^d}g(\mathbf{x})d\mathbf{x}\right|\le V(g)D_N^{\ast }$$

で抑えられる（Koksma–Hlawka の不等式）。ここで $V(g)$ は $g$ の Hardy–Krause 変動、$D_N^{\ast }$ は点列の星食い違い度である。低食い違い列では $D_N^{\ast }=O((\log N{)}^d/N)$ であり、次元 $d$ が小さいとき誤差の縮み方は $O(1/N)$ に近づく。

証明を読む（読まなくても記事は完結します）

4. 証明 ― 最適係数を分散の最小化から出す

この節は読まなくても記事は完結します。

この証明を追いたい人は、数学編 第32回「乱数生成」へ。

THM-VR-1 の証明

${\hat{\mu }}_{\mathrm{AV}}=\frac{1}{2}(f(Z)+f(-Z))$ の分散を展開します。$Z$ と $-Z$ は同じ分布に従うので $\mathrm{Var}[f(Z)]=\mathrm{Var}[f(-Z)]$ です。

$$\begin{array}{rl}& \mathrm{Var}[{\hat{\mu }}_{\mathrm{AV}}]\\ & =\frac{1}{4}(\mathrm{Var}[f(Z)]+\mathrm{Var}[f(-Z)]\\ & +2\mathrm{Cov}[f(Z),f(-Z)])\\ & =\frac{1}{2}\mathrm{Var}[f(Z)]+\frac{1}{2}\mathrm{Cov}[f(Z),f(-Z)].\end{array}\text{\hspace{1em}(4.1)}$$

第1行から第2行では、$-Z$ も標準正規分布に従うので分散が等しいことを使いました。

素朴な推定量と比較します。比較はペイオフの評価回数を揃えて行う必要があります。${\hat{\mu }}_{\mathrm{AV}}$ は1組でペイオフを2回評価するので、比べる相手は独立な2つの標本の平均 $\frac{1}{2}(f(Z_1)+f(Z_2))$ です。独立性から交差項が消えるので、その分散は $\frac{1}{2}\mathrm{Var}[f(Z)]$ です。差を取ると、(4.1) より

$$\begin{array}{rl}& \mathrm{Var}[{\hat{\mu }}_{\mathrm{AV}}]-\frac{1}{2}\mathrm{Var}[f(Z)]\\ & =\frac{1}{2}\mathrm{Cov}[f(Z),f(-Z)].\end{array}$$

したがって、対照変量法が有利（以下）になる条件は

$$\mathrm{Cov}[f(Z),f(-Z)]\le 0\text{\hspace{1em}(4.2)}$$

です。$f$ を単調増加として (4.2) を示します。$g(z)=f(-z)$ とおくと $g$ は単調減少です。$Z^{\prime }$ を $Z$ と独立で同じ分布に従う確率変数とすると、共分散は

$$\begin{array}{rl}& \mathrm{Cov}[f(Z),g(Z)]\\ & =\frac{1}{2}\mathbb{E}[(f(Z)-f(Z^{\prime }))(g(Z)-g(Z^{\prime }))]\end{array}$$

と書けます（右辺を展開すると $\mathbb{E}[f(Z)g(Z)]-\mathbb{E}[f(Z)]\mathbb{E}[g(Z)]$ が2回現れ、$\frac{1}{2}$ で1回分になります）。$Z>Z^{\prime }$ なら $f(Z)-f(Z^{\prime })\ge 0$ かつ $g(Z)-g(Z^{\prime })\le 0$、$Z<Z^{\prime }$ なら符号がそれぞれ逆になるので、期待値の中の積はいつも $0$ 以下です。よって共分散は $0$ 以下で、(4.2) が成り立ちます。$f$ が単調減少なら $f$ と $g$ の役割が入れ替わるだけで、同じ結論になります。

コールオプションのペイオフ $\max (S_T-K,0)$ は $S_T$ の単調増加関数であり、$S_T$ は $Z$ の単調増加関数です。よって合成 $f(Z)$ は $Z$ の単調増加関数であり、(4.2) が成り立ちます。 ∎

THM-VR-2 の証明

$X^{(c)}=X-c(Y-\mathbb{E}[Y])$ の分散を $c$ の関数として書きます。$Y-\mathbb{E}[Y]$ の期待値はゼロなので、$\mathbb{E}[X^{(c)}]=\mathbb{E}[X]$ であり、補正しても推定量の不偏性は保たれます。

分散を展開します。

$$\begin{array}{rl}& \mathrm{Var}[X^{(c)}]\\ & =\mathrm{Var}[X]-2c\mathrm{Cov}(X,Y)\\ & +c^2\mathrm{Var}(Y).\end{array}\text{\hspace{1em}(4.3)}$$

第1行から第2行では、分散の二次形式 $\mathrm{Var}[A-B]=\mathrm{Var}[A]-2\mathrm{Cov}(A,B)+\mathrm{Var}[B]$ を用いました。

(4.3) は $c$ の二次関数で、$\mathrm{Var}(Y)>0$ より上に凸ではなく下に凸です。最小値を与える $c$ を求めるため、$c$ で微分してゼロとおきます。

$$\frac{d}{dc}\mathrm{Var}[X^{(c)}]=-2\mathrm{Cov}(X,Y)+2c\mathrm{Var}(Y)=0.$$

これを解くと

$$c^{\ast }=\frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}(Y)}.\text{\hspace{1em}(4.4)}$$

$c=c^{\ast }$ を (4.3) に代入して最小分散を求めます。

$$\begin{array}{rl}& \mathrm{Var}[X^{(c^{\ast })}]\\ & =\mathrm{Var}[X]-\frac{2\mathrm{Cov}(X,Y{)}^2}{\mathrm{Var}(Y)}\\ & +\frac{\mathrm{Cov}(X,Y{)}^2}{\mathrm{Var}(Y)}\\ & =\mathrm{Var}[X]-\frac{\mathrm{Cov}(X,Y{)}^2}{\mathrm{Var}(Y)}.\end{array}\text{\hspace{1em}(4.5)}$$

第1行から第2行では、第2項と第3項をまとめました。第2項は $-2\mathrm{Cov}(X,Y{)}^2/\mathrm{Var}(Y)$、第3項は $+\mathrm{Cov}(X,Y{)}^2/\mathrm{Var}(Y)$ なので、合わせて $-\mathrm{Cov}(X,Y{)}^2/\mathrm{Var}(Y)$ です。

相関係数 ${\rho }_{XY}=\mathrm{Cov}(X,Y)/\sqrt{\mathrm{Var}(X)\mathrm{Var}(Y)}$ を使って書き直すと

$$\begin{array}{rl}\mathrm{Var}[X^{(c^{\ast })}]& =\mathrm{Var}[X]-{\rho }_{XY}^2\mathrm{Var}[X]\\ & =(1-{\rho }_{XY}^2)\mathrm{Var}[X].\end{array}\text{\hspace{1em}(4.6)}$$(4.6)

THM-VR-3 の証明の要点

Koksma–Hlawka の不等式の完全な証明は解析的整数論の手法を使うため、ここでは要点のみ述べます。

1. 食い違い度 $D_N^{\ast }$ は、点列 $\{{\mathbf{x}}_1,\dots ,{\mathbf{x}}_N\}$ が $[0,1{]}^d$ をどれだけ均一に埋めているかを測る量です。任意の軸に沿った矩形 $[0,a_1)\times \cdots \times [0,a_d)$ に入る点の割合と、その矩形の体積との差の上限として定義されます。

2. 被積分関数 $g$ の Hardy–Krause 変動 $V(g)$ は、$g$ の滑らかさを測る量です。$g$ が有界変動であれば有限です。

3. 不等式は「点の偏りが小さく、関数がなめらかなら、積分の近似誤差は小さい」と読めます。擬似乱数では $D_N^{\ast }=O(1/\sqrt{N})$（大数の法則の速度）ですが、Sobol 列などの低食い違い列では $D_N^{\ast }=O((\log N{)}^d/N)$ が成り立ちます。

4. 次元 $d$ が小さいとき $(\log N{)}^d/N\ll 1/\sqrt{N}$ であり、準乱数が擬似乱数より速く収束します。しかし $d$ が大きくなると $(\log N{)}^d$ の項が支配的になり、理論上の優位は薄れます。

以上が要点です。 ∎

5. 反例 ― 相関が低いと効かない

CEX-VR-1: 対照変量はいつも効くわけではない

対照変量法が有効なのは、ペイオフが $Z$ の単調関数であるときでした。では、ペイオフが対称な関数のときはどうなるでしょうか。

$|S_T-{\mathbb{E}}^{\mathbb{Q}}[S_T]|$ というペイオフを考えます。これは $S_T$ が期待値から上に外れても下に外れても同じ値を取ります。$Z$ を $-Z$ に置き換えると $S_T$ は期待値の反対側に動きますが、絶対値を取るので $f(-Z)$ は $f(Z)$ と近い値になります。共分散が正になり、(4.1) の第2項が分散を増やす方向に働きます。(4.2) が破れるので、評価回数を揃えた比較では素朴法より悪くなります。

数値で確認します（$N=100,000$）。

手法	標準誤差
素朴法	$0.0413$
対照変量法	$0.0568$

対照変量法の標準誤差は素朴法の $1.38$ 倍に増えています。ペイオフが対称な場合、対照変量法は逆効果です。

CEX-VR-2: 相関が低い変数では制御変量が効かない

制御変量法の分散減少率は ${\rho }_{XY}^2$ に比例します（(4.6)）。では、制御変量にペイオフと無相関な変数を使うとどうなるでしょうか。

独立な標準正規乱数を制御変量に使った場合の結果です。

	値
相関係数 $\rho$	$0.0017$
${\rho }^2$	$0.000003$
素朴法の標準誤差	$0.0419$
制御変量法の標準誤差	$0.0419$
分散減少率	$0.0\%$

相関がほぼゼロなので、分散はまったく減りません。制御変量法が威力を発揮するには、ペイオフと強く相関する変数を見つけることが本質的です。

CEX-VR-3: 準乱数は高次元で優位が薄れる

THM-VR-3 で見たように、Koksma–Hlawka の不等式の誤差上界には $(\log N{)}^d$ の項が含まれます。次元 $d$ を上げると、この項が大きくなり、準乱数の優位が縮まります。

次元を変えて、ヨーロピアンコールを $d$ 段階の経路で評価した結果です（$N=8,192$）。準乱数の標準誤差は、互いに独立にランダム化（スクランブル）した $16$ 組の Sobol 点列（各 $512$ 点）の推定値のばらつきから求めています。1本の Sobol 点列を区切った部分どうしは独立ではないので、区切った平均のばらつきからは正しい標準誤差が出ないためです。

次元 $d$	準乱数の標準誤差	擬似乱数の標準誤差	比（擬似/準）
1	$0.0053$	$0.1420$	$26.99$
5	$0.0130$	$0.1468$	$11.28$
10	$0.0386$	$0.1448$	$3.75$
50	$0.0364$	$0.1448$	$3.98$
100	$0.0354$	$0.1482$	$4.19$

1次元では準乱数が約27倍有利ですが、10次元以上では4倍前後まで縮まります（10次元から100次元の間の違いは、16組から見積もった標準誤差のばらつきの範囲です）。高次元の問題では準乱数だけに頼るのは危険です。

6. 数値例 ― 4つの手法の誤差を並べる

第96回と同じパラメータ（$S_0=100$, $K=105$, $r=0.05$, $\sigma =0.2$, $T=1$）で、BS 理論価格 $8.0214$ に対する各手法の結果です。標本数は $N=100,000$（準乱数は $N=131,072=2^{17}$ で、独立にスクランブルした $16$ 組の Sobol 点列（各 $8,192$ 点）の推定値のばらつきから標準誤差を求めます）。

手法	推定値	標準誤差
素朴モンテカルロ	$8.0063$	$0.0419$
対照変量法	$8.0694$	$0.0334$
制御変量法	$8.0463$	$0.0191$
準乱数（Sobol）	$8.0217$	$0.0003$

同じ精度を得るのに必要な標本数の比を、素朴法を基準に計算します。標準誤差が $\sigma /\sqrt{N}$ に比例するので、同じ標準誤差を得るのに必要な標本数の比は分散の比の逆数です。

手法	効率比（素朴法の何倍速いか）
素朴モンテカルロ	$1.0$
対照変量法	$1.6$
制御変量法	$4.8$
準乱数	$11,015$

制御変量法は素朴法の約5倍、準乱数は約1万倍の効率です。

制御変量法の最適係数と分散減少率の一致も確認します。

	値
最適係数 $c^{\ast }$	$0.5541$
相関係数 $\rho$	$0.8906$
${\rho }^2$	$0.7933$
理論上の分散比 $1-{\rho }^2$	$0.2067$
実測の分散比	$0.2067$

(4.6) の理論値と実測値が一致しています。${\rho }^2=0.79$ なので、制御変量法は分散を約 $79\%$ 削減しています。

これらの結果はすべて compute.py により再現可能です。

7. 次回へ ― 偏微分方程式を格子で解く

今回はモンテカルロ法の精度を上げる方法を見ました。対照変量法は単調なペイオフに有効であり、制御変量法は相関の高い変数が見つかれば強力です。準乱数は低次元で劇的な効果を発揮しますが、高次元では優位が薄れます。

次回（第98回）では、偏微分方程式を格子で直接解く有限差分法を扱います。モンテカルロ法とは全く異なるアプローチで、BS 方程式の数値解を求めます。

よくある質問

どの方法が一番良いのですか

問題によります。期待値が分かる近い変数があれば制御変量が強力です。

準乱数とは何ですか

空間を均一に埋めるように作った点列です。乱数ではありませんが、積分の近似に使えます。

組み合わせて使えますか

使えます。実務では複数を組み合わせるのが普通です。

参考文献

1. Glasserman P.、Monte Carlo Methods in Financial Engineering、Springer、2004年
2. Niederreiter H.、Random Number Generation and Quasi-Monte Carlo Methods、SIAM、1992年

---

数値例は compute.py によりシード固定で検算済み。本文中の小数は4桁で丸めている。