Movatterモバイル変換


[0]ホーム

URL:


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)/nhk=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)###inistializationsetting=list(B =5000,gcv =1,neighbour =1)setting$lb=floor(10/7*n^(4/15))- setting$neighboursetting$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 matrixfunction.

setting$lrvmethod=0.i=1# print(rule_of_thumb(y= hk$y, x = hk$x))for(typeinc("KPSS","RS","VS","KS")){  setting$type= typeprint(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 matrixfunction.

setting$lrvmethod=1i=1for(typeinc("KPSS","RS","VS","KS")){  setting$type= typeprint(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")
KPSSRSVSKS
plug0.28380.28960.11480.4124
diff0.68660.80660.51260.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 afactor. - Shift = 1.2 with plug-in estimator.

pvmatrix=matrix(nrow=2,ncol=4)setting$lrvmethod=0i=1for(typeinc("KPSS","RS","VS","KS")){  setting$type= typeprint(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"
setting$lrvmethod=1i=1for(typeinc("KPSS","RS","VS","KS")){  setting$type= typeprint(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"
rownames(pvmatrix)=c("plug","diff")colnames(pvmatrix)=c("KPSS","RS","VS","KS")knitr::kable(pvmatrix,type="latex")
KPSSRSVSKS
plug0.42560.36380.1180.5610
diff0.64620.69940.5380.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}
pvmatrix=matrix(nrow=2,ncol=4)setting$lrvmethod=0i=1for(typeinc("KPSS","RS","VS","KS")){  setting$type= typeprint(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"
setting$lrvmethod=1i=1for(typeinc("KPSS","RS","VS","KS")){  setting$type= typeprint(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"
rownames(pvmatrix)=c("plug","diff")colnames(pvmatrix)=c("KPSS","RS","VS","KS")knitr::kable(pvmatrix,type="latex")
KPSSRSVSKS
plug0.24580.16620.12340.2726
diff0.57100.94640.65860.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 secondyear is constant.

hk$x=as.matrix(cbind(rep(1,n), (hk_data[,1:3])))hk$y= hk_data$`Hospital Admission`setting$type=0setting$bw_set=c(0.1,0.35)setting$eta=0.2setting$lrvmethod=1setting$lb=10setting$ub=15hk1=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

[8]ページ先頭

©2009-2025 Movatter.jp