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) > suppressMessages(library(dplyr)) > library(purrr) > library(tidyr) > suppressWarnings(RNGversion("3.5.3")) > > > plot_GIRF_line_low <- tsDyn:::plot_GIRF_line_low > irf_1_shock <- tsDyn:::irf_1_shock > irf_1_shock_ave <- tsDyn:::irf_1_shock_ave > > > > add.regime = FALSE > shock_both = TRUE > n.hist = 5 > n.ahead = 10 > n.shock = 20 > > ############################ > ### Load data > ############################ > path_mod_uni <- system.file("inst/testdata/models_univariate.rds", package = "tsDyn") > if(path_mod_uni=="") path_mod_uni <- system.file("testdata/models_univariate.rds", package = "tsDyn") > > path_mod_multi <- system.file("inst/testdata/models_multivariate.rds", package = "tsDyn") > if(path_mod_multi=="") path_mod_multi <- system.file("testdata/models_multivariate.rds", package = "tsDyn") > > models_ar_setar <- readRDS(path_mod_uni) %>% + filter(model %in% c("linear", "setar")) > > models_multivariate <- readRDS(path_mod_multi) > > ############################ > ### Univariate > ############################ > > > mod_1_uni <- models_ar_setar %>% + filter(model == "setar") %>% {.$object[[1]]} > > mod_1_ar <- models_ar_setar %>% + filter(model == "linear") %>% {.$object[[1]]} > > > mod_1_uni_1_shock <- irf_1_shock(object=mod_1_uni, + shock = 1, + hist = 0, + seed = 123) > > mod_ar_1_shock <- irf_1_shock(mod_1_ar, shock = 1, hist = 0,seed = 123) > > mod_1_uni_1_shock var n.ahead sim_1 sim_2 2 lh 0 2.996618 1.996618 3 lh 1 2.155985 1.628571 4 lh 2 1.664436 2.461223 5 lh 3 3.650510 3.025997 6 lh 4 3.101592 2.772216 7 lh 5 2.823117 2.649399 8 lh 6 2.720745 2.629123 9 lh 7 1.983854 1.935531 10 lh 8 1.560084 1.534598 11 lh 9 2.506257 2.483755 12 lh 10 2.696483 2.684616 > plot_GIRF_line_low(x=mod_1_uni_1_shock) > > plot_GIRF_line_low(x=mod_ar_1_shock) > plot(irf(mod_1_ar, boot = FALSE)) > > mod_1_uni_1_shock_ave <- irf_1_shock_ave(object = mod_1_uni, + shock = 1, + hist = 0, + seed = 123) > > plot_GIRF_line_low(x=mod_1_uni_1_shock_ave) > > mod_1_uni_1_girf <- GIRF(object = mod_1_uni, + hist_li = list(1), + shock_li = list(0.1, 0.11, -0.1, -0.11), + R = 2, + seed = 123) > > head(mod_1_uni_1_girf) n_simu hist_x1_l1 shock_var1 n.ahead var sim_1 sim_2 girf 1 1 1 0.1 0 lh 2.848254 2.748254 0.100000000 2 1 1 0.1 1 lh 2.284368 2.231627 0.052741426 3 1 1 0.1 2 lh 1.788436 1.760620 0.027816580 4 1 1 0.1 3 lh 2.518312 2.498697 0.019615216 5 1 1 0.1 4 lh 3.239090 3.226137 0.012953068 6 1 1 0.1 5 lh 2.617766 2.610934 0.006831633 > plot_GIRF_line_low(mod_1_uni_1_girf, n_simu = 1:4) > > mod_1_uni_1_girf_big <- GIRF(object = mod_1_uni, + n.hist = 5, + R = 2, + seed = 123) > > plot(x=mod_1_uni_1_girf_big, plot_type = "density") > plot(x=mod_1_uni_1_girf_big, plot_type = "line", n_simu = 1:50, add_legend =FALSE) > > ## Simple, given shocks > models_ar_setar %>% + mutate(girf = map2(object, lag, ~GIRF(object=.x, n.ahead = 3, + hist_li = list(rep(1.6, .y)), + shock_li = list(0.01), R = 2, seed = 123) %>% as_tibble)) %>% + select(-object) %>% + unnest(girf) %>% + as.data.frame() %>% + head(10) %>% + print(digits=3) lag include model nthresh thDelay n_simu hist_x1_l1 shock_var1 n.ahead var 1 1 both linear NA NA 1 1.6 0.01 0 x 2 1 both linear NA NA 1 1.6 0.01 1 x 3 1 both linear NA NA 1 1.6 0.01 2 x 4 1 both linear NA NA 1 1.6 0.01 3 x 5 1 const linear NA NA 1 1.6 0.01 0 x 6 1 const linear NA NA 1 1.6 0.01 1 x 7 1 const linear NA NA 1 1.6 0.01 2 x 8 1 const linear NA NA 1 1.6 0.01 3 x 9 1 none linear NA NA 1 1.6 0.01 0 x 10 1 none linear NA NA 1 1.6 0.01 1 x sim_1 sim_2 girf hist_x1_l2 1 2.92 2.91 0.01000 NA 2 2.18 2.17 0.00529 NA 3 1.74 1.74 0.00280 NA 4 2.63 2.63 0.00148 NA 5 3.10 3.09 0.01000 NA 6 2.48 2.47 0.00586 NA 7 2.03 2.02 0.00343 NA 8 3.01 3.00 0.00201 NA 9 2.97 2.96 0.01000 NA 10 2.80 2.79 0.00984 NA > > ## Simple, random > models_ar_setar %>% + mutate(girf = map(object, ~GIRF(object=., n.ahead = 3, n.hist = 3, n.shock = 3, + R = 2, seed = 123) %>% as_tibble)) %>% + select(-object) %>% + unnest(girf) %>% + as.data.frame() %>% + head(10) %>% + print(digits=3) lag include model nthresh thDelay n_simu hist_x1_l1 shock_var1 n.ahead var 1 1 both linear NA NA 1 1.8 1.19 0 x 2 1 both linear NA NA 1 1.8 1.19 1 x 3 1 both linear NA NA 1 1.8 1.19 2 x 4 1 both linear NA NA 1 1.8 1.19 3 x 5 1 both linear NA NA 2 1.4 1.19 0 x 6 1 both linear NA NA 2 1.4 1.19 1 x 7 1 both linear NA NA 2 1.4 1.19 2 x 8 1 both linear NA NA 2 1.4 1.19 3 x 9 1 both linear NA NA 3 2.2 1.19 0 x 10 1 both linear NA NA 3 2.2 1.19 1 x sim_1 sim_2 girf hist_x1_l2 1 4.20 3.01 1.185 NA 2 2.86 2.23 0.627 NA 3 2.10 1.77 0.332 NA 4 2.82 2.65 0.175 NA 5 3.99 2.80 1.185 NA 6 2.74 2.12 0.627 NA 7 2.04 1.71 0.332 NA 8 2.79 2.62 0.175 NA 9 4.41 3.22 1.185 NA 10 2.97 2.34 0.627 NA > > > > ############################ > ### Multivariate > ############################ > > mod_TVAR <- models_multivariate %>% + filter(model == "TVAR" & lag ==2) %>% {.$object[[1]]} > > mod_VAR <- models_multivariate %>% + filter(model == "VAR"& lag ==2) %>% {.$object[[1]]} > > mod_VECM <- models_multivariate %>% + filter(model == "VECM" & lag ==2) %>% {.$object[[1]]} > > mod_TVAR_1_shock <- irf_1_shock(object = mod_TVAR, + shock = matrix(c(1, 0), nrow = 1), + hist = matrix(c(0, 0, 0, 0), nrow = 2), + seed = 123) > > mod_VAR_1_shock <- irf_1_shock(object = mod_VAR, shock = matrix(c(1, 0), nrow = 1), + hist = matrix(c(0, 0, 0, 0), nrow = 2), seed = 123) > > mod_VECM_1_shock <- irf_1_shock(object = mod_VECM, shock = matrix(c(1, 0), nrow = 1), + hist = matrix(rep(0, 6), nrow = 3), seed = 123) > > > plot_GIRF_line_low(x=mod_VAR_1_shock) > plot_GIRF_line_low(x=mod_TVAR_1_shock) > > plot_GIRF_line_low(x=mod_TVAR_1_shock, var = "cpiUSA") > plot_GIRF_line_low(x=mod_VAR_1_shock, var = "cpiUSA") > > mod_TVAR_1_shock_ave <- irf_1_shock_ave(mod_TVAR, + shock = matrix(c(1, 0), nrow = 1), + hist = matrix(c(0, 0, 0, 0), nrow = 2), + seed = 123) > > > plot_GIRF_line_low(x=mod_TVAR_1_shock_ave) > > TVAR_GIRF <- GIRF(object=mod_TVAR, + shock_li = list(matrix(c(1, 0), nrow = 1), + matrix(c(0, 1), nrow = 1)), + hist_li = list(matrix(c(0, 0, 0, 0), nrow = 2), + matrix(c(0, 1, 0, 0), nrow = 2)), + R = 2, + seed = 123) > > gi_out <- GIRF(object=mod_TVAR, seed = 123, n.hist = 40, R = 2) > plot(density(residuals(mod_TVAR)[, 1])) > head(gi_out) n_simu hist_x1_l1 hist_x1_l2 hist_x2_l1 hist_x2_l2 shock_var1 shock_var2 1 1 1.1647 1.1691 55.14 55.63 0.01642286 0.08589973 2 1 1.1647 1.1691 55.14 55.63 0.01642286 0.08589973 3 1 1.1647 1.1691 55.14 55.63 0.01642286 0.08589973 4 1 1.1647 1.1691 55.14 55.63 0.01642286 0.08589973 5 1 1.1647 1.1691 55.14 55.63 0.01642286 0.08589973 6 1 1.1647 1.1691 55.14 55.63 0.01642286 0.08589973 n.ahead var sim_1 sim_2 girf 1 0 cpiUSA 56.13770 56.05180 0.08589973 2 1 cpiUSA 56.47797 56.37500 0.10297148 3 2 cpiUSA 56.73699 56.62066 0.11633341 4 3 cpiUSA 56.93175 56.80123 0.13052292 5 4 cpiUSA 57.21324 57.06854 0.14469542 6 5 cpiUSA 57.35781 57.19943 0.15837356 > plot(x=gi_out, var = "dolcan") > plot(x=gi_out, var = "cpiUSA") > > ## Simple, random > models_multivariate %>% + # head(2) %>% + mutate(girf = map(object, ~GIRF(object=., n.ahead = 3, n.hist = 3, n.shock = 3, + R = 2, seed = 123) %>% head(2))) %>% + select(-object) %>% + unnest(girf) %>% + select(-object_vars) %>% + as.data.frame() %>% + slice(10:20) %>% + print(digits=3) lag include model nthresh n_simu hist_x1_l1 hist_x2_l1 shock_var1 shock_var2 1 2 both VAR NA 1 1.16 55.1 0.0162 0.129 2 2 const VAR NA 1 1.16 55.1 0.0164 0.127 3 2 const VAR NA 1 1.16 55.1 0.0164 0.127 4 2 none VAR NA 1 1.16 55.1 0.0161 0.122 5 2 none VAR NA 1 1.16 55.1 0.0161 0.122 6 2 trend VAR NA 1 1.16 55.1 0.0156 0.116 7 2 trend VAR NA 1 1.16 55.1 0.0156 0.116 8 1 both VECM NA 1 1.16 55.1 0.0161 0.124 9 1 both VECM NA 1 1.16 55.1 0.0161 0.124 10 1 const VECM NA 1 1.16 55.1 0.0163 0.126 11 1 const VECM NA 1 1.16 55.1 0.0163 0.126 n.ahead var sim_1 sim_2 girf hist_x1_l2 hist_x2_l2 hist_x1_l3 hist_x2_l3 1 1 cpiUSA 56.4 56.2 0.197 1.17 55.6 NA NA 2 0 cpiUSA 56.2 56.0 0.127 1.17 55.6 NA NA 3 1 cpiUSA 56.5 56.3 0.193 1.17 55.6 NA NA 4 0 cpiUSA 56.2 56.0 0.122 1.17 55.6 NA NA 5 1 cpiUSA 56.5 56.3 0.193 1.17 55.6 NA NA 6 0 cpiUSA 56.3 56.2 0.116 1.17 55.6 NA NA 7 1 cpiUSA 56.7 56.6 0.180 1.17 55.6 NA NA 8 0 cpiUSA 56.2 56.1 0.124 1.17 55.6 NA NA 9 1 cpiUSA 56.6 56.4 0.189 1.17 55.6 NA NA 10 0 cpiUSA 56.2 56.0 0.126 1.17 55.6 NA NA 11 1 cpiUSA 56.5 56.3 0.193 1.17 55.6 NA NA > > > proc.time() user system elapsed 10.26 0.45 10.71