R Under development (unstable) (2024-10-28 r87274 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")) > > roundAll.Equal <- function(x, round=8){ + isFALSE <- !isTRUE(x) + xFalse <- x[isFALSE] + # extract the number (i.e remove all the rest) + xf<- gsub("(Component ([0-9]+)?([[:punct:]][[:alnum:]]+[[:punct:]])?: )?Mean relative difference: ", + "", xFalse) + xf2<- round(as.numeric(xf),round) + x[isFALSE] <- paste("Mean relative difference at tol ", round, ": ", xf2, sep="") + x + } > > ############################ > ### 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") > > models_univariate <- readRDS(path_mod_uni) %>% + filter(model %in% c("lineaar", "setar")) > > ################ > ### run boot > ################ > > models_univariate %>% + mutate(boot = map(object, ~ tibble(n = 1:3, boot =head(setar.boot(., seed = 123), 3)))) %>% + select(-object) %>% + unnest(boot) %>% + spread(include, boot) %>% + as.data.frame() %>% + print(digits=3) lag model nthresh thDelay n both const none trend 1 1 setar 1 0 1 2.68 3.26 3.27 2.58 2 1 setar 1 0 2 2.81 3.32 3.46 2.95 3 1 setar 1 0 3 2.15 2.71 2.83 2.60 4 1 setar 2 0 1 3.26 3.31 3.30 3.26 5 1 setar 2 0 2 3.41 3.21 3.50 3.79 6 1 setar 2 0 3 3.09 2.80 2.90 3.22 7 2 setar 1 0 1 2.25 2.29 2.25 2.15 8 2 setar 1 0 2 2.44 2.35 2.45 2.47 9 2 setar 1 0 3 1.96 1.87 1.92 2.25 10 2 setar 1 1 1 2.04 2.46 2.19 2.28 11 2 setar 1 1 2 2.16 2.68 2.50 2.23 12 2 setar 1 1 3 1.86 2.34 2.14 1.59 13 2 setar 2 0 1 2.26 2.12 2.13 2.20 14 2 setar 2 0 2 2.44 2.53 2.52 2.70 15 2 setar 2 0 3 2.12 2.22 2.29 2.15 16 2 setar 2 1 1 1.99 2.20 2.21 2.06 17 2 setar 2 1 2 2.06 2.40 2.63 2.27 18 2 setar 2 1 3 1.80 1.83 1.84 1.84 > > models_univariate %>% + mutate(boot = map(object, ~ tibble(n = 1:3, boot =head(setar.boot(., seed = 123, returnStarting = TRUE), 3)))) %>% + select(-object) %>% + unnest(boot) %>% + spread(include, boot)%>% + as.data.frame() %>% + print(digits=3) lag model nthresh thDelay n both const none trend 1 1 setar 1 0 1 2.40 2.40 2.40 2.40 2 1 setar 1 0 2 2.68 3.26 3.27 2.58 3 1 setar 1 0 3 2.81 3.32 3.46 2.95 4 1 setar 2 0 1 2.40 2.40 2.40 2.40 5 1 setar 2 0 2 3.26 3.31 3.30 3.26 6 1 setar 2 0 3 3.41 3.21 3.50 3.79 7 2 setar 1 0 1 2.40 2.40 2.40 2.40 8 2 setar 1 0 2 2.40 2.40 2.40 2.40 9 2 setar 1 0 3 2.25 2.29 2.25 2.15 10 2 setar 1 1 1 2.40 2.40 2.40 2.40 11 2 setar 1 1 2 2.40 2.40 2.40 2.40 12 2 setar 1 1 3 2.04 2.46 2.19 2.28 13 2 setar 2 0 1 2.40 2.40 2.40 2.40 14 2 setar 2 0 2 2.40 2.40 2.40 2.40 15 2 setar 2 0 3 2.26 2.12 2.13 2.20 16 2 setar 2 1 1 2.40 2.40 2.40 2.40 17 2 setar 2 1 2 2.40 2.40 2.40 2.40 18 2 setar 2 1 3 1.99 2.20 2.21 2.06 > > ## add regime > models_univariate %>% + mutate(boot = map(object, ~ setar.boot(., seed = 123, add.regime = TRUE) %>% + head(3) %>% + as_tibble() %>% + mutate(n = 1:3))) %>% + select(-object) %>% + unnest(boot) %>% + gather(stat, value, res, regime) %>% + unite(stat_n, stat, n) %>% + spread(stat_n, value) %>% + as.data.frame() %>% + print(digits=3) lag include model nthresh thDelay regime_1 regime_2 regime_3 res_1 res_2 1 1 both setar 1 0 2 2 2 2.68 2.81 2 1 both setar 2 0 3 3 3 3.26 3.41 3 1 const setar 1 0 2 2 2 3.26 3.32 4 1 const setar 2 0 2 3 3 3.31 3.21 5 1 none setar 1 0 2 2 2 3.27 3.46 6 1 none setar 2 0 3 3 3 3.30 3.50 7 1 trend setar 1 0 2 2 2 2.58 2.95 8 1 trend setar 2 0 3 3 3 3.26 3.79 9 2 both setar 1 0 2 2 2 2.25 2.44 10 2 both setar 1 1 2 2 2 2.04 2.16 11 2 both setar 2 0 3 3 3 2.26 2.44 12 2 both setar 2 1 2 2 2 1.99 2.06 13 2 const setar 1 0 2 2 2 2.29 2.35 14 2 const setar 1 1 2 2 2 2.46 2.68 15 2 const setar 2 0 3 2 3 2.12 2.53 16 2 const setar 2 1 2 2 2 2.20 2.40 17 2 none setar 1 0 2 2 2 2.25 2.45 18 2 none setar 1 1 2 2 2 2.19 2.50 19 2 none setar 2 0 3 2 3 2.13 2.52 20 2 none setar 2 1 2 2 2 2.21 2.63 21 2 trend setar 1 0 2 2 2 2.15 2.47 22 2 trend setar 1 1 2 2 2 2.28 2.23 23 2 trend setar 2 0 3 3 3 2.20 2.70 24 2 trend setar 2 1 2 2 2 2.06 2.27 res_3 1 2.15 2 3.09 3 2.71 4 2.80 5 2.83 6 2.90 7 2.60 8 3.22 9 1.96 10 1.86 11 2.12 12 1.80 13 1.87 14 2.34 15 2.22 16 1.83 17 1.92 18 2.14 19 2.29 20 1.84 21 2.25 22 1.59 23 2.15 24 1.84 > > ## add regime and starting > models_univariate %>% + mutate(boot = map(object, ~ setar.boot(., seed = 123, add.regime = TRUE, returnStarting = TRUE) %>% + head(3) %>% + as_tibble() %>% + mutate(n = 1:3))) %>% + select(-object) %>% + unnest(boot) %>% + gather(stat, value, res, regime) %>% + unite(stat_n, stat, n) %>% + spread(stat_n, value) %>% + as.data.frame() %>% + print(digits=3) lag include model nthresh thDelay regime_1 regime_2 regime_3 res_1 res_2 1 1 both setar 1 0 NA 2 2 2.4 2.68 2 1 both setar 2 0 NA 3 3 2.4 3.26 3 1 const setar 1 0 NA 2 2 2.4 3.26 4 1 const setar 2 0 NA 2 3 2.4 3.31 5 1 none setar 1 0 NA 2 2 2.4 3.27 6 1 none setar 2 0 NA 3 3 2.4 3.30 7 1 trend setar 1 0 NA 2 2 2.4 2.58 8 1 trend setar 2 0 NA 3 3 2.4 3.26 9 2 both setar 1 0 NA NA 2 2.4 2.40 10 2 both setar 1 1 NA NA 2 2.4 2.40 11 2 both setar 2 0 NA NA 3 2.4 2.40 12 2 both setar 2 1 NA NA 2 2.4 2.40 13 2 const setar 1 0 NA NA 2 2.4 2.40 14 2 const setar 1 1 NA NA 2 2.4 2.40 15 2 const setar 2 0 NA NA 3 2.4 2.40 16 2 const setar 2 1 NA NA 2 2.4 2.40 17 2 none setar 1 0 NA NA 2 2.4 2.40 18 2 none setar 1 1 NA NA 2 2.4 2.40 19 2 none setar 2 0 NA NA 3 2.4 2.40 20 2 none setar 2 1 NA NA 2 2.4 2.40 21 2 trend setar 1 0 NA NA 2 2.4 2.40 22 2 trend setar 1 1 NA NA 2 2.4 2.40 23 2 trend setar 2 0 NA NA 3 2.4 2.40 24 2 trend setar 2 1 NA NA 2 2.4 2.40 res_3 1 2.81 2 3.41 3 3.32 4 3.21 5 3.46 6 3.50 7 2.95 8 3.79 9 2.25 10 2.04 11 2.26 12 1.99 13 2.29 14 2.46 15 2.12 16 2.20 17 2.25 18 2.19 19 2.13 20 2.21 21 2.15 22 2.28 23 2.20 24 2.06 > > ################ > ### test boot > ################ > > setar.boot.check <- function(object, round_digits = 10, ...) { + mod_boot <- setar.boot(object, boot.scheme = "check", round_digits = round_digits, returnStarting = TRUE) + orig_series <- as.numeric(object$str$x) + all.equal(mod_boot, orig_series, ...) + } > > linear.boot.check <- function(object, round_digits = 10) { + mod_boot <- linear.boot(object, boot.scheme = "check",round_digits = round_digits, returnStarting = TRUE) + orig_series <- as.numeric(object$str$x) + all.equal(mod_boot, orig_series) + } > > > > ## nthresh ==0 > ar_1_noInc <- linear(log(lynx), m = 1, include = "none") > ar_2_noInc <- linear(log(lynx), m = 2, include = "none") > ar_1_const <- linear(log(lynx), m = 1, include = "const") > ar_2_const <- linear(log(lynx), m = 2, include = "const") > > > setar.boot.check(object=ar_1_noInc) [1] TRUE > setar.boot.check(ar_2_noInc) [1] TRUE > setar.boot.check(ar_1_const) [1] TRUE > setar.boot.check(ar_2_const) [1] TRUE > > linear.boot.check(ar_1_noInc) [1] TRUE > linear.boot.check(ar_2_noInc) [1] TRUE > linear.boot.check(ar_1_const) [1] TRUE > linear.boot.check(ar_2_const) [1] TRUE > > > > ## nthresh ==1 > set_1th_l1 <- setar(lynx, nthresh=1, m=1) Warning message: Possible unit root in the low regime. Roots are: 0.5005 > set_1th_l2 <- setar(lynx, nthresh=1, m=2) Warning message: Possible unit root in the low regime. Roots are: 0.4681 5.0788 > set_1th_l1_tr <- setar(lynx, nthresh=1, m=1, include = "trend") Warning message: Possible unit root in the low regime. Roots are: 0.633 > > > > ## lot of problems with numerical instability on other platforms, unreliable tests > if(FALSE){ + setar.boot.check(object=set_1th_l1, round_digits = 10, tol=1e-6) + setar.boot.check(set_1th_l1, round_digits = 2) + setar.boot.check(set_1th_l2, tol=1e-1) + setar.boot.check(set_1th_l2, round_digits = 5, tol=1e-1) + setar.boot.check(set_1th_l1_tr, round_digits = 1) + } > > > > ## why difference? > if(FALSE) { + library(ggplot) + getTh(set_1th_l2) + filt_diff <- function(x, minus=2, tol =1) { + x2 <- x %>% + mutate(diff = x$boot - x$orig) + first <- which(abs(x2$diff)>tol)[1] + filter(x2, between(n_row, first -minus, first +minus)) + } + set_1th_l2_b <- setar.boot(setarObject = set_1th_l2, boot.scheme = "check", round_digits = 7) + + df_comp <- tibble(orig = lynx, boot = set_1th_l2_b) %>% + mutate(n_row = 1:n(), + th1 = getTh(set_1th_l1_tr)[1], + th2 = getTh(set_1th_l1_tr)[2], + reg = regime(set_1th_l1_tr)) + + df_comp %>% + filt_diff(tol = 0.01) + df_comp %>% + qplot(x=n_row, y = as.numeric(orig), data =., geom = "line") + + geom_point(aes(colour = as.numeric(reg) %>% factor)) + geom_line(aes(y = boot), colour = "red") + } > > ## nthresh == 2 > > ### boot > set_2th_l1 <- setar(lynx, nthresh=2, m=1) Warning message: Possible unit root in the low regime. Roots are: 0.7924 > set_2th_l2 <- setar(lynx, nthresh=2, m=2) Warning message: Possible unit root in the low regime. Roots are: 0.5739 4.6636 > set_2th_l1_tr <- setar(lynx, nthresh=2, m=1, include = "trend") Warning message: Possible unit root in the low regime. Roots are: 0.633 > > > setar.boot.check(set_2th_l1) [1] TRUE > # setar.boot.check(set_2th_l2, round_digits = 2, countEQ=TRUE, scale=1) ## big diffs > # isTRUE(setar.boot.check(set_2th_l1_tr)) ## explsive! > # setar.boot.check(set_2th_l1_tr, round_digits = 2, tol=0.0001) # explsoive too! > > ################ > ### tets sim > ################ > > ## nthresh ==0 > set.seed(123) > innov_1 <- rnorm(200) > sim_nth0 <- setar.sim(B=0.5, lag=1, nthresh=0, + include ="none", + starting= 0.4, + innov=innov_1, + show.parMat = TRUE) const trend phi_1 [1,] 0 0 0.5 > > head(sim_nth0) [1] -0.3604756 -0.4104153 1.3535007 0.7472587 0.5029171 1.9665235 > > ## nthresh ==1 > Bvals <- c(2.9,-0.4,-0.1, 2, 0.2,0.3) > sim_new <- setar.sim(B=Bvals,lag=2, nthresh=1, Thresh=2, starting=c(2.8,2.2), + innov=innov_1, show.parMat = TRUE) const_L trend_L phi_1_L phi_2_L const_H trend_H phi_1_H phi_2_H [1,] 2.9 0 -0.4 -0.1 2 0 0.2 0.3 > > head(sim_new) [1] 2.719524 2.973727 4.969311 3.956489 4.411379 5.784287 > > proc.time() user system elapsed 2.85 0.48 3.31