前言：本咩還曾經翻譯過兩篇copula相關的基礎文章，還有一個時變copula的代碼分享，都在我的微博上有。但是那個微博我現在不怎么用了。感興趣的同學可搜搜微博名：補補補牢

今天分享的這些代碼，原文請見：modelling-dependence-with-copulas-in-r，本文僅供本咩學習用途。

我不會被打倒

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

給定相關矩陣，從多元正態(tài)分布中抽樣，利用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')

這就是生成的一個簡單的依賴結構

通過上面的多元正態(tài)樣本，我們建立了一個具有預期依賴結構的樣本，但是卻具有任意的邊緣分布。

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強調極端結果：它通常適用于極端值（分布的尾部）高度相關的現象建模。?

三、對邊緣分布建模

假設兩個收益率序列服從正態(tài)分布

cree_mu <- mean(cree)

cree_sd <- sd(cree)

yahoo_mu <- mean(yahoo)

yahoo_sd <- sd(yahoo)

#看下實際數據和正態(tài)分布的情況

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)

文中為了簡化，用的正態(tài)，實際上有高峰特征

#構造函數并抽樣
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的模擬值將會發(fā)生怎樣的變化

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

看到相關系數在第二個參數發(fā)生如此大的變化后，依然沒有什么改變，而抽樣結果顯示：

自由的為1

自由度為8

一般來說，當df>10時候，就很接近高斯clopula的情況了

結語：

當它涉及到隨機變量的聯合行為建模時，copula是非常有用的工具，但是它也非常容易讓你忽略建模時候的重點。在本例中，只使用了兩個隨機變量，在繪制模擬數據與真實數據的對比圖時，一些擬合不足是顯而易見的，但是當添加更多的變量時，結果將不再可用于繪制，并且留下相關系數和配對圖，這可能會導致認為模擬值非常重要接近真實的時候，事實上你可能錯過了很多。然而，copulas的真正作用是用邊緣和依賴結構，推廣具有許多隨機變量的聯合行為建模。 另一個經常被忽略的因素是依賴結構可能不是固定的，而是隨時間而變化的。【時變copula，參考微博里的文章】。

無語。這些東西我以前就會了，又花時間搞了一遍。。沒意思。