メインコンテンツまでスキップ

「math」タグの記事が8件件あります

全てのタグを見る

-σ^2 / 2 の真実に迫る – その1

· 約3分

先ずは準備運動

元金がV0V_0、利率をrr(年率)とした場合の単利と複利の場合それぞれで、 TT年後の元利合計を改めて考えてみましょう(何を今更、、、)。

単利の場合

年間に利払いがNN回あるとすると、1回あたりの利率はr/Nr/NとなったうえでTT年後までには金利がTNTN回つくのでの元利合計は、

V(T)=V0+V0(r/N+r/N++r/N)TNtimes=V0(1+rT)V(T)=V_0+V_0\cdot\underbrace{(r/N+r/N+\cdots+r/N)}_{TN \text{times}}=V_0\cdot(1+rT)    ---(1)

うーん文句なし。まったくもってあたりまえだ。

複利の場合

上記と同様に利率はr/Nr/Nとなるので、

V(T)=V0(1+r/N)(1+r/N)(1+r/N)TNtimes=V0(1+r/N)TNV(T)=V_0\cdot\underbrace{(1+r/N)(1+r/N)\cdots(1+r/N)}_{TN \text{times}}=V_0\cdot(1+r/N)^{TN}    ---(2)

となります。

次に上記の式でNN\rightarrow\inftyとしてみましょう。世にいう連続極限と呼ばれる操作です。 これは年間の利払い回数が無限大、つまり無限に短い期間に無限に小さい利率による利払いが発生し、それが無限に積み重なった場合を考えることです。

単利の場合

式(1)の中にNが含まれていないのでNN\rightarrow\inftyとしても結果は同じです。

複利の場合

式(2)でNN\rightarrow\inftyとするとどうなるか、結果だけ書きます(高校数学!)

V(T)=V0exp(rT) V(T)=V_0\cdot\exp(rT)    ---(3)

連続極限で指数関数になります。

と、まぁこんな話はどの金融の教科書にも最初のページに書いてあります。

さて上記のケースは全て金利が一定値rrという場合です。 そうではなく、金利が乱数だったらどうなるでしょうか。それが今回の記事で書きたいことです。

その2へ続く、、、

-σ^2 / 2 の真実に迫る – その2

· 約2分

発生する金利の利率が乱数だったらどうなるか?さすがに何の仮定もないと話が進まないので、

金利(年率)は、平均mm、標準偏差σ\sigmaの正規乱数

という場合を論じましょう。

ii回目の利払金利(年率)は、

ri=σzi+mr_i= {\sigma}z_i+m   ---(4)

と与えられます(ここで ziz_iは標準正規乱数)。

式(4)は年率なので、半期分、つまり年2回の利払いがある場合の各回の利率は、

ri(2)=σ2zi+m2r_i(2)=\frac{\sigma}{\sqrt{2}}{\cdot}{z_i}+\frac{m}{2}

と与えられます。一般に年NN回の利払いがある場合の各回の利率は、

ri(N)=σNzi+mNr_i(N)=\frac{\sigma}{\sqrt{N}}{\cdot}{z_i}+\frac{m}{N}

となります。

単利の場合

V(T)=V0+i=1TNri(N)=V0[1+i=1TN(σNzi+mN)]V(T)=V_0+\sum_{i=1}^{TN}r_i(N)=V_0\cdot\left[1+\sum_{i=1}^{TN}\left(\frac{\sigma}{\sqrt{N}}{\cdot}{z_i}+\frac{m}{N}\right)\right]

となります。

ここで「正規乱数の和は正規乱数になる」という、もはや犬でも知っている事実から、 V(T)V(T)も正規乱数であることが分かります。

分布の特徴は以下の通り。

  • 平均:

    E[V(T)]=V0+iTNmN=V0+mTE[V(T)]=V_0+\sum_{i}^{TN}\frac{m}{N}=V_0+mT
  • 分散:

    V[V(T)]=σ2Ni=1TNV[zi]=σ2TV[V(T)]=\frac{\sigma^2}{N}\cdot\sum_{i=1}^{TN}V[z_i]=\sigma^2 T

以上より、

V(T)N(V0+mT,σ2T)V(T){\sim}{N}(V_0+mT,\sigma^2T)

という分布になることが判明しました。

金利が毎回変わったとしても、この式中に年間利払回数NNが含まれないので、連続極限でも同一の結果(分布)になることがすぐに分かります。

さて次は複利の場合ですが、、、これがチトややこしいので、その3に続く

-σ^2 / 2 の真実に迫る – その3

· 約4分

複利の場合

V(T)=V0i=1TNri(N)=V0i=1TN(σNzi+mN+1)V(T)=V_0\cdot\prod_{i=1}^{TN}r_i(N)=V_0\cdot\prod_{i=1}^{TN}\left(\frac{\sigma}{\sqrt{N}}{\cdot} z_i+\frac{m}{N}+1\right)    ---(5)

となります。これが離散的な場合の結果であり、議論のスタート台となります。

こいつの連続極限を考えます。 知りたいことは、連続極限でV(T)V(T)の分布はどうなるか?です。

具体的には、

  • 分布形
  • 平均
  • 分散

といったことを調べます。

式(5)、、、掛け算のままだと取扱いが難しいので、両辺の対数をとりましょう。

lnV(T)=lnV0+i=1TNln(σNzi+mN+1){\ln}V(T)={\ln}V_0+\sum_{i=1}^{TN}\ln\left(\frac{\sigma}{\sqrt{N}}{\cdot}z_i+\frac{m}{N}+1\right)    ---(6)

ここに伝家の宝刀、テイラー展開

ln(x+1)=xx22+x33\ln(x+1)=x-\frac{x^2}{2}+\frac{x^3}{3}-\cdots    ---(7)

を適用します。定石通り1次の項だけとると式(6)は、

lnV(T)=lnV0+i=1TN(σNzi+mN){\ln}V(T)= {\ln}V_0+\sum_{i=1}^{TN}\left(\frac{\sigma}{\sqrt{N}}{\cdot}z_i+\frac{m}{N}\right)    ---(8)

となって、これは正に単利の場合と同一の式でキレイに正規乱数の和に収束することが分かりました。

つまり lnV(T){\ln}V(T) は正規分布になることが判明しました(覚えていますか?正規乱数の和は正規乱数です)。

さてこれで

対数をとったヤツが正規分布になった\longrightarrow元は対数正規分布

という事実から連続極限で V(T)V(T) は対数正規分布であると結論付けられました!

では平均と分散はどうなるでしょうか。

式(8)から、

  • E[lnV(T)]=lnV0+mTE\left[{\ln}V(T)\right]={\ln}V_0+mT
  • V[lnV(T)]=σ2TV\left[{\ln}V(T)\right]=\sigma^2 T

でメデタシメデタシ、、、としたいところですが、これが大間違い!。本記事の核心がここにあります。

話はテイラー展開(7)にまで戻ります。

2次の項まで考えます。すると

lnV(T)=lnV0+i=1TN(σNzi+mN)12i=1TN(σNzi+mN)2{\ln}V(T)={\ln}V_0+\sum_{i=1}^{TN}\left(\frac{\sigma}{\sqrt{N}}{\cdot}z_i+\frac{m}{N}\right)-\frac{1}{2}\sum_{i=1}^{TN}\left(\frac{\sigma}{\sqrt{N}}{\cdot}z_i+\frac{m}{N}\right)^2

となります。

この式から平均を計算してみます。\

E[lnV(T)]=lnV0+mT12iTN[σ2NE(zi2)+2σmNE(zi)+m2N2]E\left[{\ln}V(T)\right]={\ln}V_0+mT-\frac{1}{2}\sum_{i}^{TN}\left[\frac{\sigma^2}{N}E(z_i^2)+2\sigma m\sqrt{N}E(z_i)+\frac{m^2}{N^2}\right]

となって、ここで E(zi)=0,  E(zi2)=1E(z_i)=0,\;E(z_i^2)=1 という事実を適用すると、\

E[lnV(T)]=lnV0+mTσ22Tm2T2NE\left[{\ln}V(T)\right]={\ln}V_0+mT-\frac{\sigma^2}{2}T-\frac{m^2T}{2N}

となります。そして連続極限 NN\rightarrow\infty をとってみると、

E[lnV(T)]=lnV0+mTσ22TE\left[{\ln}V(T)\right]={\ln}V_0+mT-\frac{\sigma^2}{2}T

が得られます。気が付きましたか?テーラー展開の2次の近似項から NN  に依存しない定数が発生しているんですよ!

分散も同じ要領で計算してみると、こちらは2次の近似項は NN\rightarrow\infty で消えてなくなります。従って先ほど計算した通り、

V[lnV(T)]=σ2TV\left[\ln V(T)\right]=\sigma^2 T

です。

以上まとめると、 V(T)V(T) は対数正規分布であり、つまり lnV(T){\ln}V(T) は正規分布ということになって、

lnV(T)N(lnV0+mTσ22T,σ2T){\ln}V(T){\sim} N\left({\ln}V_0+mT-\frac{\sigma^2}{2}T,\sigma^2T\right)

となります。

標準正規分布に従う確率変数 ZZ を用いると、\

V(T)=V0exp[σTZ+(mσ22)T]V(T)=V_0\exp\left[\sigma\sqrt{T}Z+\left(m-\frac{\sigma^2}{2}\right)T\right]

と書き下せます。

最後にもう一言

 lnV(T){\ln}V(T) の平均が (mσ22)T\left(m-\frac{\sigma^2}{2}\right)T である。では V(T)V(T) そのものの平均は?

対数正規分布の公式から、

E[V(T)]=V0exp(mT)E\left[V(T)\right]=V_0\exp(mT)

が得られます。どうですか?実に合理的な答えになっていると思いませんか?(その1の式(3)と比較してみてください)。

長くなりましたが、あのキモチワルイ σ22-\frac{\sigma^2}{2} の出自が理解できましたか?

「独立」と「無相関」

· 約3分

複数の確率変数があるとき、それらがてんでバラバラ好き勝手に変動しているとき「無相関」という言葉を使います。 また「独立」という言葉もあります。
どちらも何となく似たような意味合いなので、混同している(というか気にしていない)場合もあるようですが、もちろんこの2つは異なる定義があります。
では「無相関」と「独立」の関係を見てみましょう。

「独立」なら「無相関」か?

その通り!「独立」の方が強いです。これはまぁ言葉通りのイメージではないでしょうか。「独立」しているのに何らかの関連性があるとは考えにくいです。

「無相関」なら「独立」か?

これは違います!。お互いに関連せずに変動していても「独立」でない場合があります。
これを理解するにはやはり「相関」を正しく理解することが必要です。
相関の定義は、

ρxy=1Ni=1N(xixˉ)(yiyˉ)\rho_{xy}=\frac{1}{N}\sum_{i=1}^{N}(x_i-\bar{x})(y_{i}-\bar{y})

です。細かい内容はさておき、重要なのは相関はデータ全体でひとつの値、つまり平均や標準偏差と同じデータの特徴を現す量であることです。
平均は分布の位置、標準偏差は分布の広がり具合です。

相関係数が 0 つまり無相関の場合とそうでない場合のデータプロットの違いが分かりましたか?無相関だと分布が全体としてまん丸になるになりますね。もっとも無相関だからといって分布がまん丸になるとは限りませんが、少なくともまん丸なら無相関です。
では下のようなプロットはどうでしょう?

このプロット (x, y) を生成する確率変数 X と Y は無相関です(まん丸で、どっちに傾けても変化なし)。しかしこの図をみて X と Y が勝手に変動していると認識する人はいないでしょう(たとえば X が 0付近なら、Y は ±1\pm 1 付近しか値がとれない)。

STDEVP? STDEV?

· 約8分

確率や統計にかかわる者が一度は悩む...否、一生悩む問題が 「標準偏差を求めるときって、STDEV を使うのか STDEVP を使うのか?」でしょう(STDEV や STDEVP は Excel の関数です。念のため)。 この2つの関数の違いは、

  • STDEVP:((データを2乗したものの平均)-(データの平均の2乗))/ データ数
  • STDEV:((データを2乗したものの平均)-(データの平均の2乗))/ (データ数 -- 1)

となっていて、教科書には「データ数から1を引くのは、自由度が1つ減ってどうのこうの...」という説明があり、分かったような分からないような...みなさんはこの2つの標準偏差の違いをしっかり理解していますか?

データの分布:STDEVP

こちらがいわゆる普通の標準偏差です。高校まではこちらの標準偏差しか学ばないでしょう。

あるデータが与えられました。データの羅列だけでは特徴はつかめません。そこでデータを集約して特徴をとらえようと考えます。代表的な特徴値としては、

  • データの典型的な値は?
  • データはどれくらいばらついている?

といったところではないでしょうか。データの典型的な値としては平均や中央値、そしてバラツキ度合いが STDEVP で求められる標準偏差です(定義から明らか)。
このように単に与えられたデータの特徴の1つとしての標準偏差が STDEVP です。

データの背後の真の分布:STDEV

1年後の株価を知りたい!これは人間誰しも持つ欲望でしょう。しかし、残念ながら株価は確率的に変動していて正確に知る方法はなさそうです。
株価は確率的に変動する、そして株価の変動(対数収益率)は正規分布に従うというのが一般的な知見です(最近はいろいろと異論があるようですが)。
時々刻々と記録されていく株価の変動は1つの正規分布、つまりある平均 mm とある標準偏差 σ\sigma をもつ正規分布(この2つで正規分布は決定できるというのは常識です)に従っていて、きっと分布の平均 mm と標準偏差 σ\sigma は神様が決めているんでしょう。 われわれにはその平均 mm と標準偏差 σ\sigma を正確に知る術はありませんから、なんとかして推定したいと考えます。

例えば株価の変動を毎日毎日記録していきます。1年ごとにデータをまとめて、それぞれのデータセットの平均 μi\mu_{i} と "標準偏差" σˉi\bar{\sigma}_{i} を求めます。 ここで "標準偏差" σˉi\bar{\sigma}_{i} は 何らかの方法 (具体的には STDEV で求めた標準偏差か STDEVP で求めた標準偏差)に従ったデータのバラツキ具合を表す値です。

結論から言うと、このときの標準偏差は STDEV で求ます。それの根拠は、、、

過去のある1年のデータから"標準偏差" を求めたとします。これは神様が決めた本当の正規分布(無限の過去から無限の未来まですべての株価のデータの分布)からたまたま選んだ300個ほどのデータから求めた"標準偏差" です。これを σˉ1\bar{\sigma}_1としましょう。 他の1年のデータから求めた"標準偏差" もたまたま選んだ300個ほどデータから求めたもので、これをσˉ2\bar{\sigma}_2 とします。 当然 σˉ1σˉ2\bar{\sigma}_{1}\neq \bar{\sigma}_{2} となっていることでしょう。どの1年のデータを持ってきても、そこから求められる "標準偏差" はデータセット毎に異なる値になります。
こうして幾つかのバラバラの値の"標準偏差" σˉi\bar{\sigma}_{i} が得られました。

まとめると

  • 本当の標準偏差 σ\sigma は絶対に知りえない
  • 観測から得られた σˉi\bar{\sigma}_{i} はばらついている

です。そこで

ばらついている σˉi\bar{\sigma}_{i} の平均が本当の標準偏差 σ\sigma になるようにしよう

という指針を採用します。
すると、STDEVP で求めた標準偏差はこの指針に従わないことが示され、STDEV の定義式で求められる標準偏差こそが新の標準偏差の推定値としてふさわしい(この指針を満たしている)ことが分かります!(細かい計算は適当な教科書を見てください)

結論:

  • データそのもののバラツキ具合:STDEVP
  • データの抽出元の真の標準偏差の推定値:STDEV

ということです。

サンプルエクセルシート

Stdevp_Stdev.xls
ダウンロード
真の分布は 平均 3(セル B1)、標準偏差 2(セル B2) の正規分布です。
そして、この分布からたまたま選んだ10個(C列からK列)のデータセットが全部で10,000セット(5行目から10004行目)あります。 そして各データセットの STDEVP(L列) と STDEV(M列) が計算されています。
STDEVP の平均(セル C2)と STDEV の平均(セル D2)を比較してみてください。STDEV の方が真の分布の推定値にふさわしいことが分かります(もっともっとデータセットを用意すれば STDEV の平均は真の標準偏差に近づいていきます)。

特異値分解入門(1)

· 約5分

互いに相関を持った正規乱数のセット(多変量正規乱数)を生成すること、これが金融工学のモンテカルロシミュレーションにおいて最も重要な技術となります。 NN 系列の正規乱数列を生成する場合を考えます。 これらの乱数列の相関構造を記述するのが相関行列と呼ばれる N×NN\times N の正方行列 C{\boldsymbol C} です。 最終的に

C=QQt{\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}

という形に分解することが目標となります(ここで t は転置行列です)。分解して得られた行列 Q{\boldsymbol Q} を使って実際に相関乱数を生成する方法については、これも一般的な数値計算の教科書に載っていますのでそちらをご覧ください。

すぐ後で説明するように、この分解にはコレスキー分解を使えと、多くの教科書で教えていますが、実務上コレスキー分解ではうまくいかない場合があります。
NtRand あるいは Numerical Technologies 社の製品では相関行列の分解に特異値分解という方法を採用しています。
これから3回にわたって、なぜコレスキー分解ではダメなのか、そしてなぜ特異値分解だとうまくいくのかを説明します。

相関行列について

相関行列 C{\boldsymbol C} は、その ij 成分が確率変数 XiX_{i} と XjX_{j} との相関係数 ρij\rho_{ij} である行列です。 この定義から対角成分 ρii\rho_{ii} は確率変数 XiX_{i} と自分自身との相関、つまり 1 となります。また、ρij=ρji\rho_{ij}=\rho_{ji} も定義から明らかです。
つまりこの行列は、

  • N×NN\times N正方行列
  • 成分は実数
  • 対角成分は 1
  • 対角成分以外は 1ρij1-1\leq\rho_{ij}\leq 1
  • 対称(ρij=ρji\rho_{ij}=\rho_{ji}

という性質を持っています。
更に相関係数の定義から、この行列は半正定値であることが示されます(半正定値って何だ?それは行列の固有値が全て0以上であること。固有値って何かって?それは勉強してください)。 つまり相関行列はその定義から半正定値実対称行列であるということです。

コレスキー分解 (Cholesky decomposition)

C=QQt{\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}

という分解を行う技として、どの教科書にも書いてあるのが**コレスキー分解という方法です。 どの教科書にも載っているので、ここでは説明しません。
この手法は正定値実対称行列を目的の形に分解してくれるので、正定値実対称行列である相関行列に対して使えば目標達成!メデタシ、メデタシとなる筋書きです。
しかし...
実はコレスキー分解は容易に失敗します。その代表的な状況が
観測データが少ない場合です。 確率変数の数が100個で、その観測データの数が100個未満の場合は、定義に従った正しい相関行列であってもコレスキー分解は失敗します。 なぜなら、この状況で作られる相関行列は必然的に正定値行列になるからです。コレスキー分解は正定値行列専用**です。
それでは、少なくとも正しく作られた相関行列では間違いなく

C=QQt{\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}

と分解してくれるより強力な方法はないでしょうか?

特異値分解入門(2)

· 約7分

ここで相関行列 C{\boldsymbol C} が満たしているはずの性質を今一度列挙してみましょう。

  • 成分 ρij\rho_{ij} は実数
  • 対称行列 ρij=ρji\rho_{ij}=\rho_{ji}
  • 半正定値
  • 対角成分 ρii=1\rho_{ii}=1
  • 対角成分以外 1ρij1-1\leq\rho_{ij}\leq 1

この相関行列 C{\boldsymbol C} を、

C=QQt{\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}

という形に分解するのが目的です(覚えてましたか?)
前回その方法としてコレスキー分解を紹介しましたが、どうも安心して使えなさそうだとわかりました。
次に紹介するのは固有値分解です。これは半正定値実対称行列を分解してくれます。
ではここからは線形代数の勉強です(証明などの詳細は省きます。教科書に必ず載っていますので調べてみてください)。

事実1

任意の実正方行列 A{\boldsymbol A} は、ある正方行列 U{\boldsymbol U} とこれまた正方行列、ただし対角成分のみしか値がない対角行列 W{\boldsymbol W} を使って、

A=UWU1{\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}

という形に分解される(これを行列のスペクトル分解といいます)。
ここで行列 U{\boldsymbol U} は行列 A{\boldsymbol A} の固有ベクトルを横に並べたもので、対角行列 W{\boldsymbol W} の成分は行列 A{\boldsymbol A} の固有値です。

事実2

更に行列 A{\boldsymbol A} が対称行列であった場合、行列の固有ベクトルが全て直交する。すると、

U1=Ut{\boldsymbol U}^{-1}={\boldsymbol U}^\text{t}

となる。というわけで

A=UWU1=UWUt{\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^\text{t}

事実3

更に更に行列 A{\boldsymbol A} が半正定値だとすれば、行列の固有値は全て0以上。つまり行列 W{\boldsymbol W} は、

(λ1000λ20vdots00λN)=(λ1000λ20vdots00λN)(λ1000λ20vdots00λN)=WW{\small \begin{pmatrix}\lambda_1&0&\cdots&0\\0&\lambda_2&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\lambda_{N}\end{pmatrix}=\begin{pmatrix}\sqrt{\lambda_1}&0&\cdots&0\\0&\sqrt{\lambda_2}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{\lambda_{N}}\end{pmatrix}\begin{pmatrix}\sqrt{\lambda_1}&0&\cdots&0\\0&\sqrt{\lambda_2}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{\lambda_{N}}\end{pmatrix}}={\boldsymbol W}^\prime {\boldsymbol W}^\prime

と出来るではないか!すると、

A=UWU1=UWUt=UWWUt{\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime {\boldsymbol W}^\prime {\boldsymbol U}^\text{t}

となります。ここで (AB)t=BtAt({\boldsymbol A}{\boldsymbol B})^\text{t}={\boldsymbol B}^\text{t}{\boldsymbol A}^\text{t} という事実(積の順番が入れ替わる)と、対角行列は転置しても変化しないこと(この場合 (W)t=W({\boldsymbol W}^\prime)^\text{t}={\boldsymbol W}^\prime)に注意すると、

A=UWU1=UWUt=UWWUt=UW(UW)t{\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime {\boldsymbol W}^\prime {\boldsymbol U}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime ({\boldsymbol U}{\boldsymbol W}^\prime)^\text{t}

やった!! UW=Q{\boldsymbol U}{\boldsymbol W}^\prime = {\boldsymbol Q} とおけば目標達成じゃないか!

A=QQt{\boldsymbol A}={\boldsymbol Q}{\boldsymbol Q}^\text{t}

しかし...
実務上(ここがポイント)、相関行列がいつもいつもタチの良い形であるとは限らないのです!もちろん実対称行列であることはどんな場合も間違いないですが(複素数の相関や、100×234100\times 234の相関行列を作れるものなら作ってみろ)、 半正定値を満たさない場合は往々にしてありえるのです。例えば...

  • データに欠損があって、欠損データを適当に埋める
  • 乱数列のうち、いくつかのペアについては相関が(観測とは無関係に)与えられている
  • 測定誤差

などの場合、相関係数に矛盾が出ることがあるのです。相関係数が矛盾するとはどういうことか、ひとつ最も大げさな例を示してみましょう。
以下の3変数の相関行列を見てください。

C=(111111111){\boldsymbol C}=\begin{pmatrix}1&1&-1\\1&1&1\\-1&1&1\end{pmatrix}

実数であるとか、対称であるとか相関行列の性質は確かに満たしているんですが...何がおかしいか分かりますか?

  • 第1変数と第2変数の相関係数 ρ12\rho_{12} は 1。第1変数と第2変数は完全相関で全く同じ値になる
  • 第2変数と第3変数の相関係数 ρ23\rho_{23} は 1。第2変数と第3変数は完全相関で全く同じ値になる

ってことは必然的に第1変数と第3変数は同じ値になるはずです。しかし与えられた相関係数 ρ13\rho_{13} は -1 になっている!つまり矛盾しているわけです。
要するに相関係数がありえない組み合わせになっている場合がありえるのです(さすがに例のような大げさなのはないでしょうけど)。このとき、相関行列は半正定値ではなくなってコレスキー分解が機能しなくなります。

実務上、相関行列が与えられたら(それに矛盾があろうとなかろうと)計算がストップしてしまうことは望ましくありません。とにかく安定して何らかの答えを出すことが求められます。 もちろん正しい相関行列(半正定値実対称行列)は正しく分解されることが前提です。

そこで NtRand では相関行列の分解に特異値分解(Singular Value Decomposition,SVD)という手法を採用しています。

特異値分解入門(3)

· 約7分

不正な相関行列 C{\boldsymbol C} が与えられると何が起きるが考えてみましょう。
前回の解説の "事実3" で悲劇が起きます。つまり、

(λ1000λ20vdots00λN)=(λ1000λ20vdots00λN)(λ1000λ20vdots00λN){\small \begin{pmatrix}\lambda_1&0&\cdots&0\\0&\lambda_2&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\lambda_{N}\end{pmatrix}=\begin{pmatrix}\sqrt{\lambda_1}&0&\cdots&0\\0&\sqrt{\lambda_2}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{\lambda_{N}}\end{pmatrix}\begin{pmatrix}\sqrt{\lambda_1}&0&\cdots&0\\0&\sqrt{\lambda_2}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{\lambda_{N}}\end{pmatrix}}

とする時になって、ハッと気が付くのです...「固有値が負だ!」と。
そしてコンピューターは無情にも「負の平方根は計算できません」というメッセージを残して息絶えます。
ちなみに第1回に載せた不正な相関行列、

C=(111111111){\boldsymbol C}=\begin{pmatrix}1&1&-1\\1&1&1\\-1&1&1\end{pmatrix}

の固有値は 1,2,2-1,2,2 となります。

この困難を打破してとにかく答えを出したい!そこでいよいよ特異値分解の登場です。

先ずは特異値分解で任意の行列(なんと正方行列でなくてもいい!)A{\boldsymbol A} はどのように分解されるか見てみましょう。

A=UWVt{\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol V}^\text{t}

となります。とりあえずここでは A{\boldsymbol A} として実対称正方行列であるところの相関行列 C{\boldsymbol C} だけを対象にして順に説明していきましょう。

行列 C{\boldsymbol C} から M=CCt{\boldsymbol M}={\boldsymbol C}{\boldsymbol C}^\text{t} なる行列を計算します。この行列は自動的に、

  • 正方行列
  • 対称行列
  • 半正定値行列

となるんです!(もちろんもとの行列が実数の行列ならば M{\boldsymbol M} も実数)。最後の特徴に注目してください。簡単に言うと、元の行列(C{\boldsymbol C})2乗することで負だった固有値があったとしてもそれが正になるということです。
というわけで行列 M{\boldsymbol M} は互いに直交する固有ベクトルを持ち、固有値は全て0以上となります。この固有ベクトルを並べた行列が U{\boldsymbol U} です。
次に、N=CtC{\boldsymbol N}={\boldsymbol C}^\text{t}{\boldsymbol C} なる行列を作ります。この行列は M{\boldsymbol M} の場合と全く同様に

  • 正方行列
  • 対称行列
  • 半正定値行列

となります。N{\boldsymbol N} の固有ベクトルを並べた行列が V{\boldsymbol V} です。
最後に W{\boldsymbol W} は行列 C{\boldsymbol C} の固有値の絶対値を対角に並べた対角行列になります。

ここまでの手順で分かるように、各行列 U{\boldsymbol U}V{\boldsymbol V}W{\boldsymbol W} は元の行列 C{\boldsymbol C} がどんなものであれ安定的に得ることができます(計算は止まらない!)。 前にも述べましたがこれが最も重要な点です。

とにもかくにも行列 C{\boldsymbol C} は、

C=UWVt{\boldsymbol C}={\boldsymbol U}{\boldsymbol W}{\boldsymbol V}^\text{t}

という形に分解されました。次に

W=(λ1000λ20vdots00λN)=(λ1000λ20vdots00λN)(λ1000λ20vdots00λN)=WW{\small {\boldsymbol W}=\begin{pmatrix}|\lambda_1|&0&\cdots&0\\0&|\lambda_2|&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&|\lambda_{N}|\end{pmatrix}=\begin{pmatrix}\sqrt{|\lambda_1|}&0&\cdots&0\\0&\sqrt{|\lambda_2|}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{|\lambda_{N}|}\end{pmatrix}\begin{pmatrix}\sqrt{|\lambda_1|}&0&\cdots&0\\0&\sqrt{|\lambda_2|}&\cdots&0\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{|\lambda_{N}|}\end{pmatrix}}={\boldsymbol W}^\prime {\boldsymbol W}^\prime

とします(対角成分は固有値の絶対値なので安心して平方根がとれる)。すると、

C=UWVt=UWWVt=UW(VW)t{\boldsymbol C}={\boldsymbol U}{\boldsymbol W}{\boldsymbol V}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime {\boldsymbol W}^\prime {\boldsymbol V}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime({\boldsymbol V}{\boldsymbol W}^\prime)^\text{t}

となります。

あと一息です

さてここで行列 C{\boldsymbol C} が正しい相関行列、つまり半正定値だとどうなるか。その場合は U=V{\boldsymbol U}={\boldsymbol V} となって特異値分解は前回の解説の結果と同じになります。 つまり Q=VW{\boldsymbol Q}={\boldsymbol V}{\boldsymbol W}^\prime とすることで、相関行列は C=QQt{\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t} となって目標達成です。
しかしそうでない場合はどうでしょう。その場合もやはり同じく Q=VW{\boldsymbol Q}={\boldsymbol V}{\boldsymbol W}^\prime とすることで手をうちませんか?もちろんこの場合は CQQt{\boldsymbol C}\neq {\boldsymbol Q}{\boldsymbol Q}^\text{t} ではありますが(UV{\boldsymbol U}\neq{\boldsymbol V}だから)...でも、

  • 安定的に答えが出る
  • 正しい相関行列の場合は正しい答えが出る

という望ましい性質を持ってます。
因みに正しい相関行列ではない場合にも、QQt{\boldsymbol Q}{\boldsymbol Q}^\text{t} の各成分の符号は元の行列 C{\boldsymbol C} のそれに等しくなります(正の相関が負になることはない)。

再び

C=(111111111){\boldsymbol C}=\begin{pmatrix}1&1&-1\\1&1&1\\-1&1&1\end{pmatrix}

を例に計算してみましょう。特異値分解の結果は、

C=(0.5770.7390.3480.5770.6710.4660.5770.0680.814)U(100020002)W(0.5770.5770.5770.7390.6710.0680.3480.4660.814)Vt{\boldsymbol C}=\underbrace{\begin{pmatrix}0.577&-0.739&-0.348\\-0.577&-0.671&0.466\\0.577&0.068&0.814\end{pmatrix}}_{{\boldsymbol U}}\underbrace{\begin{pmatrix}1&0&0\\0&2&0\\0&0&2\end{pmatrix}}_{{\boldsymbol W}^\prime}\underbrace{\begin{pmatrix}-0.577&0.577&-0.577\\-0.739&-0.671&0.068\\-0.348&0.466&0.814\end{pmatrix}}_{{\boldsymbol V}^\text{t}}

となります(UV{\boldsymbol U}\neq {\boldsymbol V}であることが分かりますか?第1固有ベクトルの符号が逆転しています!こうすることで、固有値の方を正に揃えているという訳です)。そして、

Q=VW=(0.5770.5770.5770.7390.6710.0680.3480.4660.814)(100020002){\boldsymbol Q}={\boldsymbol V}{\boldsymbol W}^\prime=\begin{pmatrix}-0.577&0.577&-0.577\\-0.739&-0.671&0.068\\-0.348&0.466&0.814\end{pmatrix}\begin{pmatrix}1&0&0\\0&\sqrt{2}&0\\0&0&\sqrt{2}\end{pmatrix}

となります。ここから QQt{\boldsymbol Q}{\boldsymbol Q}^\text{t} を計算すると、

QQt=(5/31/31/31/35/31/31/31/35/3){\boldsymbol Q}{\boldsymbol Q}^\text{t}=\begin{pmatrix}5/3&1/3&-1/3\\1/3&5/3&1/3\\-1/3&1/3&5/3\end{pmatrix}

となります(うーん、確かにC{\boldsymbol C}とはちょっと違うなぁ。残念)。

以上で相関行列の分解に特異値分解を用いる理由とその方法の解説を終了します。