基于copula熵的时序因果发现及对北京空气污染问题的分析代码示例

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

因果关系发现在科学研究中非常重要。传递熵（Transfer Entropy）是广泛采用的因果关系度量，基于条件独立的概念，较之其他经验式因果关系建模方法更科学合理。这里分享一个利用copula熵从数据中估计因果关系的例子，利用copula熵估计传递熵，具体见如下文献：

Jian Ma. Estimating Transfer Entropy via Copula Entropy. arXiv:1910.04375, 2019.
     URL: https://arxiv.org/abs/1910.04375

文中给出了利用copula熵估计传递熵的非参数方法，方法简便易行，普遍适用。下面的代码例子利用这个方法做北京地区气象因素对PM2.5的因果关系估计的分析。其中，第11行代码（3个copula熵的运算）即为估计传递熵的方法。代码在github的共享仓库网址为:     https://github.com/majianthu/transferentropy

1. library(copent)
   
2. 
   
3. prsa2010data = read.csv("https://archive.ics.uci.edu/ml/machine-learning-databases/00381/PRSA_data_2010.1.1-2014.12.31.csv")
   
4. data = prsa2010data[2200:2700,c(6,9)]
   
5. tslag = 0
   
6. for (lag in 1:24){
   
7.   pm25a = data[1:(501-lag),1]
   
8.   pm25b = data[(lag+1):501,1]
   
9.   v1 = data[1:(501-lag),2]
   
10.   data1 = cbind(pm25a, pm25b, v1)
    
11.   tslag[lag] = copent(data1) - copent(data1[,c(1,2)]) - copent(data1[,c(1,3)])
    
12. }
    
13. x11()
    
14. plot(tslag, xlab = "lag (hours)", ylab = "Transfer Entropy", main = "Pressure")
    
15. lines(tslag)

复制代码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏7 回帖

关键词：Copula opula 因果关系 空气污染 北京地区 因果发现 大气污染 传递熵 copula熵 因果发现 大气污染 传递熵 copula熵 因果发现 大气污染 传递熵 copula熵 因果发现 大气污染 传递熵 copula熵 因果发现 大气污染 传递熵 copula熵

majianthu 发表于 2020-8-27 09:27
因果关系发现在科学研究中非常重要。传递熵（Transfer Entropy）是广泛采用的因果关系度量，基于条件独立的 ...

谢谢楼主分享

适合时间序列数据？

楼主是马健老师本尊吗？能介绍一下详细的如何利用COPULA熵相关，因果分析的实操步骤吗？

lonestone 发表于 2020-8-28 13:31
楼主是马健老师本尊吗？能介绍一下详细的如何利用COPULA熵相关，因果分析的实操步骤吗？

就是将传递熵转换成3个copula熵的运算，上面的代码（第11行）已经演示的很明白了。运行这个代码需要安装copent包，其中的copent函数是用于估计copula熵的。如对更多相关理论感兴趣，可阅读论文。

majianthu 发表于 2020-8-27 09:27
因果关系发现在科学研究中非常重要。传递熵（Transfer Entropy）是广泛采用的因果关系度量，基于条件独立的 ...

收到，谢谢

论文代码已在github共享，网址：
    https://github.com/majianthu/transferentropy
可参考示例做时序因果分析。

github代码更新：增加了python版本的示例代码   网址：https://github.com/majianthu/transferentropy

代码如下：

1. from copent import copent
   
2. from pandas import read_csv
   
3. import numpy as np
   
4. import matplotlib.pyplot as plt
   
5. 
   
6. url = "https://archive.ics.uci.edu/ml/machine-learning-databases/00381/PRSA_data_2010.1.1-2014.12.31.csv"
   
7. prsa2010 = read_csv(url)
   
8. # index: 5(PM2.5),6(Dew Point),7(Temperature),8(Pressure),10(Cumulative Wind Speed)
   
9. data = prsa2010.iloc[2200:2700,[5,8]].values
   
10. 
    
11. te = np.zeros(24)
    
12. for lag in range(1,25):
    
13.         pm25a = data[0:(500-lag),0]
    
14.         pm25b = data[lag:500,0]
    
15.         cause = data[0:(500-lag),1]
    
16.         data123 = np.vstack((pm25a,pm25b,cause)).T
    
17.         data12 = np.vstack((pm25a,pm25b)).T
    
18.         data13 = np.vstack((pm25a,cause)).T
    
19.         # estimating transfer entropy
    
20.         te[lag-1] = copent(data123) - copent(data12) - copent(data13)
    
21.         str = "TE from pressure to PM2.5 at %d hours lag : %f" %(lag,te[lag-1])
    
22.         print(str)
    
23.         
    
24. plt.plot(te)
    
25. plt.title("Pressure on PM2.5")
    
26. plt.xlabel('lag(hours)')
    
27. pos = np.array([4,9,14,19])
    
28. plt.xticks(pos,pos+1)
    
29. plt.ylabel('Transfer Entropy')
    
30. plt.show()

复制代码

示例代码更新，加入了与Kernel-based 条件独立测试和中山大学开发的conditional distance correlation的对比。采用的R包是CondIndTests和cdcsis。

1. library(copent) # copula entropy
   
2. library(CondIndTests) # kernel-based CI test
   
3. library(cdcsis) # conditional distance correlation
   
4. 
   
5. url = "https://archive.ics.uci.edu/ml/machine-learning-databases/00381/PRSA_data_2010.1.1-2014.12.31.csv"
   
6. prsa2010data = read.csv(url)
   
7. # id: 6(PM2.5), 7(Dew Point), 8(Temperature), 9(Pressure), 11(Cumulative Wind Speed)
   
8. idx = c(6,9)
   
9. data = prsa2010data[2200:2700, idx]
   
10. 
    
11. te1 = 0
    
12. kci1 = 0
    
13. cdc1 = 0
    
14. for (lag in 1:24){
    
15.   pm25a = data[1:(501-lag),1]
    
16.   pm25b = data[(lag+1):501,1]
    
17.   v1 = data[1:(501-lag),2]
    
18.   data1 = cbind(pm25a, pm25b, v1)
    
19.   # estimating transfer entropy via copula entropy
    
20.   te1[lag] = copent(data1) - copent(data1[,c(1,2)]) - copent(data1[,c(1,3)])
    
21.   kci1[lag] = KCI(pm25b,v1,pm25a)$testStatistic
    
22.   cdc1[lag] = cdcor(pm25b,v1,pm25a)
    
23. }
    
24. 
    
25. x11()
    
26. plot(te1, xlab = "lag (hours)", ylab = "Transfer Entropy", main = "Pressure")
    
27. lines(te1)
    
28. x11()
    
29. plot(kci1, xlab = "lag (hours)", ylab = "KCI", main = "Pressure")
    
30. lines(kci1)
    
31. x11()
    
32. cdc1 = unlist(cdc1)
    
33. plot(cdc1, xlab = "lag (hours)", ylab = "CDC", main = "Pressure")
    
34. lines(cdc1)

复制代码

马老师怎么不写 vignetti 文件？  PDF 格式的比较难写， 但是 .html 格式的非常容易。也可以写出数学公式。

论坛支持 LATEX 的。 我给马老师回帖发几个数学公式，这里能支持的，  Rstudio 也支持。