本文原版地址：
http://firsttimeprogrammer.blogspot.ca/2015/02/how-to-fit-copula-model-in-r.html
本文是一系列文章中的一部分，翻譯主要是為了自己理解copula是什么，不精確到句子和詞纸镊。
我從事copula的相關(guān)工作很久了,老實說,R關(guān)于此問題的文檔,比起那些python來說,差得太遠(yuǎn)了.由于copula本身就很難,這就使得問題更難被解決了.其中最大的問題還是在于例子少,用戶也少.如果你有好資源推薦,歡迎在下面留言,我也會留下一部分資源.
如果你剛聽說copulas,或許可以先去看看introduction to the Gumbel copula in R here.
我接下來要用的package是copula package,一個很好的工具,安裝也很便利.
Dataset
這里我用一個簡單數(shù)據(jù),股票x和股票y,可以在此下載.
首先我們要加載數(shù)據(jù)到矩陣.記得要先加載copula package.
x <- read.table("/home/kevin/Downloads/x.txt")
y <- read.table("/home/kevin/Downloads/y.txt")
mat <- matrix(nrow = 100, ncol = 2)
for(i in 1:100)
{
mat[i,1] <- x[,1][i]
mat[i,2] <- y[,1][i]
}
plot(mat[,1],mat[,2],main="Returns",xlab="x",ylab="y",col="blue")
# cor(mat[,1],mat[,2])
用plot函數(shù)作圖之后如下.
下一步就是要進(jìn)行copula了.為了讓數(shù)據(jù)可用,我們需要選擇合適的copula模型,該模型應(yīng)該基于給出的數(shù)據(jù)及一些其他因素.我們可以初步估計出這個數(shù)據(jù)集有一定的正相關(guān),因此能夠反映這種關(guān)系的copula都是可以的.要小心的是,copula模型會越來越復(fù)雜,這種靠眼睛看的方式很難一直奏效.我選擇用normal copula.當(dāng)然下列過程對其他copula也適用.
我們來擬合以下數(shù)據(jù).
# normal copula
normal.cop <- normalCopula(dim = 2)
fit.cop <- fitCopula(normal.cop,pobs(mat),method="ml")
# Coefficients
rho <- coef(fit.cop)
print(rho)
注意數(shù)據(jù)要通過pobs()輸送朵逝，該參數(shù)會將真實數(shù)據(jù)轉(zhuǎn)換為偽數(shù)據(jù)以符合[0,1]區(qū)間要求.
另外,這里我們選用的是ml方法－－極大似然方法踩寇，其他方法也可以唯蝶，比如itau.
我們可以plot偽數(shù)據(jù)和模擬數(shù)據(jù),可以發(fā)現(xiàn)copula模擬很好的匹配了偽數(shù)據(jù).
# Pseudo observations
p_obs <- pobs(mat)
plot(p_obs[,1],p_obs[,2],main="Pseudo/simulated observations:
BLUE/RED", xlab="u",ylab="v",col="blue")
# Simulate data
set.seed(100)
u1 = rCopula(500,normalCopula(coef(fit.cop),dim = 2))
points(u1[,1],u1[,2],col="red")
我們選用的這個copula還可以,但不是最好的,因為它有較大的尾部相關(guān),而數(shù)據(jù)本身并沒有那么大.不過這只是我們的一個嘗試而已.
我們還可以選擇用柱狀圖顯示隨機(jī)變量的分布,如下:
# Plot data with histograms
xhist <- hist(mat[,1],breaks = 30,plot = FALSE)
yhist <- hist(mat[,2],breaks = 30,plot = FALSE)
top <- max(c(xhist$counts, yhist$counts))
xrange <- c(-4,4)
yrange <- c(-6,6)
nf <- layout(matrix(c(2,0,1,3),2,2,byrow=TRUE),c(3,1),c(1,3),
TRUE)
# for the matrix, it simplely means put position of #n of plot;
# for c(3,1), it means col1 is 3x width more than col2, vice versa
par(mar=c(3,3,1,1))
plot(mat[,1],mat[,2],xlim=xrange,ylim=yrange,xlab="",ylab="")
par(mar=c(0,3,1,1))
barplot(xhist$counts, axes = FALSE, ylim = c(0,top),space=0)
par(mar=c(3,0,1,1))
barplot(yhist$counts, axes = FALSE, xlim = c(0,top),space=0,
horiz = TRUE)
得到如下的圖像:
