R Under development (unstable) (2024-10-26 r87273 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(tsDyn) > library(mnormt) > suppressWarnings(RNGversion("3.5.3")) > > TVECM.boot.check <- tsDyn:::TVECM.boot.check > options(useFancyQuotes=FALSE) # useful for all.equal comparison > > ################################################################ > ######### From man file: > ################################################################ > #see that: > a<-matrix(c(-0.2, 0.2), ncol=1) > b<-matrix(c(1,-1), nrow=1) > a%*%b [,1] [,2] [1,] -0.2 0.2 [2,] 0.2 -0.2 > > set.seed(123) > innov<-rmnorm(100, varcov=diag(2)) > vecm1<-TVECM.sim(B=rbind(c(-0.2, 0,0), c(0.2, 0,0)), nthresh=0, beta=1,n=100, lag=1,include="none", innov=innov) > ECT<-vecm1[,1]-vecm1[,2] > ECT[1:5] [1] -0.3302982 1.2900210 -0.8117646 1.2389187 0.5021603 > > #add an intercept as in panel B > B2 <- rbind(c(-0.2, 0.1,0,0), c(0.2,0.4, 0,0)) > b<- TVECM.sim(B=B2, nthresh=0, n=100,beta=1, lag=1,include="const", innov=innov, show.parMat=TRUE) ECT Const Trend L{x1}{1} L{x2}{1} Equ x1: -0.2 0.1 0 0 0 Equ x2: 0.2 0.4 0 0 0 > b_2 <- VECM.sim(B=B2, n=100,beta=1, lag=1,include="const", innov=innov, show.parMat=TRUE) ECT Const Trend L{x1}{1} L{x2}{1} Equ x1: -0.2 0.1 0 0 0 Equ x2: 0.2 0.4 0 0 0 > b[1:5,] [,1] [,2] [1,] -0.4604756 0.1698225 [2,] 1.3242923 0.5142713 [3,] 1.3915758 2.7913405 [4,] 2.2324450 1.6463263 [5,] 1.5283684 1.7178881 > all.equal(b, b_2) [1] TRUE > > ## other ways to input beta: > b_beta_vec <- TVECM.sim(B=B2, nthresh=0, n=100,beta=c(1,-1), lag=1,include="const", innov=innov, show.parMat=TRUE) ECT Const Trend L{x1}{1} L{x2}{1} Equ x1: -0.2 0.1 0 0 0 Equ x2: 0.2 0.4 0 0 0 > b_beta_vec_2 <- VECM.sim(B=B2, n=100,beta=c(1,-1), lag=1,include="const", innov=innov, show.parMat=TRUE) ECT Const Trend L{x1}{1} L{x2}{1} Equ x1: -0.2 0.1 0 0 0 Equ x2: 0.2 0.4 0 0 0 > b_beta_vec[1:5,] [,1] [,2] [1,] -0.4604756 0.1698225 [2,] 1.3242923 0.5142713 [3,] 1.3915758 2.7913405 [4,] 2.2324450 1.6463263 [5,] 1.5283684 1.7178881 > all.equal(b,b_beta_vec) [1] TRUE > all.equal(b_beta_vec, b_beta_vec_2) [1] TRUE > > beta_mat <- matrix(c(1,-1), ncol=1) > b_beta_mat <- TVECM.sim(B=B2, nthresh=0, n=100,beta=beta_mat, lag=1,include="const", innov=innov, show.parMat=TRUE) ECT Const Trend L{x1}{1} L{x2}{1} Equ x1: -0.2 0.1 0 0 0 Equ x2: 0.2 0.4 0 0 0 > b_beta_mat[1:5,] [,1] [,2] [1,] -0.4604756 0.1698225 [2,] 1.3242923 0.5142713 [3,] 1.3915758 2.7913405 [4,] 2.2324450 1.6463263 [5,] 1.5283684 1.7178881 > all.equal(b,b_beta_mat) [1] TRUE > > > ######################## > ######## TVECM > ######################## > > ##Bootstrap a TVECM with 1 threshold (two regimes) > data(zeroyld) > dat<-zeroyld > TVECMobject<-TVECM(dat, nthresh=1, lag=1, ngridBeta=20, th1=list(exact=-1.195), plot=FALSE) > tv_1_boot <-TVECM.boot(TVECMobject, seed=123, show.parMat=TRUE) ECT_low Const_low Trend_low L{x1}{1}_low L{x2}{1}_low ECT_upr Equ x1: 0.2025264 0.4857667 0 0.2252923 -0.09926611 0.001291099 Equ x2: 0.9575983 1.5929831 0 0.8703592 0.02252604 0.058816293 Const_upr Trend_upr L{x1}{1}_upr L{x2}{1}_upr Equ x1: 0.004290751 0 -0.02328976 0.06741675 Equ x2: -0.007954949 0 0.16032265 0.14668443 > head(tv_1_boot) short.run long.run [1,] 2.1830000 1.57500000 [2,] 2.2460000 1.54500000 [3,] 2.2373746 1.40316599 [4,] 1.6102958 0.94029117 [5,] 1.4410637 0.65045167 [6,] 0.6969609 -0.06715804 > > ##Check the bootstrap > TVECM.boot.check(TVECMobject) [1] TRUE > > ## check correspondence bootstrap/simul: > tv_1_sim <-TVECM.sim(B=tsDyn:::coefMat.nlVar(TVECMobject),beta=TVECMobject$model.specific$beta, + Thresh=getTh(TVECMobject), show.parMat=TRUE, innov=matrix(0,200,2)) ECT_low Const_low Trend_low L{x1}{1}_low L{x2}{1}_low ECT_upr Equ x1: 0.2025264 0.4857667 0 0.2252923 -0.09926611 0.001291099 Equ x2: 0.9575983 1.5929831 0 0.8703592 0.02252604 0.058816293 Const_upr Trend_upr L{x1}{1}_upr L{x2}{1}_upr Equ x1: 0.004290751 0 -0.02328976 0.06741675 Equ x2: -0.007954949 0 0.16032265 0.14668443 > head(tv_1_boot) short.run long.run [1,] 2.1830000 1.57500000 [2,] 2.2460000 1.54500000 [3,] 2.2373746 1.40316599 [4,] 1.6102958 0.94029117 [5,] 1.4410637 0.65045167 [6,] 0.6969609 -0.06715804 > > tv_1_sim <-TVECM.sim(B=tsDyn:::coefMat.nlVar(TVECMobject), + beta=TVECMobject$model.specific$beta, + Thresh=getTh(TVECMobject), show.parMat=TRUE) ECT_low Const_low Trend_low L{x1}{1}_low L{x2}{1}_low ECT_upr Equ x1: 0.2025264 0.4857667 0 0.2252923 -0.09926611 0.001291099 Equ x2: 0.9575983 1.5929831 0 0.8703592 0.02252604 0.058816293 Const_upr Trend_upr L{x1}{1}_upr L{x2}{1}_upr Equ x1: 0.004290751 0 -0.02328976 0.06741675 Equ x2: -0.007954949 0 0.16032265 0.14668443 > head(tv_1_boot) short.run long.run [1,] 2.1830000 1.57500000 [2,] 2.2460000 1.54500000 [3,] 2.2373746 1.40316599 [4,] 1.6102958 0.94029117 [5,] 1.4410637 0.65045167 [6,] 0.6969609 -0.06715804 > > ##Bootstrap a TVAR with two threshold (three regimes) > tv_2_const <- TVECM(dat, lag=1, nthresh=2, plot=FALSE, trace=FALSE, th1=list(exact=-1.312), + th2=list(exact=0.774)) > tv_2_none <- TVECM(dat, lag=1, nthresh=2, plot=FALSE, trace=FALSE, th1=list(exact=-1.312), + th2=list(exact=0.774), include="none") > tv_2_trend <- TVECM(dat, lag=1, nthresh=2, plot=FALSE, trace=FALSE, th1=list(exact=-1.312), + th2=list(exact=0.774), include="trend") > tv_2_both <- TVECM(dat, lag=1, nthresh=2, plot=FALSE, trace=FALSE, th1=list(exact=-1.312), + th2=list(exact=0.774), ngridBeta=5, include="both") > > tv_2_const_common <- TVECM(dat, lag=1, nthresh=2, plot=FALSE, trace=FALSE, th1=list(exact=-1.451), + th2=list(exact=0.754), include="const", common="only_ECT") > tv_2_const_l2 <- TVECM(dat, nthresh=2, lag=2, ngridBeta=5, plot=FALSE, include="none", + th1=list(exact=-1.312),th2=list(exact=0.774), trace=FALSE) > > TVECM.boot(tv_2_const, show.parMat=TRUE, seed=456)[1:5,] ECT_low Const_low Trend_low L{x1}{1}_low L{x2}{1}_low ECT_mid Equ x1: 0.2343191 0.5627267 0 0.2133872 -0.13045497 -0.01031875 Equ x2: 1.1492230 1.9445867 0 0.7652367 0.02490849 0.01991474 Const_mid Trend_mid L{x1}{1}_mid L{x2}{1}_mid ECT_upr Equ x1: 0.007322243 0 0.09944158 0.06907443 -0.04920032 Equ x2: -0.005735839 0 0.32699066 0.23687489 0.10732160 Const_upr Trend_upr L{x1}{1}_upr L{x2}{1}_upr Equ x1: 0.06691590 0 -0.11216037 0.04637838 Equ x2: -0.09079342 0 0.08524092 -0.02873180 short.run long.run [1,] 2.183000 1.575000 [2,] 2.246000 1.545000 [3,] 2.255825 1.664895 [4,] 2.289260 1.701369 [5,] 2.151199 1.412521 > TVECM.boot(tv_2_none, show.parMat=TRUE, seed=456)[1:5,] ECT_low Const_low Trend_low L{x1}{1}_low L{x2}{1}_low ECT_mid Equ x1: -0.03252245 0 0 0.3508268 -0.2229371 -0.01855284 Equ x2: 0.17029747 0 0 1.2574671 -0.3689984 -0.02128513 Const_mid Trend_mid L{x1}{1}_mid L{x2}{1}_mid ECT_upr Const_upr Equ x1: 0 0 0.09248489 0.05361497 0.008401237 0 Equ x2: 0 0 0.25119037 0.19618437 0.061539636 0 Trend_upr L{x1}{1}_upr L{x2}{1}_upr Equ x1: 0 -0.18027015 0.1024449 Equ x2: 0 0.01649033 0.1198955 short.run long.run [1,] 2.183000 1.575000 [2,] 2.246000 1.545000 [3,] 2.228665 1.598876 [4,] 2.254419 1.609751 [5,] 2.198737 1.570296 > TVECM.boot(tv_2_trend,show.parMat=TRUE, seed=456)[1:5,] ECT_low Const_low Trend_low L{x1}{1}_low L{x2}{1}_low ECT_mid Equ x1: 0.3391893 0 0.002350986 0.07245793 -0.05503711 0 Equ x2: 1.4764723 0 0.007927246 0.30567992 0.26850260 0 Const_mid Trend_mid L{x1}{1}_mid L{x2}{1}_mid ECT_upr Equ x1: -0.009627084 -1.272937e-05 0.0967326 0 0.07017931 Equ x2: 0.019137028 -2.862608e-05 0.3246202 0 0.23703865 Const_upr Trend_upr L{x1}{1}_upr L{x2}{1}_upr Equ x1: -0.004106742 1.678226e-05 -0.11558062 0.04875773 Equ x2: 0.080441375 -1.830732e-04 0.08127482 -0.02692949 short.run long.run [1,] 2.183000 1.575000 [2,] 2.246000 1.545000 [3,] 2.251418 1.638299 [4,] 2.273142 1.657000 [5,] 2.133514 1.377453 > TVECM.boot(tv_2_both, show.parMat=TRUE, seed=456)[1:5,] ECT_low Const_low Trend_low L{x1}{1}_low L{x2}{1}_low ECT_mid Equ x1: 0.3344445 -0.1764324 0.002892201 0.06269098 -0.04669461 -0.01913800 Equ x2: 1.4605964 -0.2530547 0.008833043 0.30094440 0.27609545 0.02432369 Const_mid Trend_mid L{x1}{1}_mid L{x2}{1}_mid ECT_upr Equ x1: 0.0329355255 -1.274667e-04 0.06886216 0.07908553 -0.05978909 Equ x2: 0.0003496974 -2.007436e-05 0.29465882 0.24365050 0.12349362 Const_upr Trend_upr L{x1}{1}_upr L{x2}{1}_upr Equ x1: 0.09243691 -3.868788e-05 -0.11591998 0.04050465 Equ x2: -0.06383229 -1.475473e-04 0.07572181 -0.02512837 short.run long.run [1,] 2.183000 1.575000 [2,] 2.246000 1.545000 [3,] 2.265084 1.629432 [4,] 2.303443 1.670040 [5,] 2.166667 1.409473 > try(TVECM.boot(tv_2_const_common, show.parMat=TRUE, seed=456)[1:5,], silent=TRUE) > > TVECM.boot.check(tv_2_const) [1] TRUE > TVECM.boot.check(tv_2_none) [1] TRUE > TVECM.boot.check(tv_2_const_l2) [1] TRUE > TVECM.boot.check(tv_2_both) [1] TRUE > > ## does not work: > TVECM.boot.check(tv_2_trend) [1] "Component \"short.run\": Mean relative difference: 0.2180451" [2] "Component \"long.run\": Mean relative difference: 0.1916152" > > ############### > #### p>2 > ############### > > data(barry) > > ve_r1_l1 <- VECM(barry, lag=1) > ve_r1_l3 <- VECM(barry, lag=3) > ve_r2_l3 <- VECM(barry, lag=3, estim="ML", r=2) > > TVECM.boot.check(ve_r1_l1) [1] TRUE > TVECM.boot.check(ve_r1_l3) [1] TRUE > TVECM.boot.check(ve_r2_l3) [1] TRUE > > VECM.boot(ve_r1_l1, show.parMat=TRUE, seed=234)[1:5,] ECT Const Trend L{x1}{1} L{x2}{1} L{x3}{1} Equ x1: -0.002537673 0.002202355 0 0.1587853 -0.003002494 0.0001286435 Equ x2: 0.055320874 0.135043876 0 0.5426890 0.483673616 -0.0353661496 Equ x3: 0.179057392 0.157670998 0 0.2897487 0.327296468 -0.0224027426 dolcan cpiUSA cpiCAN [1,] 0.9993000 27.98000 25.86000 [2,] 0.9956000 28.17000 26.00000 [3,] 0.9836973 28.47662 25.99240 [4,] 1.0075128 28.86076 25.35688 [5,] 1.0154009 28.92388 25.50034 > VECM.boot(ve_r1_l3, show.parMat=TRUE, seed=234)[1:5,] ECT Const Trend L{x1}{1} L{x2}{1} L{x3}{1} Equ x1: -0.001338343 0.003662911 0 0.1611940 -0.002242126 0.0003953115 Equ x2: 0.029339845 0.108353256 0 0.8711430 0.468684873 -0.0510243367 Equ x3: 0.133099298 0.100383436 0 0.6925224 0.272239507 -0.0494533196 L{x1}{2} L{x2}{2} L{x3}{2} L{x1}{3} L{x2}{3} Equ x1: -0.04139589 0.002012077 -0.006607018 0.0202477 -0.002153848 Equ x2: -0.49404335 -0.023189036 0.071589503 -0.3943149 0.007943649 Equ x3: -0.13395315 0.009117886 0.164960151 -0.6081994 0.022337674 L{x3}{3} Equ x1: -0.0003941692 Equ x2: 0.0928680566 Equ x3: 0.1324489696 dolcan cpiUSA cpiCAN [1,] 0.9993000 27.98000 25.86000 [2,] 0.9956000 28.17000 26.00000 [3,] 0.9967000 28.44000 26.07000 [4,] 1.0007000 28.63000 26.36000 [5,] 0.9916491 28.93368 26.32038 > VECM.boot(ve_r2_l3, show.parMat=TRUE, seed=234)[1:5,] ECT ECT Const Trend L{x1}{1} L{x2}{1} Equ x1: -0.02128737 0.0005968341 0.02324814 0 0.1622498 -0.003419996 Equ x2: -0.05886126 0.0019115449 0.19169805 0 0.8619558 0.462197612 Equ x3: 0.13234280 -0.0025774055 0.10561308 0 0.6225139 0.266045773 L{x3}{1} L{x1}{2} L{x2}{2} L{x3}{2} L{x1}{3} Equ x1: 0.0006640453 -0.03681695 0.001110734 -0.006387489 0.02526483 Equ x2: -0.0485913611 -0.48826864 -0.028412522 0.073825879 -0.38942643 Equ x3: -0.0476408945 -0.19364163 0.004214840 0.166485918 -0.68004445 L{x2}{3} L{x3}{3} Equ x1: -0.0036154386 -0.0002821687 Equ x2: -0.0005587153 0.0944265013 Equ x3: 0.0141090720 0.1337360507 dolcan cpiUSA cpiCAN [1,] 0.9993000 27.98000 25.86000 [2,] 0.9956000 28.17000 26.00000 [3,] 0.9967000 28.44000 26.07000 [4,] 1.0007000 28.63000 26.36000 [5,] 0.9939164 28.94688 26.34743 > > ################################################################ > ######### Check error message when matrix badly specified: > ################################################################ > > B <- matrix(rnorm(14), byrow=TRUE,ncol=7) > > ## 0 thresh > try(a<-TVECM.sim(B=B, beta=1, nthresh=0, n=100, lag=1,show.parMat=TRUE, include="none")) Matrix B badly specified: should have 3 columns, but has 7 ECT L{x1}{1} L{x2}{1} Equ x1 NA NA NA Equ x2 NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > try(a<-TVECM.sim(B=B, beta=1, nthresh=0, n=100, lag=1,show.parMat=TRUE, include="const")) Matrix B badly specified: should have 4 columns, but has 7 ECT const L{x1}{1} L{x2}{1} Equ x1 NA NA NA NA Equ x2 NA NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > try(a<-TVECM.sim(B=B, beta=1, nthresh=0, n=100, lag=1,show.parMat=TRUE, include="both")) Matrix B badly specified: should have 5 columns, but has 7 ECT const trend L{x1}{1} L{x2}{1} Equ x1 NA NA NA NA NA Equ x2 NA NA NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > try(a<-TVECM.sim(B=B, beta=1, nthresh=0, n=100, lag=1,show.parMat=TRUE, include="trend")) Matrix B badly specified: should have 4 columns, but has 7 ECT trend L{x1}{1} L{x2}{1} Equ x1 NA NA NA NA Equ x2 NA NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > > try(a<-TVECM.sim(B=B, beta=1, nthresh=0, n=100, lag=2,show.parMat=TRUE, include="none")) Matrix B badly specified: should have 5 columns, but has 7 ECT L{x1}{1} L{x2}{1} L{x1}{2} L{x2}{2} Equ x1 NA NA NA NA NA Equ x2 NA NA NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > try(a<-TVECM.sim(B=B, beta=1, nthresh=0, n=100, lag=2,show.parMat=TRUE, include="const")) Matrix B badly specified: should have 6 columns, but has 7 ECT const L{x1}{1} L{x2}{1} L{x1}{2} L{x2}{2} Equ x1 NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > > ## 1 thresh > try(a<-TVECM.sim(B=B, beta=1, nthresh=1, Thresh=0, n=100, lag=1,show.parMat=TRUE, include="none")) Matrix B badly specified: should have 6 columns, but has 7 ECT_low L{x1}{1}_low L{x2}{1}_low ECT_upr L{x1}{1}_upr L{x2}{1}_upr Equ x1 NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > try(a<-TVECM.sim(B=B, beta=1, nthresh=1, Thresh=0, n=100, lag=1,show.parMat=TRUE, include="const")) Matrix B badly specified: should have 8 columns, but has 7 ECT_low const_low L{x1}{1}_low L{x2}{1}_low ECT_upr const_upr Equ x1 NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA L{x1}{1}_upr L{x2}{1}_upr Equ x1 NA NA Equ x2 NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > try(a<-TVECM.sim(B=B, beta=1, nthresh=1, Thresh=0, n=100, lag=1,show.parMat=TRUE, include="both")) Matrix B badly specified: should have 10 columns, but has 7 ECT_low const_low trend_low L{x1}{1}_low L{x2}{1}_low ECT_upr const_upr Equ x1 NA NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA NA trend_upr L{x1}{1}_upr L{x2}{1}_upr Equ x1 NA NA NA Equ x2 NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > try(a<-TVECM.sim(B=B, beta=1, nthresh=1, Thresh=0, n=100, lag=1,show.parMat=TRUE, include="trend")) Matrix B badly specified: should have 8 columns, but has 7 ECT_low trend_low L{x1}{1}_low L{x2}{1}_low ECT_upr trend_upr Equ x1 NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA L{x1}{1}_upr L{x2}{1}_upr Equ x1 NA NA Equ x2 NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > > try(a<-TVECM.sim(B=B, beta=1, nthresh=1, Thresh=0, n=100, lag=2,show.parMat=TRUE, include="const")) Matrix B badly specified: should have 12 columns, but has 7 ECT_low const_low L{x1}{1}_low L{x2}{1}_low L{x1}{2}_low L{x2}{2}_low Equ x1 NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA ECT_upr const_upr L{x1}{1}_upr L{x2}{1}_upr L{x1}{2}_upr L{x2}{2}_upr Equ x1 NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > > ## 2 thresh > try(a<-TVECM.sim(B=B, beta=1, nthresh=2, Thresh=0, n=100, lag=1,show.parMat=TRUE, include="none")) Matrix B badly specified: should have 9 columns, but has 7 ECT_low L{x1}{1}_low L{x2}{1}_low ECT_mid L{x1}{1}_mid L{x2}{1}_mid Equ x1 NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA ECT_upr L{x1}{1}_upr L{x2}{1}_upr Equ x1 NA NA NA Equ x2 NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > try(a<-TVECM.sim(B=B, beta=1, nthresh=2, Thresh=0, n=100, lag=2,show.parMat=TRUE, include="const")) Matrix B badly specified: should have 18 columns, but has 7 ECT_low const_low L{x1}{1}_low L{x2}{1}_low L{x1}{2}_low L{x2}{2}_low Equ x1 NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA ECT_mid const_mid L{x1}{1}_mid L{x2}{1}_mid L{x1}{2}_mid L{x2}{2}_mid Equ x1 NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA ECT_upr const_upr L{x1}{1}_upr L{x2}{1}_upr L{x1}{2}_upr L{x2}{2}_upr Equ x1 NA NA NA NA NA NA Equ x2 NA NA NA NA NA NA Error in TVECM.gen(B = B, n = n, lag = lag, include = include, beta = beta, : > > > proc.time() user system elapsed 2.31 0.43 2.73