Data analysis in the paper of Bai and Wu (2023b).
Hong Kong circulatory and respiratory data.
library(mlrv)
library(foreach)
library(magrittr)
load("../data/hk_data.RData")
# data(hk_data)
colnames(hk_data) = c("SO2","NO2","Dust","Ozone","Temperature",
"Humidity","num_circu","num_respir","Hospital Admission",
"w1","w2","w3","w4","w5","w6")
n = nrow(hk_data)
t = (1:n)/n
hk = list()
hk$x = as.matrix(cbind(rep(1,n), scale(hk_data[,1:3])))
hk$y = hk_data$`Hospital Admission`
pvmatrix = matrix(nrow=2, ncol=4)
###inistialization
setting = list(B = 5000, gcv = 1, neighbour = 1)
setting$lb = floor(10/7*n^(4/15)) - setting$neighbour
setting$ub = max(floor(25/7*n^(4/15))+ setting$neighbour,
setting$lb+2*setting$neighbour+1)
setting$lrvmethod =0.
i=1
# print(rule_of_thumb(y= hk$y, x = hk$x))
for(type in c("KPSS","RS","VS","KS")){
setting$type = type
print(type)
result_reg = heter_covariate(list(y= hk$y, x = hk$x), setting, mvselect = -2)
print(paste("p-value",result_reg))
pvmatrix[1,i] = result_reg
i = i + 1
}
## [1] "KPSS"
## [1] "p-value 0.2862"
## [1] "RS"
## [1] "p-value 0.3022"
## [1] "VS"
## [1] "p-value 0.126"
## [1] "KS"
## [1] "p-value 0.4046"
setting$lrvmethod =1
i=1
for(type in c("KPSS","RS","VS","KS"))
{
setting$type = type
print(type)
result_reg = heter_covariate(list(y= hk$y, x = hk$x), setting, mvselect = -2)
print(paste("p-value",result_reg))
pvmatrix[2,i] = result_reg
i = i + 1
}
## [1] "KPSS"
## [1] "p-value 0.7034"
## [1] "RS"
## [1] "p-value 0.8002"
## [1] "VS"
## [1] "p-value 0.7702"
## [1] "KS"
## [1] "p-value 0.8612"
rownames(pvmatrix) = c("plug","diff")
colnames(pvmatrix) = c("KPSS","RS","VS","KS")
knitr::kable(pvmatrix,type="latex")
KPSS | RS | VS | KS | |
---|---|---|---|---|
plug | 0.2862 | 0.3022 | 0.1260 | 0.4046 |
diff | 0.7034 | 0.8002 | 0.7702 | 0.8612 |
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Tue Dec 26 10:15:14 2023
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
## \hline
## & KPSS & RS & VS & KS \\
## \hline
## plug & 0.286 & 0.302 & 0.126 & 0.405 \\
## diff & 0.703 & 0.800 & 0.770 & 0.861 \\
## \hline
## \end{tabular}
## \end{table}
Using parameter `shift’ to multiply the GCV selected bandwidth by a factor. - Shift = 1.2 with plug-in estimator.
pvmatrix = matrix(nrow=2, ncol=4)
setting$lrvmethod = 0
i=1
for(type in c("KPSS","RS","VS","KS")){
setting$type = type
print(type)
result_reg = heter_covariate(list(y= hk$y, x = hk$x),
setting,
mvselect = -2, shift = 1.2)
print(paste("p-value",result_reg))
pvmatrix[1,i] = result_reg
i = i + 1
}
## [1] "KPSS"
## [1] "p-value 0.418"
## [1] "RS"
## [1] "p-value 0.3628"
## [1] "VS"
## [1] "p-value 0.1198"
## [1] "KS"
## [1] "p-value 0.569"
setting$lrvmethod =1
i=1
for(type in c("KPSS","RS","VS","KS"))
{
setting$type = type
print(type)
result_reg = heter_covariate(list(y= hk$y, x = hk$x),
setting,
mvselect = -2, verbose_dist = TRUE, shift = 1.2)
print(paste("p-value",result_reg))
pvmatrix[2,i] = result_reg
i = i + 1
}
## [1] "KPSS"
## [1] "gcv 0.193398841583897"
## [1] "m 8 tau_n 0.332134206312301"
## [1] "test statistic: 141.654657280933"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 11.78 69.74 142.39 235.86 302.58 3490.68
## [1] "p-value 0.5018"
## [1] "RS"
## [1] "gcv 0.193398841583897"
## [1] "m 15 tau_n 0.332134206312301"
## [1] "test statistic: 1067.76713443354"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 494.1 1027.5 1232.8 1304.9 1512.1 3531.3
## [1] "p-value 0.7044"
## [1] "VS"
## [1] "gcv 0.193398841583897"
## [1] "m 8 tau_n 0.332134206312301"
## [1] "test statistic: 103.342038019402"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 6.723 41.592 71.645 104.158 129.037 994.853
## [1] "p-value 0.336"
## [1] "KS"
## [1] "gcv 0.193398841583897"
## [1] "m 17 tau_n 0.332134206312301"
## [1] "test statistic: 671.676091515897"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 337.2 719.1 920.8 1003.3 1208.8 3178.5
## [1] "p-value 0.812"
rownames(pvmatrix) = c("plug","diff")
colnames(pvmatrix) = c("KPSS","RS","VS","KS")
knitr::kable(pvmatrix,type="latex")
KPSS | RS | VS | KS | |
---|---|---|---|---|
plug | 0.4180 | 0.3628 | 0.1198 | 0.569 |
diff | 0.5018 | 0.7044 | 0.3360 | 0.812 |
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Tue Dec 26 10:16:38 2023
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
## \hline
## & KPSS & RS & VS & KS \\
## \hline
## plug & 0.418 & 0.363 & 0.120 & 0.569 \\
## diff & 0.502 & 0.704 & 0.336 & 0.812 \\
## \hline
## \end{tabular}
## \end{table}
pvmatrix = matrix(nrow=2, ncol=4)
setting$lrvmethod =0
i=1
for(type in c("KPSS","RS","VS","KS")){
setting$type = type
print(type)
result_reg = heter_covariate(list(y= hk$y, x = hk$x),
setting,
mvselect = -2, shift = 0.8)
print(paste("p-value",result_reg))
pvmatrix[1,i] = result_reg
i = i + 1
}
## [1] "KPSS"
## [1] "p-value 0.1714"
## [1] "RS"
## [1] "p-value 0.123"
## [1] "VS"
## [1] "p-value 0.1156"
## [1] "KS"
## [1] "p-value 0.2482"
setting$lrvmethod =1
i=1
for(type in c("KPSS","RS","VS","KS"))
{
setting$type = type
print(type)
result_reg = heter_covariate(list(y= hk$y, x = hk$x),
setting,
mvselect = -2, verbose_dist = TRUE, shift = 0.8)
print(paste("p-value",result_reg))
pvmatrix[2,i] = result_reg
i = i + 1
}
## [1] "KPSS"
## [1] "gcv 0.128932561055931"
## [1] "m 9 tau_n 0.382134206312301"
## [1] "test statistic: 166.543448031107"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18.63 101.76 202.70 335.57 421.49 3033.43
## [1] "p-value 0.5762"
## [1] "RS"
## [1] "gcv 0.128932561055931"
## [1] "m 17 tau_n 0.382134206312301"
## [1] "test statistic: 998.08124125936"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 652.5 1239.0 1498.5 1569.0 1826.4 3783.7
## [1] "p-value 0.9314"
## [1] "VS"
## [1] "gcv 0.128932561055931"
## [1] "m 9 tau_n 0.382134206312301"
## [1] "test statistic: 78.0587445148255"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 12.79 65.05 109.38 154.78 193.22 1514.37
## [1] "p-value 0.6688"
## [1] "KS"
## [1] "gcv 0.128932561055931"
## [1] "m 9 tau_n 0.332134206312301"
## [1] "test statistic: 709.345279801765"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 288.9 695.7 912.5 991.8 1194.7 2857.1
## [1] "p-value 0.7358"
rownames(pvmatrix) = c("plug","diff")
colnames(pvmatrix) = c("KPSS","RS","VS","KS")
knitr::kable(pvmatrix,type="latex")
KPSS | RS | VS | KS | |
---|---|---|---|---|
plug | 0.1714 | 0.1230 | 0.1156 | 0.2482 |
diff | 0.5762 | 0.9314 | 0.6688 | 0.7358 |
## % latex table generated in R 4.3.2 by xtable 1.8-4 package
## % Tue Dec 26 10:17:49 2023
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
## \hline
## & KPSS & RS & VS & KS \\
## \hline
## plug & 0.171 & 0.123 & 0.116 & 0.248 \\
## diff & 0.576 & 0.931 & 0.669 & 0.736 \\
## \hline
## \end{tabular}
## \end{table}
Test if the coefficient function of “SO2”,“NO2”,“Dust” of the second year is constant.
hk$x = as.matrix(cbind(rep(1,n), (hk_data[,1:3])))
hk$y = hk_data$`Hospital Admission`
setting$type = 0
setting$bw_set = c(0.1, 0.35)
setting$eta = 0.2
setting$lrvmethod = 1
setting$lb = 10
setting$ub = 15
hk1 = list()
hk1$x = hk$x[366:730,]
hk1$y = hk$y[366:730]
p1 <- heter_gradient(hk1, setting, mvselect = -2, verbose = T)
## [1] "m 11 tau_n 0.414293094094381"
## [1] 10464.35
## V1
## Min. : 1343
## 1st Qu.: 3343
## Median : 4328
## Mean : 4674
## 3rd Qu.: 5633
## Max. :13336
## [1] 0.0058