-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathRestricted_t_FGM_and_CA_n=200.R
101 lines (94 loc) · 4.53 KB
/
Restricted_t_FGM_and_CA_n=200.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
rm(list=ls(all=TRUE))
library(xtable)
#########################################################################
# R functions for testing independence versus positive quadrant #
# dependence corresponding to the manuscript titled, #
# "Testing for positive quadrant dependence." #
# Date: 06/17/2018 #
#########################################################################
library(Rcpp)
library(copula)
library(FDGcopulas)
library(GFGM.copula)
source("https://raw.githubusercontent.com/cftang9/PQD/master/EL_PQD_Library.R")
#########################################################################
# Generating critical values for EL, KS, CvM, AD, Spearman, and Kendall
# tests for independence versus PQD from uniform (0,1)^2 with sample
# size n and 10000 Monte Carlo replications.
# (It may take some time to generate critical for large sample sizes,
# for example, 100.)
n = 200; #sample size;
alpha = 0.05; #significance level
#CV = c(1.425291, 0.6719852, 0.08987011, 4.724722, 0.2342857, 0.1591837) #n=50 set.seed(50)
#CV = c(1.402599, 0.6740648, 0.07325512, 3.728709, 0.1646115, 0.1119192) #n=100 set.seed(100)
CV = c(1.393776, 0.6801778, 0.06367364, 2.948650, 0.1178105, 0.07880402) #n=200 set.seed(200)
#########################################################################
# Perform simulation using conditional 2D Student's t distribution on the first quadrant
# with sample size n and B Monte Carlo replications.
B = 10000#number of Monte Carlo replications
dt = c(1,2,3,4,5,10,50,100); ndt = length(dt);
PowerT = array(0,c(6,ndt));
set.seed(1600)
for(i in 1:ndt){
myCop.t = tCopula(0,dim=2,df=dt[i]);
for(b in 1:B){
X = rCopula(n, myCop.t); X = FlipQ(X);
PowerT[1,i] = c(Rcpp_IP_FF(X,n)>CV[1])/B + PowerT[1,i];
temp = Rcpp_KCA_01(X,n);
PowerT[2,i] = c(temp[1]>CV[2])/B + PowerT[2,i];
PowerT[3,i] = c(temp[2]>CV[3])/B + PowerT[3,i];
PowerT[4,i] = c(temp[3]>CV[4])/B + PowerT[4,i];
PowerT[5,i] = c(cor(X,method="spearman")[2]>CV[5])/B + PowerT[5,i]
PowerT[6,i] = c(cor(X,method="kendall")[2]>CV[6])/B + PowerT[6,i]
#PowerT[7,i] = c(Vexler2014(X[,1], X[,2]) > CV_Vexler2014(n))/B + PowerT[7,i]
}
}
colnames(PowerT) = dt; rownames(PowerT) = c("EL","KS","CvM","AD","spearman_rho","kendall_tau")
#########################################################################
# Perform simulation using conditional 2D FGM copula on the first quadrant
# with sample size n and B Monte Carlo replications.
B = 10000#number of Monte Carlo replications
dt = seq(0,1,by=0.2); ndt = length(dt);
PowerFGM = array(0,c(6,ndt));
set.seed(1700)
for(i in 1:ndt){
fgm.cop <- fgmCopula(param=dt[i],dim=2)
for(b in 1:B){
X <- rCopula(n, fgm.cop)
PowerFGM[1,i] = c(Rcpp_IP_FF(X,n)>CV[1])/B + PowerFGM[1,i];
temp = Rcpp_KCA_01(X,n);
PowerFGM[2,i] = c(temp[1]>CV[2])/B + PowerFGM[2,i];
PowerFGM[3,i] = c(temp[2]>CV[3])/B + PowerFGM[3,i];
PowerFGM[4,i] = c(temp[3]>CV[4])/B + PowerFGM[4,i];
PowerFGM[5,i] = c(cor(X,method="spearman")[2]>CV[5])/B + PowerFGM[5,i]
PowerFGM[6,i] = c(cor(X,method="kendall")[2]>CV[6])/B + PowerFGM[6,i]
#PowerFGM[7,i] = c(Vexler2014(X[,1], X[,2]) > CV_Vexler2014(n))/B + PowerFGM[7,i]
}
}
colnames(PowerFGM) = dt; rownames(PowerFGM) = c("EL","KS","CvM","AD","spearman_rho","kendall_tau")
#########################################################################
# Perform simulation using Cuadras-Aug'e copula distribution on the first quadrant
# with sample size n and B Monte Carlo replications.
B = 10000#number of Monte Carlo replications
dt = seq(0,1,by=0.2); ndt = length(dt);
PowerCA = array(0,c(6,ndt));
set.seed(1800)
for(i in 1:ndt){
theta = c(dt[i], dt[i])
myFDGcopula <- FDGcopula("cuadrasauge", theta)
for(b in 1:B){
X = rFDG(n, myFDGcopula)
PowerCA[1,i] = c(Rcpp_IP_FF(X,n)>CV[1])/B + PowerCA[1,i];
temp = Rcpp_KCA_01(X,n);
PowerCA[2,i] = c(temp[1]>CV[2])/B + PowerCA[2,i];
PowerCA[3,i] = c(temp[2]>CV[3])/B + PowerCA[3,i];
PowerCA[4,i] = c(temp[3]>CV[4])/B + PowerCA[4,i];
PowerCA[5,i] = c(cor(X,method="spearman")[2]>CV[5])/B + PowerCA[5,i]
PowerCA[6,i] = c(cor(X,method="kendall")[2]>CV[6])/B + PowerCA[6,i]
#PowerCA[7,i] = c(Vexler2014(X[,1], X[,2]) > CV_Vexler2014(n))/B + PowerCA[7,i]
}
}
colnames(PowerCA) = dt; rownames(PowerCA) = c("EL","KS","CvM","AD","spearman_rho","kendall_tau")
print("Restricted Students' t"); print(PowerT,digits=3)
print("FGM"); print(PowerFGM,digits=3)
print("CA"); print(PowerCA,digits=3)