附上论文的R代码，供参考。
------
library(mnormt)
library(copent)

## for section{Implementation}
data("airquality")
x1 = airquality[,1:4]
xc1 = construct_empirical_copula(x1)

data = list(x = x1, xc = xc1)
dat <- as.data.frame(do.call(rbind,data))
dat$type <- rep(c('x','xc'),each = 153)

foo.upper <- function(x,y,...){
  points(x[dat\$type == 'x'],y[dat\$type == 'x'],col = 'red',...)
}
foo.lower <- function(x,y,...){
  plot.window(xlim = c(0,1), ylim = c(0,1))
  points(x[dat\$type == 'xc'],y[dat\$type == 'xc'],col = 'blue', ...)
}

pairs(dat[,-5],lower.panel = foo.lower,upper.panel = foo.upper)

entknn(xc1)
copent(x1)

## for Section{Two Examples}
## Simulated Example
rho = 0.75
sigma = matrix(c(1,rho,rho,1),2,2)
x = rmnorm(500,c(0,0),sigma)
truevalue = -0.5 * log(1- rho^2)
est = copent(x)

# Example on Causal Discovery
# Data URL: http://archive.ics.uci.edu/ml/datasets/Beijing+PM2.5+Data
prsa2010data = read.csv('~/Rworks/beijingair/PRSA2010.csv')
data = prsa2010data[2200:2700,c(6,9)]
tslag = 0
for (lag in 1:24){
        pm25a = data[1:(501-lag),1]
        pm25b = data[(lag+1):501,1]
        var1 = data[1:(501-lag),2]
        data1 = cbind(pm25a, pm25b, var1)
        tslag[lag] = copent(data1) - copent(data1[,c(1,2)]) - copent(data1[,c(1,3)])
}
plot(tslag, xlab = "lag (hours)", ylab = "Transfer Entropy", main = "Pressure")
lines(tslag)

怎么删自己重复发的贴？求版主帮忙。