Using mlrv to anaylze data

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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.2886"
## [1] "RS"
## [1] "p-value 0.2792"
## [1] "VS"
## [1] "p-value 0.1458"
## [1] "KS"
## [1] "p-value 0.4296"

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.516"
## [1] "RS"
## [1] "p-value 0.8568"
## [1] "VS"
## [1] "p-value 0.5028"
## [1] "KS"
## [1] "p-value 0.8134"

Output

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

	KPSS	RS	VS	KS
plug	0.2886	0.2792	0.1458	0.4296
diff	0.5160	0.8568	0.5028	0.8134

xtable::xtable(pvmatrix, digits = 3)

## % latex table generated in R 4.6.0 by xtable 1.8-8 package
## % Thu Apr  2 16:13:34 2026
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
##   \hline
##  & KPSS & RS & VS & KS \\ 
##   \hline
## plug & 0.289 & 0.279 & 0.146 & 0.430 \\ 
##   diff & 0.516 & 0.857 & 0.503 & 0.813 \\ 
##    \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.4406"
## [1] "RS"
## [1] "p-value 0.3938"
## [1] "VS"
## [1] "p-value 0.1212"
## [1] "KS"
## [1] "p-value 0.5818"

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 8 tau_n 0.382134206312301"
## [1] "test statistic: 141.654657280933"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   11.29   71.65  135.17  221.78  280.41 2459.83 
## [1] "p-value 0.483"
## [1] "RS"
## [1] "gcv 0.193398841583897"
## [1] "m 17 tau_n 0.382134206312301"
## [1] "test statistic: 1067.76713443354"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   568.1  1067.8  1277.4  1336.0  1551.1  2908.2 
## [1] "p-value 0.75"
## [1] "VS"
## [1] "gcv 0.193398841583897"
## [1] "m 15 tau_n 0.382134206312301"
## [1] "test statistic: 103.342038019402"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   14.22   65.84  109.81  153.69  193.18 1254.33 
## [1] "p-value 0.5286"
## [1] "KS"
## [1] "gcv 0.193398841583897"
## [1] "m 8 tau_n 0.382134206312301"
## [1] "test statistic: 671.676091515897"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   223.4   558.0   723.3   787.2   950.5  2761.9 
## [1] "p-value 0.5712"

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.4406	0.3938	0.1212	0.5818
diff	0.4830	0.7500	0.5286	0.5712

xtable::xtable(pvmatrix, digits = 3)

## % latex table generated in R 4.6.0 by xtable 1.8-8 package
## % Thu Apr  2 16:14:16 2026
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
##   \hline
##  & KPSS & RS & VS & KS \\ 
##   \hline
## plug & 0.441 & 0.394 & 0.121 & 0.582 \\ 
##   diff & 0.483 & 0.750 & 0.529 & 0.571 \\ 
##    \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.1644"
## [1] "RS"
## [1] "p-value 0.1564"
## [1] "VS"
## [1] "p-value 0.1166"
## [1] "KS"
## [1] "p-value 0.2656"

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.543448031107"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15.19  103.87  201.95  323.50  409.03 3972.71 
## [1] "p-value 0.5742"
## [1] "RS"
## [1] "gcv 0.128932561055931"
## [1] "m 9 tau_n 0.332134206312301"
## [1] "test statistic: 998.08124125936"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     454    1001    1221    1284    1507    3716 
## [1] "p-value 0.7526"
## [1] "VS"
## [1] "gcv 0.128932561055931"
## [1] "m 15 tau_n 0.382134206312301"
## [1] "test statistic: 78.0587445148255"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   18.49  100.50  172.31  242.58  301.05 2240.43 
## [1] "p-value 0.843"
## [1] "KS"
## [1] "gcv 0.128932561055931"
## [1] "m 17 tau_n 0.332134206312301"
## [1] "test statistic: 709.345279801765"
## [1] "Bootstrap distribution"
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   378.5   820.3  1064.2  1151.7  1398.4  3816.3 
## [1] "p-value 0.8596"

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.1644	0.1564	0.1166	0.2656
diff	0.5742	0.7526	0.8430	0.8596

xtable::xtable(pvmatrix, digits = 3)

## % latex table generated in R 4.6.0 by xtable 1.8-8 package
## % Thu Apr  2 16:14:51 2026
## \begin{table}[ht]
## \centering
## \begin{tabular}{rrrrr}
##   \hline
##  & KPSS & RS & VS & KS \\ 
##   \hline
## plug & 0.164 & 0.156 & 0.117 & 0.266 \\ 
##   diff & 0.574 & 0.753 & 0.843 & 0.860 \\ 
##    \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 13 tau_n 0.414293094094381"
## [1] 10464.35
##        V1       
##  Min.   : 1843  
##  1st Qu.: 3990  
##  Median : 5110  
##  Mean   : 5467  
##  3rd Qu.: 6622  
##  Max.   :14819

p1

## [1] 0.0178

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.