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("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.01 0.64 10.61