前言：本咩还曾经翻译过两篇copula相关的基础文章，还有一个时变copula的代码分享，都在我的微博上有。但是那个微博我现在不怎么用了。感兴趣的同学可搜搜微博名：补补补牢

今天分享的这些代码，原文请见：modelling-dependence-with-copulas-in-r，本文仅供本咩学习用途。

我不会被打倒

第一部分：copula是如何运行的

给定相关矩阵，从多元正态分布中抽样，利用MASS包。

library(MASS) #加载MASS包

set.seed(100)#随机数种子

m <- 3 #3维变量

n <- 2000 #样本量

sigma <- matrix(c(1, 0.4, 0.2,

                  0.4, 1, -0.8,

                  0.2, -0.8, 1),

                nrow=3)#相关矩阵，可以看出相关系数分别为0.4，0.2和0.8

z <- mvrnorm(n,mu=rep(0, m),Sigma=sigma,empirical=T)#抽样

library(psych)

cor(z,method='spearman')#在椭圆copula下，这个可以用什么相关系数都不打紧

pairs.panels(z)#画图

样本的特征符合我们的预期

u <- pnorm(z)#因为已经知道分布了，所以分容易得到F(x)

pairs.panels(u)

注意这里的每个变量都处于【0，1】之间，相关系数没变

library(rgl)#加载rgl包，可以画出变量u的3d形式    

plot3d(u[,1],u[,2],u[,3],pch=20,col='navyblue')

三个变量都在0-1之间

x1 <- qgamma(u[,1],shape=2,scale=1)#选择边际分布，得到向相应分布的分位点

x2 <- qbeta(u[,2],2,2)

x3 <- qt(u[,3],df=5)

plot3d(x1,x2,x3,pch=20,col='blue')

这就是生成的一个简单的依赖结构

通过上面的多元正态样本，我们建立了一个具有预期依赖结构的样本，但是却具有任意的边缘分布。

df <- cbind(x1,x2,x3)

pairs.panels(df)

cor(df,meth='spearman')

不同边缘分布

作者想强调的是他这里的边缘分布的选择是任意的，无论怎样的分布都能构造出他想要的相依结构。这就是copula函数的具体运行过程。

上面的所有过程都可以用【copula】包简单完成。

library(copula)

set.seed(100)

myCop <- normalCopula(param=c(0.4,0.2,-0.8), dim = 3, dispstr = "un")

myMvd <- mvdc(copula=myCop, margins=c("gamma", "beta", "t"),

              paramMargins=list(list(shape=2, scale=1),

                                list(shape1=2, shape2=2),

                                list(df=5)) )

这里通过指定copula确定了相依性，并指定了边缘分布。

Z2 <- rmvdc(myMvd, 2000)#抽样

colnames(Z2) <- c("x1", "x2", "x3")

pairs.panels(Z2)

这里的模拟数据当然非常接近之前模拟的数据，并显示在下面的pairplot中： 

与之前的结果类似

第二部分：一个小例子

这里用到的数据文件可以在开头那个链接里找到。

#读入变量

cree <- read.csv('cree_r.csv',header=F)$V2

yahoo <- read.csv('yahoo_r.csv',header=F)$V2

#在直接使用copula模型之前，查看一下两个变量的相关性。

plot(cree,yahoo,pch='.')

abline(lm(yahoo~cree),col='red',lwd=1)

cor(cree,yahoo,method='spearman')

[1] 0.4023584

有一个非常轻微的正相关，收益率序列看着有点像随机的

至于copula模型的选择，一般要根据数据特征，金融数据一般都会选择t-copula模型进行拟合，但是也可以使用R的函数BiCopSelect让它帮你选择，可以根据不同的标准（最大似然，AIC or BIC）

library(VineCopula)

u <- pobs(as.matrix(cbind(cree,yahoo)))[,1]

v <- pobs(as.matrix(cbind(cree,yahoo)))[,2]

selectedCopula <- BiCopSelect(u,v,familyset=NA)

selectedCopula

$family

[1] 2 #确实是t-copula比较合适

$par#tcopula是有两个参数的

[1] 0.4356302 

$par2

[1] 3.844534

作者想再用copula模型check一下，再看一下参数的拟合情况【实际上是一模一样的】

t.cop <- tCopula(dim=2)

set.seed(500)

m <- pobs(as.matrix(cbind(cree,yahoo)))#将数据变换为伪观测值【0，1】上的

fit <- fitCopula(t.cop,m,method='ml')

coef(fit)

  rho.1      df

0.43563 3.84i4e53'g 

最后这个部分是结果显示，查看一下密度函数，t-copula是对称的

rho <- coef(fit)[1]

df <- coef(fit)[2]

persp(tCopula(dim=2,rho,df=df),dCopula)

t-copula的密度函数

从t-copula中抽样画图

u <- rCopula(3965,tCopula(dim=2,rho,df=df))#抽样3965个

plot(u[,1],u[,2],pch='.',col='blue')

cor(u,method='spearman')

          [,1]      [,2]

[1,] 1.0000000 0.3972454

[2,] 0.3972454 1.0000000

抽样图，两个变量看着像独立情况，因为相关系数确实比较低

t-copula强调极端结果：它通常适用于极端值（分布的尾部）高度相关的现象建模。 

三、对边缘分布建模

假设两个收益率序列服从正态分布

cree_mu <- mean(cree)

cree_sd <- sd(cree)

yahoo_mu <- mean(yahoo)

yahoo_sd <- sd(yahoo)

#看下实际数据和正态分布的情况

hist(cree,breaks=80,main='Cree returns',freq=F,density=30,col='cyan',ylim=c(0,20),xlim=c(-0.2,0.3))

lines(seq(-0.5,0.5,0.01),dnorm(seq(-0.5,0.5,0.01),cree_mu,cree_sd),col='red',lwd=2)

legend('topright',c('Fitted normal'),col=c('red'),lwd=2)

hist(yahoo,breaks=80,main='Yahoo returns',density=30,col='cyan',freq=F,ylim=c(0,20),xlim=c(-0.2,0.2))

lines(seq(-0.5,0.5,0.01),dnorm(seq(-0.5,0.5,0.01),yahoo_mu,yahoo_sd),col='red',lwd=2)

legend('topright',c('Fitted normal'),col=c('red'),lwd=2)

文中为了简化，用的正态，实际上有高峰特征

#构造函数并抽样
copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=df), margins=c("norm","norm"),

                    paramMargins=list(list(mean=cree_mu, sd=cree_sd),

                                      list(mean=yahoo_mu, sd=yahoo_sd)))

sim <- rmvdc(copula_dist, 3965)

#画图对比

plot(cree,yahoo,main='Returns')

points(sim[,1],sim[,2],col='red')

legend('bottomright',c('Observed','Simulated'),col=c('black','red'),pch=21)

模拟和观察结果

从上图看到，tcopula的模拟结果已经很接近真实的观察样本，除了样本有一些极端值没有被覆盖到，这可能要通过进一步的校准。

三、通过改变t-copula的参数值，copula的模拟值将会发生怎样的变化

set.seed(4258)

copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=1), margins=c("norm","norm"),

                    paramMargins=list(list(mean=cree_mu, sd=cree_sd),

                                      list(mean=yahoo_mu, sd=yahoo_sd)))

sim <- rmvdc(copula_dist, 3965)

plot(cree,yahoo,main='Returns')

points(sim[,1],sim[,2],col='red')

legend('bottomright',c('Observed','Simulated df=1'),col=c('black','red'),pch=21)

cor(sim[,1],sim[,2],method='spearman')

copula_dist <- mvdc(copula=tCopula(rho,dim=2,df=8), margins=c("norm","norm"),

                    paramMargins=list(list(mean=cree_mu, sd=cree_sd),

                                      list(mean=yahoo_mu, sd=yahoo_sd)))

sim <- rmvdc(copula_dist, 3965)

plot(cree,yahoo,main='Returns')

points(sim[,1],sim[,2],col='red')

legend('bottomright',c('Observed','Simulated df=8'),col=c('black','red'),pch=21)

cor(sim[,1],sim[,2],method='spearman')

[1] 0.3929509

[1] 0.3911127

看到相关系数在第二个参数发生如此大的变化后，依然没有什么改变，而抽样结果显示：

自由的为1

自由度为8

一般来说，当df>10时候，就很接近高斯clopula的情况了

结语：

当它涉及到随机变量的联合行为建模时，copula是非常有用的工具，但是它也非常容易让你忽略建模时候的重点。在本例中，只使用了两个随机变量，在绘制模拟数据与真实数据的对比图时，一些拟合不足是显而易见的，但是当添加更多的变量时，结果将不再可用于绘制，并且留下相关系数和配对图，这可能会导致认为模拟值非常重要接近真实的时候，事实上你可能错过了很多。然而，copulas的真正作用是用边缘和依赖结构，推广具有许多随机变量的联合行为建模。 另一个经常被忽略的因素是依赖结构可能不是固定的，而是随时间而变化的。【时变copula，参考微博里的文章】。

无语。这些东西我以前就会了，又花时间搞了一遍。。没意思。