-
Notifications
You must be signed in to change notification settings - Fork 13
/
Copy pathte1.r
101 lines (96 loc) · 3.11 KB
/
te1.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
library(copent) # transfer entropy (TE)
library(CondIndTests) # kernel-based CI test (KCI)
library(RCIT) # Randomized conditional Correlation Test (RCoT)
library(cdcsis) # conditional distance correlation (CDC)
library(FOCI) # conditional dependence coefficient (CODEC)
library(GeneralisedCovarianceMeasure) # Generalised Covariance Measure (GCM)
library(weightedGCM) # weighted GCM
library(KPC) # Kernel Partial Correlation (KPC)
library(ppcor) # Partial Correlation (pcor)
library(CondCopulas) # Conditional Kendall's Tau (CKT)
library(EDMeasure) # Conditional Mean Dependence (CMD)
source("https://raw.githubusercontent.com/lassepetersen/partial-copula-CI-test/main/parCopCITest.R") # partial copula based CI test (pcop)
prsa2010data = read.csv("https://archive.ics.uci.edu/ml/machine-learning-databases/00381/PRSA_data_2010.1.1-2014.12.31.csv")
# id: 6(PM2.5), 7(Dew Point), 8(Temperature), 9(Pressure), 11(Cumulative Wind Speed)
idx = c(6,9)
data = prsa2010data[2200:2700, idx]
te1 = 0
kci1 = 0
rcot1 = 0
cdc1 = 0
codec1 = 0
gcm1 = 0
wgcm1 = 0
kpc1 = 0
pcor1 = 0
ckt1 = 0
cmd1 = 0
pcop1 = 0
for (lag in 1:24){
pm25a = data[1:(501-lag),1]
pm25b = data[(lag+1):501,1]
v1 = data[1:(501-lag),2]
# conditional independence test with copula entropy
te1[lag] = ci(pm25b,v1,pm25a)
# or estimating transfer entropy directly as below
# te1[lag] = transent(data[,1],data[,2],lag)
kci1[lag] = KCI(pm25b,v1,pm25a)$testStatistic
rcot1[lag] = RCIT(pm25b,v1,pm25a)$Sta
cdc1[lag] = cdcor(pm25b,v1,pm25a)$statistic
codec1[lag] = codec(pm25b,v1,pm25a)
gcm1[lag] = gcm.test(pm25b,v1,pm25a)$test.statistic
wgcm1[lag] = wgcm.est(pm25b,v1,pm25a, regr.meth = "xgboost")
kpc1[lag] = KPCRKHS(pm25b,pm25a,v1)
pcor1[lag] = pcor.test(pm25b,v1,pm25a)$statistic
ckt1[lag] = CKT.estimate(pm25b,v1,pm25a, methodEstimation = "tree", h = 0.1)
cmd1[lag] = cmdm_test(pm25b,v1,pm25a)$stat
pcop1[lag] = test_CI(pm25b,v1,pm25a)$statistic
}
# TE via CE
x11()
plot(te1, xlab = "lag (hours)", ylab = "Transfer Entropy", main = "Pressure")
lines(te1)
# KCI
x11()
plot(kci1, xlab = "lag (hours)", ylab = "KCI", main = "Pressure")
lines(kci1)
# RCoT
x11()
plot(rcot1, xlab = "lag (hours)", ylab = "RCoT", main = "Pressure")
lines(rcot1)
# CDC
x11()
plot(cdc1, xlab = "lag (hours)", ylab = "CDC", main = "Pressure")
lines(cdc1)
# CODEC
x11()
plot(codec1, xlab = "lag (hours)", ylab = "CODEC", main = "Pressure")
lines(codec1)
# GCM
x11()
plot(gcm1, xlab = "lag (hours)", ylab = "GCM", main = "Pressure")
lines(gcm1)
# wGCM
x11()
plot(wgcm1, xlab = "lag (hours)", ylab = "wGCM", main = "Pressure")
lines(wgcm1)
# KPC
x11()
plot(kpc1, xlab = "lag (hours)", ylab = "KPC", main = "Pressure")
lines(kpc1)
# Partial Correlation
x11()
plot(pcor1, xlab = "lag (hours)", ylab = "Partial Correlation", main = "Pressure")
lines(pcor1)
# Conditional Kendall's Tau
x11()
plot(ckt1, xlab = "lag (hours)", ylab = "CKT", main = "Pressure")
lines(ckt1)
# Conditional Mean Dependence
x11()
plot(cmd1, xlab = "lag (hours)", ylab = "CMD", main = "Pressure")
lines(cmd1)
# partial copula based CI test
x11()
plot(pcop1, xlab = "lag (hours)", ylab = "PartialCopula", main = "Pressure")
lines(pcop1)