アルキメデス型のコピュラは生成素という関数Φ（）を使って以下のように表現されている。
クレイトンコピュラの生成素(generator)は
その逆関数は
となる。
逆関数の中に生成素の関数を入れると
と定義し、その逆関数を
で表している場合がある。理由は分からないが、生成素と逆関数を逆転させて使っている。その場合にはアルキメデス型のコピュラの表記は以下のようにφ(φ-1)の順になる。
最初はどちらの表記が正しいのか分からなくて
非常に混乱してしまったが、生成素の定義をよく確かめれば実質的には同一の内容と分かった。
以下では一番最初の方の定義にしたがって生成素を
その逆関数を
で表記している。
ここからはMAXIMAを使って生成素とその逆関数の特徴をグラフで見てみる。 ０＜θ＜∞ の範囲に限定する。
生成素をw(t,θ）で表してグラフで示すと
単調減少の凸関数になっている。w(0,θ）＝∞
w(１,θ）＝０
w(t,θ)の逆関数 winv(t,θ)のグラフも見てみると
こちらも、単調減少の凸関数になっている。
winv(0,θ)=１，winv(∞,θ)=０
ガンマ分布で形状(shape)パラメータ＝１／θ、スケールパラメータ＝１とした標準ガンマ分布をラプラス変換するとクレイトンコピュラの生成素の逆関数と同じ数式が選られる。
クレイトンコピュラの乱数を発生させるときにこの関係が使われている。
に従って2変数ｘ、ｙについて、逆関数の中に生成素を代入してクレイトンコピュラを作ると
と表せる。これは2変数のクレイトンコピュラの同時分布を表している。
同時分布をｘ、ｙで偏微分してクレイトンコピュラの密度関数を求めると
合成関数の微分となるが、念のためにMAXIMAで答合わせをすると
実質的に
同じ結果を得た。
以上でクレイトンコピュラの２変数の同時分布と密度関数が導かれた。
相関係数は周知のように２変数ｘ、ｙの共分散covを各変数の標準偏差σで除して
cov(x,y)/（σｘ・σy）で計算する。これは便利であるが係数という単一の数値でｘとｙの関係性を単純化して表している。また、よく言われているようにｘとｙが無関係で独立していれば相関係数はゼロとなるが、逆は必ずしも真ならずで、ｘとｙに確かな関係性があっても相関係数がゼロとなることがある。相関係数は線形の関係性を発見するには便利であるが非線形の関係性については弱点がある。多くのテキストの例題でよく取り上げられるが,xが標準正規変数（平均０、分散１）で,ｙ=ｘ2
であるときにはｘとｙの相関係数を計算するとゼロとなる。計算過程は以下のとおりである。
ｘとｙの共分散がゼロとなるため、相関係数はゼロとなる。
Ｅ（ｘ3
）は奇関数の積分となり、うまい具合にゼロとなる。要するに、相関係数がゼロであったとしても安易に無関係と判断してはいけないと注意を促す有名な事例である。自分はすっかり忘れていたが、学生時代に使っていたテキストにも全く同じ例が書かれていたので驚いた。このテキストは絶版になっていると思うが、森田優三 統計数理入門 日本評論社 １９６８
（p158-159の例4.9参照）である。当時では文系向けの統計学と理系数学科で使う数理統計学の中間辺りを行くようなテキストで今でも大変に重宝している。相関係数は昔からなじんでいる指標で分かりやすくて便利であるが、金融の世界でよく見られるようなテールリスクといった確率分布の裾部分での依存関係を説明するのに不向きであり、２次積率である分散が無限に近くなるような裾の長い確率分布では相関係数を定義しがたいといった弱点が指摘もされている。コピュラはこのような問題を克服する一つの方法として利用されている。クレイトンコピュラを例にすれば左側（下方側）での依存関係をモデル化するのに適していると言われている。
以下は簡単な数値例でクレイトンコピュラを検討してみる。
まずｘ、ｙの２変数からなるクレイトンコピュラの乱数を生成する。乱数生成法は多数あるようだが、A Copula-Based Approach to
Accommodate Residential Self-Selection Effects in Travel Behavior
Modeling Chandra R. Bhat
The University of Texas at
Austin で紹介されている方法をそのままＲコードにして使ってみた。小生にはこの論文を吟味するほどの能力は無いが、自分なりのレベルで検証が出来たので納得している。それについては後ほど触れることにする。クレイトンコピュラではパラメータ θを指定する必要があるのでθ＝５の場合の乱数のペアを５０個、生成して散布図にすると以下のようになる。
θ＝５の場合、全体がラッパ状に散布して左側下方では依存関係が高くなる傾向が見られる。
先ほど２変数のクレイトンコピュラの密度関数を求めたので、それを使って密度関数をグラフにすると以下のようになる。左側は等高線を描いている。
２変数のクレイトンコピュラの同時分布をグラフにすると以下のようになる。左側は等高線を描いている。
統計ソフトＲにはコピュラに係るパッケージが多数あるがライブラリーcopulaを使うと上記の作業は簡単にできる。
最初にθ＝５の乱数のペアを５０個、生成する。
このパッケージを利用して生成した乱数を先ほど作った乱数と比較すると同じよになっている。
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
多分、パッケージの乱数は先の論文で紹介されていた方法と同じ方法で生成されているのではないかと推測した。
密度関数をグラフにすると
同時分布は以下のようになる。
パラメータθをゼロに近い０．０１、と５，１０の３つのケースについて分布を比べてみる。θが大きくなると依存度が大きくなり∞になると完全依存となる。逆に０に近づくほど依存度は減少していく。θが５の場合は左側（下方部）での依存関係は大きいが右側（上方部）では依存関係は小さくなっている。
θ＝０．０１のケース
θ＝５ のケース
θ＝１０ のケース
拙文の作成に使ったＲコードは以下のとおりである。
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 コピュラ乱数 散布図 θ =５ ")))),
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 コピュラ乱数 散布図 θ =５ ")))),
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))