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.

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)

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`

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.2862"
## [1] "RS"
## [1] "p-value 0.3022"
## [1] "VS"
## [1] "p-value 0.126"
## [1] "KS"
## [1] "p-value 0.4046"

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.7034"
## [1] "RS"
## [1] "p-value 0.8002"
## [1] "VS"
## [1] "p-value 0.7702"
## [1] "KS"
## [1] "p-value 0.8612"

Output

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 
xtable::xtable(pvmatrix, digits = 3)
## % 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}

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.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 
xtable::xtable(pvmatrix, digits = 3)
## % 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 
xtable::xtable(pvmatrix, digits = 3)
## % 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 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 11 tau_n 0.414293094094381"
## [1] 10464.35
##        V1       
##  Min.   : 1343  
##  1st Qu.: 3343  
##  Median : 4328  
##  Mean   : 4674  
##  3rd Qu.: 5633  
##  Max.   :13336
p1
## [1] 0.0058

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.