From 279892a422d9dc5c57597eab52fcd09977cca3c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dainius=20Masili=C5=ABnas?= Date: Wed, 22 Dec 2021 19:12:35 +0100 Subject: [PATCH 1/8] Fix not honouring numeric breaks in confint bfast() will now also calculate the correct number of breaks when it is requested. Links to the fixes also in strucchangeRcpp, see https://github.com/bfast2/strucchangeRcpp/issues/40. --- R/bfast.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/bfast.R b/R/bfast.R index 941e061..a9a10c1 100644 --- a/R/bfast.R +++ b/R/bfast.R @@ -167,7 +167,7 @@ bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"), ci.Vt <- NA } else { fm1 <- lm(Vt[which(!is.na(Yt))] ~ breakfactor(bp.Vt)/ti[which(!is.na(Yt))] ) - ci.Vt <- confint(bp.Vt, het.err = FALSE) + ci.Vt <- confint(bp.Vt, het.err = FALSE, breaks = breaks) Vt.bp <- ci.Vt$confint[, 2] # Define empty copy of original time series Tt <- ts(data=NA,start = ti[1], end = ti[length(ti)],frequency = f) @@ -208,7 +208,7 @@ bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"), sm1 <- lm(Wt[!is.na(Wt)] ~ ( co[!is.na(Wt)] + si[!is.na(Wt)] + co2[!is.na(Wt)] + si2[!is.na(Wt)] + co3[!is.na(Wt)] + si3[!is.na(Wt)] ) %in% breakfactor(bp.Wt)) - ci.Wt <- confint(bp.Wt, het.err = FALSE) + ci.Wt <- confint(bp.Wt, het.err = FALSE, breaks = breaks) Wt.bp <- ci.Wt$confint[, 2] # Define empty copy of original time series From 1753daccdbada1c0ccd29e1567794a0ea2a93de4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dainius=20Masili=C5=ABnas?= Date: Wed, 22 Dec 2021 19:13:11 +0100 Subject: [PATCH 2/8] Bump version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index bed4c69..a821647 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,5 +1,5 @@ Package: bfast -Version: 1.7.0 +Version: 1.7.0.9999 Title: Breaks for Additive Season and Trend Authors@R: c(person(given = "Jan", family = "Verbesselt", role = c("aut"), email = "Jan.Verbesselt@wur.nl"), person(given = "Dainius", family = "Masiliunas", role = c("aut", "cre"), From 4d953b93fd637888fec315bebc23ef26d7cbae61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dainius=20Masili=C5=ABnas?= Date: Wed, 22 Dec 2021 19:15:42 +0100 Subject: [PATCH 3/8] Rebump version, 1.7.0 is not yet upstreamed --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index a821647..f2486e5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,5 +1,5 @@ Package: bfast -Version: 1.7.0.9999 +Version: 1.6.1.9999 Title: Breaks for Additive Season and Trend Authors@R: c(person(given = "Jan", family = "Verbesselt", role = c("aut"), email = "Jan.Verbesselt@wur.nl"), person(given = "Dainius", family = "Masiliunas", role = c("aut", "cre"), From 1bd4416ba68e500dd8e9b6722387f0b00635d9e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dainius=20Masili=C5=ABnas?= Date: Wed, 22 Dec 2021 19:16:32 +0100 Subject: [PATCH 4/8] Update NEWS --- NEWS | 3 +++ 1 file changed, 3 insertions(+) diff --git a/NEWS b/NEWS index 2c8eae3..bbe52d3 100644 --- a/NEWS +++ b/NEWS @@ -5,6 +5,9 @@ Changes in Version 1.7.0 anything equal to or below 0 in 'level' will skip the sctest. o Fixed a bug where bfastpp() would throw an error when annual data is input. + + o bfast() now honours a numeric number of 'breaks' also in calculating + confidence intervals. Changes in Version 1.6.1 From 02ad305b232fe0e90cb1c7457b8cf76d3fc8a8e9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dainius=20Masili=C5=ABnas?= Date: Wed, 22 Dec 2021 19:30:02 +0100 Subject: [PATCH 5/8] Honour breaks in checking if there are any breakpoints Also update the description of the breaks parameter, as it represents something different than currently stated --- R/bfast.R | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/R/bfast.R b/R/bfast.R index a9a10c1..a3005e8 100644 --- a/R/bfast.R +++ b/R/bfast.R @@ -29,8 +29,10 @@ #' time series is 1) "none" can be selected to avoid fitting a seasonal model. #' @param max.iter maximum amount of iterations allowed for estimation of #' breakpoints in seasonal and trend component. -#' @param breaks integer specifying the maximal number of breaks to be -#' calculated. By default the maximal number allowed by h is used. +#' @param breaks either an integer specifying the number of breaks to compute +#' for both components, or a string specifying a statistic with which to compute +#' the optimal number of breakpoints (see [strucchangeRcpp::breakpoints()] for +#' more information) #' @param hpc A character specifying the high performance computing support. #' Default is "none", can be set to "foreach". Install the "foreach" package #' for hpc support. @@ -152,7 +154,7 @@ bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"), p.Vt <- sctest(efp(Vt ~ ti, h = h, type = type)) if (p.Vt$p.value <= level[1]) { bp.Vt <- breakpoints(Vt ~ ti, h = h, breaks = breaks, na.action=na.exclude,hpc = hpc) - nobp.Vt <- is.na(breakpoints(bp.Vt)[1]) + nobp.Vt <- is.na(breakpoints(bp.Vt, breaks = breaks)[1]) } else { nobp.Vt <- TRUE bp.Vt <- NA @@ -187,7 +189,7 @@ bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"), if (p.Wt$p.value <= level[2]) { bp.Wt <- breakpoints(smod, h = h, breaks = breaks, hpc = hpc) - nobp.Wt <- is.na(breakpoints(bp.Wt)[1]) + nobp.Wt <- is.na(breakpoints(bp.Wt, breaks = breaks)[1]) } else { nobp.Wt <- TRUE From e2d9411d70c9ab123be0b1cf12a2e9b1c7582ab5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dainius=20Masili=C5=ABnas?= Date: Wed, 22 Dec 2021 19:30:50 +0100 Subject: [PATCH 6/8] Update NEWS --- NEWS | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS b/NEWS index bbe52d3..25dedc4 100644 --- a/NEWS +++ b/NEWS @@ -7,7 +7,7 @@ Changes in Version 1.7.0 o Fixed a bug where bfastpp() would throw an error when annual data is input. o bfast() now honours a numeric number of 'breaks' also in calculating - confidence intervals. + confidence intervals and checking if any breaks have been detected. Changes in Version 1.6.1 From 42f9ef166003a12b83e1c74b13d0cf8b10ca2fc8 Mon Sep 17 00:00:00 2001 From: Dainius Date: Thu, 23 Dec 2021 13:57:01 +0100 Subject: [PATCH 7/8] Add the ability to pass breaks to individual components in bfast() --- R/bfast.R | 28 +++++++++++++++++++--------- 1 file changed, 19 insertions(+), 9 deletions(-) diff --git a/R/bfast.R b/R/bfast.R index a3005e8..5bbe81e 100644 --- a/R/bfast.R +++ b/R/bfast.R @@ -29,10 +29,12 @@ #' time series is 1) "none" can be selected to avoid fitting a seasonal model. #' @param max.iter maximum amount of iterations allowed for estimation of #' breakpoints in seasonal and trend component. -#' @param breaks either an integer specifying the number of breaks to compute -#' for both components, or a string specifying a statistic with which to compute +#' @param breaks either an integer specifying the number of breaks to compute, +#' or a string specifying a statistic with which to compute #' the optimal number of breakpoints (see [strucchangeRcpp::breakpoints()] for -#' more information) +#' more information). If one value is given, it is used for both the trend and +#' seasonal components, and if two are given, the first one is considered to be +#' for the trend and the second for the seasonal component. #' @param hpc A character specifying the high performance computing support. #' Default is "none", can be set to "foreach". Install the "foreach" package #' for hpc support. @@ -95,6 +97,14 @@ bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"), } ## Get Arguments season <- match.arg(season) + if (length(breaks) > 1) + { + breaks.Vt <- breaks[[1]] + breaks.Wt <- breaks[[2]] + } else { + breaks.Vt <- breaks + breaks.Wt <- breaks + } level <- rep(level, length.out = 2) ti <- time(Yt) f <- frequency(Yt) # on cycle every f time points (seasonal cycle) @@ -153,8 +163,8 @@ bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"), Vt <- Yt - St # Deasonalized Time series p.Vt <- sctest(efp(Vt ~ ti, h = h, type = type)) if (p.Vt$p.value <= level[1]) { - bp.Vt <- breakpoints(Vt ~ ti, h = h, breaks = breaks, na.action=na.exclude,hpc = hpc) - nobp.Vt <- is.na(breakpoints(bp.Vt, breaks = breaks)[1]) + bp.Vt <- breakpoints(Vt ~ ti, h = h, breaks = breaks.Vt, na.action=na.exclude,hpc = hpc) + nobp.Vt <- is.na(breakpoints(bp.Vt, breaks = breaks.Vt)[1]) } else { nobp.Vt <- TRUE bp.Vt <- NA @@ -169,7 +179,7 @@ bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"), ci.Vt <- NA } else { fm1 <- lm(Vt[which(!is.na(Yt))] ~ breakfactor(bp.Vt)/ti[which(!is.na(Yt))] ) - ci.Vt <- confint(bp.Vt, het.err = FALSE, breaks = breaks) + ci.Vt <- confint(bp.Vt, het.err = FALSE, breaks = breaks.Vt) Vt.bp <- ci.Vt$confint[, 2] # Define empty copy of original time series Tt <- ts(data=NA,start = ti[1], end = ti[length(ti)],frequency = f) @@ -187,9 +197,9 @@ bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"), Wt <- Yt - Tt p.Wt <- sctest(efp(smod, h = h, type = type)) # preliminary test if (p.Wt$p.value <= level[2]) { - bp.Wt <- breakpoints(smod, h = h, breaks = breaks, + bp.Wt <- breakpoints(smod, h = h, breaks = breaks.Wt, hpc = hpc) - nobp.Wt <- is.na(breakpoints(bp.Wt, breaks = breaks)[1]) + nobp.Wt <- is.na(breakpoints(bp.Wt, breaks = breaks.Wt)[1]) } else { nobp.Wt <- TRUE @@ -210,7 +220,7 @@ bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"), sm1 <- lm(Wt[!is.na(Wt)] ~ ( co[!is.na(Wt)] + si[!is.na(Wt)] + co2[!is.na(Wt)] + si2[!is.na(Wt)] + co3[!is.na(Wt)] + si3[!is.na(Wt)] ) %in% breakfactor(bp.Wt)) - ci.Wt <- confint(bp.Wt, het.err = FALSE, breaks = breaks) + ci.Wt <- confint(bp.Wt, het.err = FALSE, breaks = breaks.Wt) Wt.bp <- ci.Wt$confint[, 2] # Define empty copy of original time series From d387e24a0afac78830fdc6049c1a14294cac6068 Mon Sep 17 00:00:00 2001 From: Dainius Date: Thu, 23 Dec 2021 14:01:11 +0100 Subject: [PATCH 8/8] Update NEWS item --- NEWS | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/NEWS b/NEWS index 25dedc4..d98fa24 100644 --- a/NEWS +++ b/NEWS @@ -6,8 +6,9 @@ Changes in Version 1.7.0 o Fixed a bug where bfastpp() would throw an error when annual data is input. - o bfast() now honours a numeric number of 'breaks' also in calculating - confidence intervals and checking if any breaks have been detected. + o Rework of the 'breaks' argument in bfast(), it now works correctly with a + numeric input and gives the opportunity to specify the breaks/statistic for + computing the optimal number of breaks separately per component. Changes in Version 1.6.1