本文原始地址:
https://www.r-bloggers.com/how-to-fit-a-copula-model-in-r-heavily-revised-part-1-basic-tools/
本文只是自学copula的记录文件,并非逐字逐句翻译,详情请见原文.
作者一年前发过个帖子,因为当时相关的文档少,难度大,对新手非常之不友好.
现在要更新一下这个老帖子,原因有二:
- 对原帖不满意
- 通过观察原帖的访问者,似乎大家都想知道的是该如何入门,如何生成一个copula模型.那好吧,这个可以有.
此新帖包含两个部分:
给新手介绍使用R进行copulas的基本工具和函数.
介绍如何选择copula,拟合过程,评估及其它问题,最后是一个案例.
本文是介绍如何用R处理copulas,言下之意是对copulas,假设你已经知道起码椭圆和阿基米德copulas了.如果你不知道,那推荐看看这本书.此书对新手很合适,条理和解释清晰.
椭圆copulas
如果定义$F$为椭圆分布的多元累积分布函数,$F_i$为第i个边缘分布,$F^{-1}$为其反函数,那么就可以构建如下椭圆copula:
$$ C(u_1, … , u_n) = F (F_1^{-1}(u_1) + … + F_1^{-1}(u_1)) $$
阿基米德copulas
阿基米德copula可以用$phi$函数(生成器函数)来生成:
$$ C(u_1, … , u_n) = phi ^{-1} (phi(u_1) + … + phi(u_n)) $$
比如,可以用如下生成器函数来生成一个Gumbel copula:
$$-( ln x )^{theta}$$
其反函数$$ exp( – y^{frac{1}{theta}} ) $$
最常见的阿基米德copulas包括Frank,Gumbel和Clayton.这三个也是本文所要用到的.
加载必要的库
我们首先加载需要用到的库文件,设置好random seed以复用.
# Copula package
library(copula)
# Fancy 3D plain scatterplots
library(scatterplot3d)
# ggplot2
library(ggplot2)
# Useful package to set ggplot plots one next to the other
library(grid)
set.seed(235)
椭圆copulas和阿基米德copulas热身--用已知参数生成copula
如你所知,椭圆与阿基米德couplas有很好的数学性质,形状和公式.用它们来热身再好不过.假设你已经知道了参数,下面的步骤将告诉你如何生成二元的正态copula和t copula.
# Generate a bivariate normal copula with rho = 0.7
normal <- normalCopula(param = 0.7, dim = 2)
# Generate a bivariate t-copula with rho = 0.8 and df = 2
stc <- tCopula(param = 0.8, dim = 2, df = 2)
这里的参数已经很明显表达了含义.正态copula有两个参数:copula的维度,这里等于2,以及可以从数据中估计出来的$rho$参数,关于它我会在第二部分讲到.
t copula除了有上述的两个参数以外,还包括一个df参数,它能决定copula的形状,其值和$rho$一样也可从数据中估计出.copula包提供了相应的函数来生成需要的copula模型:一般的形式为"name"+"Copula()",用于生成阿基米德copulas.阿基米德copulas只需要一个参数$theta$,具体如下:
# Build a Frank, a Gumbel and a Clayton copula
frank <- frankCopula(dim = 2, param = 8)
gumbel <- gumbelCopula(dim = 3, param = 5.6)
clayton <- claytonCopula(dim = 4, param = 19)
# Print information on the Frank copula
print(frank)
用mvdc在指定边缘分布下生成copula
如你所知，copulas的吸引力在于，能将边缘分布从相关性结构中剥离出来，独立的形成模型。这一特性能够极大地简化对复杂相关性行为进行建模的过程，因为边缘分布比联合分布要简单得多。比如，你要对20多个对象进行建模（比如对电器中的20多个电容），你会发现直接建模很困难，而用边缘分布就简单多了。
copula包提供了一系列好用的函数（mvdc，dMvdc，pMvdc，rMvdc）用于建立多元分布的copula.如下:
# Select the copula
cp <- claytonCopula(param = c(3.4), dim = 2)
# Generate the multivariate distribution (in this case it just bivariate) with normal and t marginals)
multivariate_dist <- mvdc(copula = cp,
margins = c("norm","t"),
paramMargins = list(list(mean = 2, sd = 3),
list(df = 2)))
print(multivariate_dist)
从copula随机抽样
一旦生成了copula,就可以很容易的从中生成随机数.针对single copula用的是rCopula方法,针对多元分布用的是rMvdc方法,如下:
# Generate random samples
fr <- rCopula(2000,frank)
gu <- rCopula(2000,gumbel)
cl <- rCopula(2000,clayton)
# Plot the samples
p1 <- qplot(fr[,1], fr[,2], colour = fr[,1], main = "Frank copula random samples theta = 8", xlab = "u", ylab = "v")
p2 <- qplot(gu[,1], gu[,2], colour = gu[,1], main = "Gumbel copula random samples theta = 5.6", xlab = "u", ylab = "v")
p3 <- qplot(cl[,1], cl[,2], colour = cl[,1], main = "Clayton copula random samples theta = 19", xlab = "u", ylab = "v")
# Define grid layout to locate plots and print each graph^(1)
pushViewport(viewport(layout = grid.layout(1, 3)))
print(p1, vp = viewport(layout.pos.row = 1, layout.pos.col = 1))
print(p2, vp = viewport(layout.pos.row = 1, layout.pos.col = 2))
print(p3, vp = viewport(layout.pos.row = 1, layout.pos.col = 3))
samples <- rMvdc(2000, multivariate_dist)
scatterplot3d(samples[,1], samples[,2], color = "blue", pch = ".")
PDF和CDF
计算copula的PDF(概率密度函数)和CDF(累积分布函数)可能对将来比较有用.R包提供了计算密度和CDF的函数,命名方式也保持了一致性:pCopula计算累积分布,dCopula计算密度.注意rCopula是用来抽样的.针对多元分布的也是类似命名习惯,pMvdc,dMvdc,rMvdc.如下:
# Generate the normal copula and sample some observations
coef_ <- 0.8
mycopula <- normalCopula(coef_, dim = 2)
u <- rCopula(2000, mycopula)
# Compute the density
pdf_ <- dCopula(u, mycopula)
# Compute the CDF
cdf <- pCopula(u, mycopula)
# Generate random sample ovservations from the multivariate distributions
v <- rMvdc(2000, multivariate_dist)
# Compute the density
pdf_mvd <- dMvdc(v, multivariate_dist)
# Compute the CDF
cdf_mvd <- pMvdc(v, multivariate_dist)
现在你只需要plot就能显示密度和CDF了.
图像
当需要作演示时,图像就很重要了.不过,如果涉及到高维的情况,要作图是很有挑战性的.而copula经常就要与高维相关性模型打交道.假设你现在的copula是二维的,你可以用不同方法对密度和累积分布函数作图.比如你可以画一个密度散点图或轮廓图.
copula包包含了一些有用的绘制轮廓图的方法,可以作为传统3d绘图的很好替代.尤其是在mvdc控件和自定义边缘分布的条件下,轮廓图可以很好的反映密度函数的形状.如下:
\par(mfrow = c(1, 3))
# 3D plain scatterplot of the density, plot of the density and contour plot
scatterplot3d(u[,1], u[,2], pdf_, color="red", main="Density", xlab ="u1", ylab="u2", zlab="dCopula", pch=".")
persp(mycopula, dCopula, main ="Density")
contour(mycopula, dCopula, xlim = c(0, 1), ylim=c(0, 1), main = "Contour plot")
par(mfrow = c(1, 3))
# 3D plain scatterplot of the CDF, plot of the CDF and contour plot
scatterplot3d(u[,1], u[,2], cdf, color="red", main="CDF", xlab = "u1", ylab="u2", zlab="pCopula",pch=".")
persp(mycopula, pCopula, main = "CDF")
contour(mycopula, pCopula, xlim = c(0, 1), ylim=c(0, 1), main = "Contour plot")
# 3D plain scatterplot of the multivariate distribution
par(mfrow = c(1, 2))
scatterplot3d(v[,1],v[,2], pdf_mvd, color="red", main="Density", xlab = "u1", ylab="u2", zlab="pMvdc",pch=".")
scatterplot3d(v[,1],v[,2], cdf_mvd, color="red", main="CDF", xlab = "u1", ylab="u2", zlab="pMvdc",pch=".")
persp(multivariate_dist, dMvdc, xlim = c(-4, 4), ylim=c(0, 2), main = "Density")
contour(multivariate_dist, dMvdc, xlim = c(-4, 4), ylim=c(0, 2), main = "Contour plot")
persp(multivariate_dist, pMvdc, xlim = c(-4, 4), ylim=c(0, 2), main = "CDF")
contour(multivariate_dist, pMvdc, xlim = c(-4, 4), ylim=c(0, 2), main = "Contour plot")
对多元分布的密度和CDF作图可能更有趣.注意散点图和密度图可能会比较误导.这也是为什么很多时候我宁可看轮廓图.
常见阿基米德copulas的密度和CDF的图像比较
把Frank，Gumbel和Clayton的密度拿来比较看看应该很有趣。通过copula包内的方法可以很容易做到这一点。注意，我对参数的选择是随机的。
frank <- frankCopula(dim = 2, param = 3)
clayton <- claytonCopula(dim = 2, param = 1.2)
gumbel <- gumbelCopula(dim = 2, param = 1.5)
par(mfrow = c(1, 3))
# Density plot
persp(frank, dCopula, main ="Frank copula density")
persp(clayton, dCopula, main ="Clayton copula density")
persp(gumbel, dCopula, main ="Gumbel copula density")
# Contour plot of the densities
contour(frank, dCopula, xlim = c(0, 1), ylim=c(0, 1), main = "Contour plot Frank")
contour(clayton, dCopula, xlim = c(0, 1), ylim=c(0, 1), main = "Contour plot Clayton")
contour(gumbel, dCopula, xlim = c(0, 1), ylim=c(0, 1), main = "Contour plot Gumbel")
