Using mlrv to anaylze data

Data analysis in the paper of Bai and Wu (2023b).

Loading data

Hong Kong circulatory and respiratory data.

 library(mlrv)
 library(foreach)
 library(magrittr)
 
 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`

Test for long memory

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)

Using the plug-in estimator for long-run covariance matrix function.

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.2838"
## [1] "RS"
## [1] "p-value 0.2896"
## [1] "VS"
## [1] "p-value 0.1148"
## [1] "KS"
## [1] "p-value 0.4124"

Debias difference-based estimator for long-run covariance matrix function.

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.6866"
## [1] "RS"
## [1] "p-value 0.8066"
## [1] "VS"
## [1] "p-value 0.5126"
## [1] "KS"
## [1] "p-value 0.8346"

Output

 rownames(pvmatrix) = c("plug","diff")
 colnames(pvmatrix) = c("KPSS","RS","VS","KS")
knitr::kable(pvmatrix,type="latex")

KPSS	RS	VS	KS
plug	0.2838	0.2896	0.1148	0.4124
diff	0.6866	0.8066	0.5126	0.8346

xtable::xtable(pvmatrix, digits = 3)

## % latex table generated in R 4.4.1 by xtable 1.8-4 package
## % Tue Jul 30 21:23:49 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
## \hline
## & KPSS & RS & VS & KS \\ 
## \hline
## plug & 0.284 & 0.290 & 0.115 & 0.412 \\ 
## diff & 0.687 & 0.807 & 0.513 & 0.835 \\ 
## \hline
## \end{tabular}
## \end{table}

Sensitivity Check

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.4256"
## [1] "RS"
## [1] "p-value 0.3638"
## [1] "VS"
## [1] "p-value 0.118"
## [1] "KS"
## [1] "p-value 0.561"

Shift = 1.2 with debias difference-based estimator. Choosing verbose_dist = TRUE can print more intermediate results, the exact values of smoothing parameters, the test statistic, the summary of bootstrap distribution and the p-value.

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 14 tau_n 0.332134206312301"
## [1] "test statistic: 141.654657280933"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max. 
## 13.81 107.33 215.16 368.80 460.94 5186.84 
## [1] "p-value 0.6462"
## [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. 
## 535.2 1023.8 1232.0 1302.6 1508.2 3383.7 
## [1] "p-value 0.6994"
## [1] "VS"
## [1] "gcv 0.193398841583897"
## [1] "m 17 tau_n 0.332134206312301"
## [1] "test statistic: 103.342038019402"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max. 
## 11.94 70.04 110.26 152.36 187.41 1230.74 
## [1] "p-value 0.538"
## [1] "KS"
## [1] "gcv 0.193398841583897"
## [1] "m 17 tau_n 0.382134206312301"
## [1] "test statistic: 671.676091515897"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max. 
## 352.6 720.4 923.3 1007.3 1213.0 3363.9 
## [1] "p-value 0.8108"

Results of shift = 1.2

 rownames(pvmatrix) = c("plug","diff")
 colnames(pvmatrix) = c("KPSS","RS","VS","KS")
knitr::kable(pvmatrix,type="latex")

KPSS	RS	VS	KS
plug	0.4256	0.3638	0.118	0.5610
diff	0.6462	0.6994	0.538	0.8108

xtable::xtable(pvmatrix, digits = 3)

## % latex table generated in R 4.4.1 by xtable 1.8-4 package
## % Tue Jul 30 21:25:13 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
## \hline
## & KPSS & RS & VS & KS \\ 
## \hline
## plug & 0.426 & 0.364 & 0.118 & 0.561 \\ 
## diff & 0.646 & 0.699 & 0.538 & 0.811 \\ 
## \hline
## \end{tabular}
## \end{table}

Shift = 0.8 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 = 0.8)
 print(paste("p-value",result_reg))
 pvmatrix[1,i] = result_reg
 i = i + 1
}

## [1] "KPSS"
## [1] "p-value 0.2458"
## [1] "RS"
## [1] "p-value 0.1662"
## [1] "VS"
## [1] "p-value 0.1234"
## [1] "KS"
## [1] "p-value 0.2726"

Shift = 0.8 with the debias difference-based estimator.

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 8 tau_n 0.382134206312301"
## [1] "test statistic: 166.543448031107"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max. 
## 12.78 99.10 200.38 331.58 421.82 3152.45 
## [1] "p-value 0.571"
## [1] "RS"
## [1] "gcv 0.128932561055931"
## [1] "m 18 tau_n 0.382134206312301"
## [1] "test statistic: 998.08124125936"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max. 
## 699.6 1276.2 1535.9 1614.1 1874.8 3730.1 
## [1] "p-value 0.9464"
## [1] "VS"
## [1] "gcv 0.128932561055931"
## [1] "m 9 tau_n 0.332134206312301"
## [1] "test statistic: 78.0587445148255"
## [1] "Bootstrap distribution"
## Min. 1st Qu. Median Mean 3rd Qu. Max. 
## 14.07 63.32 110.68 157.21 199.91 1498.85 
## [1] "p-value 0.6586"
## [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. 
## 290.6 695.8 910.2 986.8 1200.7 2953.0 
## [1] "p-value 0.7344"

Results of shift = 0.8

 rownames(pvmatrix) = c("plug","diff")
 colnames(pvmatrix) = c("KPSS","RS","VS","KS")
knitr::kable(pvmatrix,type="latex")

KPSS	RS	VS	KS
plug	0.2458	0.1662	0.1234	0.2726
diff	0.5710	0.9464	0.6586	0.7344

xtable::xtable(pvmatrix, digits = 3)

## % latex table generated in R 4.4.1 by xtable 1.8-4 package
## % Tue Jul 30 21:26:26 2024
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
## \hline
## & KPSS & RS & VS & KS \\ 
## \hline
## plug & 0.246 & 0.166 & 0.123 & 0.273 \\ 
## diff & 0.571 & 0.946 & 0.659 & 0.734 \\ 
## \hline
## \end{tabular}
## \end{table}

Test for structure stability

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 12 tau_n 0.364293094094381"
## [1] 10464.35
## V1 
## Min. : 1537 
## 1st Qu.: 3809 
## Median : 4838 
## Mean : 5239 
## 3rd Qu.: 6296 
## Max. :15552

p1

## [1] 0.0158