クレイトン コピュラ(Clayton copula)の自由研究

数式が平易そうに見えたのでクレイトン コピュラの初歩を少し学習してみた。しかし実際に取り組んでみたものの、けっこう難渋してしまった。最初のつまずきはアルキメデス型コピュラの数式の表記法が論者によって異なっていたりして頭が混乱したことだ。しかし、家庭教師のMAXIMA先生の支援を得ながら、なんとか初歩の初歩といえるところをまとめてみた。

コピュラは一様分布である周辺分布を結合させて同時分布を導き出す接合関数であるが、どんな関数でも接合関数の役割を果たせるわけではなく、連続性とか共単調(comonotone)とか単調変換の下での不変性とか様々な数学的な特性を備えていることが必要となる。このような数学的特性を備えた接合関数としては正規コピュラ、tコピュラ、アルキメデス型コピュラなどが知られている。クレイトンコピュラはアルキメデス型コピュラに属するもので統計学でなじみ深いガンマ分布とも関係があるので、その辺を手がかりに学習してみた。

クレイトンコピュラについて

数学的に単純な形式で表せる接合関数のクラスの一つで、数学者たちにより開発されてきたのがアルキメデス型コピュラと呼ばれているようだ。クレイトンコピュラはアルキメデス型コピュラに属するコピュラで生成素(generator)とその逆関数で表すことが出来る。アルキメデス型コピュラの話になると生成素(generator)という用語が登場してくるが、文字通りに、コピュラを生成する関数で、連続性とか凸関数とか数学的に望ましい特性を備えているようだ。

アルキメデス型のコピュラは生成素という関数Φ()を使って以下のように表現されている。

アルキメデス型コピュラ

クレイトンコピュラの生成素(generator)は

クレイトンコピュラの生成素

その逆関数は
クレイトンコピュラの生成素の逆関数
となる。

逆関数の中に生成素の関数を入れると

逆関数と元の関数の関係

となり、恒等写像の関係が確かめられる。

様々な資料で アルキメデス型コピュラの解説を読むと、論者によっては生成素(generator)を

と定義し、その逆関数を

で表している場合がある。理由は分からないが、生成素と逆関数を逆転させて使っている。その場合にはアルキメデス型のコピュラの表記は以下のようにφ(φ-1)の順になる。



最初はどちらの表記が正しいのか分からなくて 非常に混乱してしまったが、生成素の定義をよく確かめれば実質的には同一の内容と分かった。

以下では一番最初の方の定義にしたがって生成素を

クレイトンコピュラの生成素
その逆関数を

クレイトンコピュラの生成素の逆関数

で表記している。


クレイトンコピュラの生成素(generator)について

ここからはMAXIMAを使って生成素とその逆関数の特徴をグラフで見てみる。 0<θ<∞  の範囲に限定する。


生成素をw(t,θ)で表してグラフで示すと



単調減少の凸関数になっている。w(0,θ)=∞

w(1,θ)=0

w(t,θ)の逆関数 winv(t,θ)のグラフも見てみると



こちらも、単調減少の凸関数になっている。


winv(0,θ)=1,winv(∞,θ)=0



ガンマ分布で形状(shape)パラメータ=1/θ、スケールパラメータ=1とした標準ガンマ分布をラプラス変換するとクレイトンコピュラの生成素の逆関数と同じ数式が選られる。

ガンマ分布

クレイトンコピュラの乱数を発生させるときにこの関係が使われている。


2変数クレイトンコピュラの分布関数と密度関数

アルキメデス型コピュラ

に従って2変数x、yについて、逆関数の中に生成素を代入してクレイトンコピュラを作ると


と表せる。これは2変数のクレイトンコピュラの同時分布を表している。

同時分布をx、yで偏微分してクレイトンコピュラの密度関数を求めると


合成関数の微分となるが、念のためにMAXIMAで答合わせをすると



実質的に 同じ結果を得た。

以上でクレイトンコピュラの2変数の同時分布と密度関数が導かれた。


相関係数の弱点について

相関係数は周知のように2変数x、yの共分散covを各変数の標準偏差σで除して
cov(x,y)/(σx・σy)で計算する。これは便利であるが係数という単一の数値でxとyの関係性を単純化して表している。また、よく言われているようにxとyが無関係で独立していれば相関係数はゼロとなるが、逆は必ずしも真ならずで、xとyに確かな関係性があっても相関係数がゼロとなることがある。相関係数は線形の関係性を発見するには便利であるが非線形の関係性については弱点がある。多くのテキストの例題でよく取り上げられるが,xが標準正規変数(平均0、分散1)で,=2 であるときにはxとyの相関係数を計算するとゼロとなる。計算過程は以下のとおりである。


相関係数の弱点

xとyの共分散がゼロとなるため、相関係数はゼロとなる。
E(x
3  )は奇関数の積分となり、うまい具合にゼロとなる。要するに、相関係数がゼロであったとしても安易に無関係と判断してはいけないと注意を促す有名な事例である。自分はすっかり忘れていたが、学生時代に使っていたテキストにも全く同じ例が書かれていたので驚いた。このテキストは絶版になっていると思うが、森田優三  統計数理入門 日本評論社 1968 (p158-159の例4.9参照)である。当時では文系向けの統計学と理系数学科で使う数理統計学の中間辺りを行くようなテキストで今でも大変に重宝している。相関係数は昔からなじんでいる指標で分かりやすくて便利であるが、金融の世界でよく見られるようなテールリスクといった確率分布の裾部分での依存関係を説明するのに不向きであり、2次積率である分散が無限に近くなるような裾の長い確率分布では相関係数を定義しがたいといった弱点が指摘もされている。コピュラはこのような問題を克服する一つの方法として利用されている。クレイトンコピュラを例にすれば左側(下方側)での依存関係をモデル化するのに適していると言われている。


Rを使って2変数クレイトンコピュラを計算してみる


以下は簡単な数値例でクレイトンコピュラを検討してみる。
まずx、yの2変数からなるクレイトンコピュラの乱数を生成する。乱数生成法は多数あるようだが、A Copula-Based Approach to Accommodate Residential Self-Selection Effects in Travel Behavior Modeling Chandra R. Bhat
The University of Texas at Austin で紹介されている方法をそのままRコードにして使ってみた。小生にはこの論文を吟味するほどの能力は無いが、自分なりのレベルで検証が出来たので納得している。それについては後ほど触れることにする。クレイトンコピュラではパラメータ θを指定する必要があるのでθ=5の場合の乱数のペアを50個、生成して散布図にすると以下のようになる。

2変数 クレイトンコピュラ 散布図

θ=5の場合、全体がラッパ状に散布して左側下方では依存関係が高くなる傾向が見られる。

 

先ほど2変数のクレイトンコピュラの密度関数を求めたので、それを使って密度関数をグラフにすると以下のようになる。左側は等高線を描いている。

クレイトンコピュラ確率密度


2変数のクレイトンコピュラの同時分布をグラフにすると以下のようになる。左側は等高線を描いている。

クレイトンコピュラ同時分布

 

統計ソフトRにはコピュラに係るパッケージが多数あるがライブラリーcopulaを使うと上記の作業は簡単にできる。
最初にθ=5の乱数のペアを50個、生成する。

クレイトンコピュラ散布図

 

このパッケージを利用して生成した乱数を先ほど作った乱数と比較すると同じよになっている。

 r1 u2
[1,] 0.4512674 0.4867697
[2,] 0.7837798 0.8383791
[3,] 0.7096822 0.4675296
[4,] 0.3817443 0.4097074
[5,] 0.6363238 0.5917164

> head(spc)
[,1] [,2]
[1,] 0.4512674 0.4867697
[2,] 0.7837798 0.8383791
[3,] 0.7096822 0.4675296
[4,] 0.3817443 0.4097074
[5,] 0.6363238 0.5917164

 

多分、パッケージの乱数は先の論文で紹介されていた方法と同じ方法で生成されているのではないかと推測した。

 

密度関数をグラフにすると

クレイトンコピュラ密度

 

同時分布は以下のようになる。

クレイトンコピュラ同時分布

 

パラメータθをゼロに近い0.01、と5,10の3つのケースについて分布を比べてみる。θが大きくなると依存度が大きくなり∞になると完全依存となる。逆に0に近づくほど依存度は減少していく。θが5の場合は左側(下方部)での依存関係は大きいが右側(上方部)では依存関係は小さくなっている。

 

θ=0.01のケース

θ=0.01のケース

 

θ=5 のケース

θ=5 のケース

 

θ=10 のケース

 

θ=10 のケース

 

 

拙文の作成に使ったRコードは以下のとおりである。

library(copula)

set.seed(2021)
theta<-5
r1<-runif(50)
r2<-runif(50)

#####clayton copula 乱数の生成

u2<-(r1^(-theta)*(r2^(-theta/(1+theta))-1)+1)^(-1/theta)
rru2<-cbind(r1,u2)

##コピュラ 散布図
plot(rru2,main=eval(substitute(expression(paste("Clayton コピュラ乱数 散布図 θ =5 ")))),
col=rainbow(50),pch=20,cex=3)

#コピュラ 密度関数
par(mfrow=c(1,2))
x <- seq(0,1,length=50)
y <- x
t <- theta
claydist <- function(x,y,t) {
(1/t+1)*t*(x^(-t-1))*(y^(-t-1))*((1/y^t)+(1/x^t)-1)^(-1/t-2)
}
z <- outer(x,y,claydist,t)
persp(x, y, z, col=cm.colors(1200), tick='detailed', theta=10, phi=10,zlim=c(0,70))
contour(x,y,z,nlev=50 ,method = "flattest",labcex=0.5)

##コピュラ 分布
cprru2<-(r1^(-theta)+u2^(-theta)-1)^(-1/theta)
cprru2

x <- seq(0,1,length=50)
y <- x
t <- theta
claycumdist <- function(x,y,t) {
(x^(-t)+y^(-t)-1)^(-1/t)}

zz <- outer(x,y,claycumdist,t)
persp(x, y, zz, col=cm.colors(1200), tick='detailed', theta=10, phi=10)
contour(x,y,zz,nlev=30)

par(mfrow=c(1,1))

##パッケージ copula 利用
theta=5
cx<-claytonCopula(theta,dim=2)
set.seed(2021)
spc<-rCopula(50,cx)

par(mfrow=c(1,1))

##コピュラ 散布図
plot(spc,main=eval(substitute(expression(paste("Clayton コピュラ乱数 散布図 θ =5 ")))),
col=rainbow(50),pch=20,cex=3)

par(mfrow=c(1,2))

persp(claytonCopula(theta), dCopula,phi=10,theta=60,col=cm.colors(1200))
contour(cx, dCopula,nlev=40)

persp(claytonCopula(theta), pCopula,phi=10,theta=20,col=cm.colors(1200))
contour(cx, pCopula,nlev=30)

##生成した乱数の比較
head(rru2)
head(spc)

## パッケージ copula で作成された乱数と同じなので多分
#上記の方法と同じ方法を使っているのではないかと推測できる。

par(mfrow=c(1,3))

theta<-c(0.01,5,10)
rannum<-100 #生成する乱数の数
for(i in 1:3){
set.seed(2021)
cxx<-claytonCopula(theta[i],dim=2)
spc2<-rCopula(rannum,cxx)

plot(spc2,main=eval(substitute(expression(paste("Clayton コピュラ乱数散布図 ",theta," = ",j)),
list(j=as.character(theta[i])))),col=rainbow(100),pch=20,cex=2)
persp(cxx,dCopula,main=eval(substitute(expression(paste("Clayton コピュラ密度関数 ", theta," = ",j)),
list(j=as.character(theta[i])))),col=cm.colors(1200),phi=20,theta=10)
persp(cxx,pCopula,main=eval(substitute(expression(paste("Clayton コピュラ分布関数 ",theta," = ",j)),
list(j=as.character(theta[i])))),col=cm.colors(1200),phi=20,theta=10)

}

par(mfrow=c(1,1))

 

コピュラ(接合関数)の初歩への道

 

金融工学を初等数学で 目次