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.2966"
## [1] "RS"
## [1] "p-value 0.2508"
## [1] "VS"
## [1] "p-value 0.1548"
## [1] "KS"
## [1] "p-value 0.3528"

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.6974"
## [1] "RS"
## [1] "p-value 0.8246"
## [1] "VS"
## [1] "p-value 0.6196"
## [1] "KS"
## [1] "p-value 0.844"

Output

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

	KPSS	RS	VS	KS
plug	0.2966	0.2508	0.1548	0.3528
diff	0.6974	0.8246	0.6196	0.8440

xtable::xtable(pvmatrix, digits = 3)

## % latex table generated in R 4.6.1 by xtable 1.8-8 package
## % Sat Jul  4 19:06:06 2026
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
##   \hline
##  & KPSS & RS & VS & KS \\ 
##   \hline
## plug & 0.297 & 0.251 & 0.155 & 0.353 \\ 
##   diff & 0.697 & 0.825 & 0.620 & 0.844 \\ 
##    \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.4318"
## [1] "RS"
## [1] "p-value 0.3856"
## [1] "VS"
## [1] "p-value 0.1246"
## [1] "KS"
## [1] "p-value 0.5672"

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 16 tau_n 0.382134206312301"
## [1] "test statistic: 141.654657280932"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   19.47  112.25  225.82  393.07  497.31 5186.53 
## [1] "p-value 0.6678"
## [1] "RS"
## [1] "gcv 0.193398841583897"
## [1] "m 17 tau_n 0.332134206312301"
## [1] "test statistic: 1067.76713443354"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   555.6  1069.1  1273.8  1329.5  1525.3  3302.2 
## [1] "p-value 0.7518"
## [1] "VS"
## [1] "gcv 0.193398841583897"
## [1] "m 14 tau_n 0.332134206312301"
## [1] "test statistic: 103.342038019402"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.70   64.68  108.15  153.22  194.44 1381.83 
## [1] "p-value 0.5222"
## [1] "KS"
## [1] "gcv 0.193398841583897"
## [1] "m 15 tau_n 0.382134206312301"
## [1] "test statistic: 671.676091515896"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   307.1   686.3   914.5   989.2  1219.0  2983.4 
## [1] "p-value 0.7628"

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.4318	0.3856	0.1246	0.5672
diff	0.6678	0.7518	0.5222	0.7628

xtable::xtable(pvmatrix, digits = 3)

## % latex table generated in R 4.6.1 by xtable 1.8-8 package
## % Sat Jul  4 19:06:39 2026
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
##   \hline
##  & KPSS & RS & VS & KS \\ 
##   \hline
## plug & 0.432 & 0.386 & 0.125 & 0.567 \\ 
##   diff & 0.668 & 0.752 & 0.522 & 0.763 \\ 
##    \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.2402"
## [1] "RS"
## [1] "p-value 0.1752"
## [1] "VS"
## [1] "p-value 0.1026"
## [1] "KS"
## [1] "p-value 0.2512"

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 9 tau_n 0.382134206312301"
## [1] "test statistic: 166.543448031108"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   17.36  102.16  199.95  326.93  417.35 3380.05 
## [1] "p-value 0.565"
## [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. 
##   706.2  1267.1  1525.8  1596.8  1858.7  3958.0 
## [1] "p-value 0.9416"
## [1] "VS"
## [1] "gcv 0.128932561055931"
## [1] "m 14 tau_n 0.382134206312301"
## [1] "test statistic: 78.0587445148257"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.08   95.16  162.12  236.15  293.23 2290.23 
## [1] "p-value 0.8314"
## [1] "KS"
## [1] "gcv 0.128932561055931"
## [1] "m 18 tau_n 0.382134206312301"
## [1] "test statistic: 709.345279801765"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   390.8   851.2  1104.9  1206.1  1477.7  3394.8 
## [1] "p-value 0.881"

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.2402	0.1752	0.1026	0.2512
diff	0.5650	0.9416	0.8314	0.8810

xtable::xtable(pvmatrix, digits = 3)

## % latex table generated in R 4.6.1 by xtable 1.8-8 package
## % Sat Jul  4 19:07:09 2026
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
##   \hline
##  & KPSS & RS & VS & KS \\ 
##   \hline
## plug & 0.240 & 0.175 & 0.103 & 0.251 \\ 
##   diff & 0.565 & 0.942 & 0.831 & 0.881 \\ 
##    \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.414293094094381"
## [1] 10464.35
##        V1       
##  Min.   : 1769  
##  1st Qu.: 3824  
##  Median : 4836  
##  Mean   : 5176  
##  3rd Qu.: 6208  
##  Max.   :14062

p1

## [1] 0.0108