Nothing
#' Break Detection in the Seasonal and Trend Component of a Univariate Time
#' Series
#'
#' Iterative break detection in seasonal and trend component of a time series.
#' Seasonal breaks is a function that combines the iterative decomposition of
#' time series into trend, seasonal and remainder components with significant
#' break detection in the decomposed components of the time series.
#'
#' The algorithm decomposes the input time series `Yt` into three components:
#' trend, seasonality and remainder, using the function defined by the `decomp`
#' parameter. Then each component is checked for at least one significant
#' break using [strucchangeRcpp::efp()], and if there is one, [strucchangeRcpp::breakpoints()]
#' is run on the component. The result allows differentiating between breaks in
#' trend and seasonality.
#'
#' @param Yt univariate time series to be analyzed. This should be an object of
#' class "ts" with a frequency greater than one.
#' @param h minimal segment size between potentially detected breaks in the
#' trend model given as fraction relative to the sample size (i.e. the minimal
#' number of observations in each segment divided by the total length of the
#' timeseries).
#' @param season the seasonal model used to fit the seasonal component and
#' detect seasonal breaks (i.e. significant phenological change). There are
#' three options: "dummy", "harmonic", or "none" where "dummy" is the model
#' proposed in the first Remote Sensing of Environment paper and "harmonic" is
#' the model used in the second Remote Sensing of Environment paper (See paper
#' for more details) and where "none" indicates that no seasonal model will be
#' fitted (i.e. St = 0 ). If there is no seasonal cycle (e.g. frequency of the
#' 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 hpc A character specifying the high performance computing support.
#' Default is "none", can be set to "foreach". Install the "foreach" package
#' for hpc support.
#' @param level numeric; threshold value for the \link[strucchangeRcpp]{sctest.efp}
#' test; if a length 2 vector is passed, the first value is used for the trend,
#' the second for the seasonality
#' @param decomp "stlplus" or "stl": the function to use for decomposition.
#' \code{stl} can handle sparse time series (1 < frequency < 4), \code{stlplus}
#' can handle \code{NA} values in the time series.
#' @param type character, indicating the type argument to
#' \link[strucchangeRcpp]{efp}
#' @param \dots additional arguments passed to [stlplus::stlplus()], if
#' its usage has been enabled, otherwise ignored.
#' @return An object of the class "bfast" is a list with the following
#' elements: \item{Yt}{ equals the Yt used as input.} \item{output}{ is a list
#' with the following elements (for each iteration): \tabular{ll}{ \code{Tt}
#' \tab the fitted trend component\cr \code{St} \tab the fitted seasonal
#' component\cr \code{Nt} \tab the noise or remainder component\cr \code{Vt}
#' \tab equals the deseasonalized data \code{Yt - St} for each iteration\cr
#' \code{bp.Vt} \tab output of the \code{\link[strucchangeRcpp]{breakpoints}}
#' function for the trend model. Note that the output breakpoints are index
#' numbers of \code{na.omit(as.numeric(Vt))}.\cr \code{ci.Vt} \tab output of the
#' \code{\link[strucchangeRcpp]{breakpoints}} confint function for the trend
#' model\cr \code{Wt} \tab equals the detrended data \code{Yt - Tt} for each
#' iteration\cr \code{bp.Wt} \tab output of the
#' \code{\link[strucchangeRcpp]{breakpoints}} function for the seasonal model.
#' Note that the output breakpoints are index numbers of \code{na.omit(as.numeric(Wt))}. \cr
#' \code{ci.Wt} \tab output of the \code{\link[strucchangeRcpp]{breakpoints}}
#' confint function for the seasonal model }} \item{nobp}{ is a list with the
#' following elements: \tabular{ll}{ \code{nobp.Vt} \tab logical, TRUE if there
#' are no breakpoints detected\cr \code{nobp.Wt} \tab logical, TRUE if there
#' are no breakpoints detected }} \item{Magnitude}{ magnitude of the biggest
#' change detected in the trend component} \item{Time}{ timing of the biggest
#' change detected in the trend component}
#' @author Jan Verbesselt
#' @seealso \code{\link[bfast]{plot.bfast}} for plotting of bfast() results.
#' \cr \code{\link[strucchangeRcpp]{breakpoints}} for more examples and background
#' information about estimation of breakpoints in time series.
#' @references \insertRef{janbfast}{bfast}
#'
#' \insertRef{janbfastpheno}{bfast}
#' @keywords ts
#' @example examples/bfast.r
#'
#' @export bfast
bfast <- function (Yt, h = 0.15, season = c("dummy", "harmonic", "none"),
max.iter = 10, breaks = NULL, hpc = "none", level = 0.05,
decomp = c("stl", "stlplus"),
type = "OLS-MOSUM", ...)
{
# Error catching
decomp = match.arg(decomp)
if (decomp == "stl" && any(is.na(Yt))) {
warning("The stl() function cannot deal with missing values in the time series, falling back to decomp='stlplus'.")
decomp <- "stlplus"
}
if(decomp == "stlplus" && !requireNamespace("stlplus", quietly = TRUE)) {
warning("stlplus package could not be loaded, falling back to decomp='stl'.")
decomp <- "stl"
}
## Get Arguments
season <- match.arg(season)
level <- rep(level, length.out = 2)
ti <- time(Yt)
f <- frequency(Yt) # on cycle every f time points (seasonal cycle)
if (class(Yt) != "ts")
stop("Not a time series object")
if (f < 1)
stop("Time series frequency needs to be more than 1")
if (decomp == "stlplus" && f < 4) {
warning("Not enough seasons (frequency < 4) for using stlplus decomposition, falling back to decomp = 'stl'")
decomp <- "stl"
}
output <- list()
Tt <- 0
if (season == "harmonic") {
w <- 1/f
tl <- 1:length(Yt)
co <- cos(2 * pi * tl * w)
si <- sin(2 * pi * tl * w)
co2 <- cos(2 * pi * tl * w * 2)
si2 <- sin(2 * pi * tl * w * 2)
co3 <- cos(2 * pi * tl * w * 3)
si3 <- sin(2 * pi * tl * w * 3)
smod <- Wt ~ co + si + co2 + si2 + co3 + si3
# Start the iterative procedure and for first iteration St=decompose result
St <- if (decomp == "stlplus") {
stlplus::stlplus(Yt, t=ti, n.p = f, s.window = "periodic", ...)$data[, "seasonal"]
} else {
stl (Yt, "periodic")$time.series[, "seasonal"]
}
}
else if (season == "dummy") {
# Start the iterative procedure and for first iteration St=decompose result
St <- if (decomp == "stlplus") {
stlplus::stlplus(Yt, t = ti, n.p = f, s.window = "periodic", ...)$data[, "seasonal"]
} else {
stl(Yt, "periodic")$time.series[, "seasonal"]
}
D <- seasonaldummy(Yt)
D[rowSums(D) == 0, ] <- -1
smod <- Wt ~ -1 + D
}
else if (season == "none") {
print("No seasonal model will be fitted!")
St <- 0
}
else stop("Not a correct seasonal model is selected ('harmonic' or 'dummy') ")
Vt.bp <- 0
Wt.bp <- 0
CheckTimeTt <- 1
CheckTimeSt <- 1
i <- 0
while ((!identical(CheckTimeTt, Vt.bp) | !identical(CheckTimeSt, Wt.bp)) & i < max.iter) {
CheckTimeTt <- Vt.bp
CheckTimeSt <- Wt.bp
### Change in trend component
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)[1])
} else {
nobp.Vt <- TRUE
bp.Vt <- NA
}
if (nobp.Vt) {
## No Change detected
fm0 <- lm(Vt ~ ti)
Vt.bp <- 0 # no breaks times
Tt <- ts(data=NA,start = ti[1], end = ti[length(ti)],frequency = f) # Data minus trend
Tt[which(!is.na(Yt))] <- fitted(fm0) # Overwrite non-missing with fitted values
tsp(Tt) <- tsp(Yt)
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)
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)
Tt[which(!is.na(Yt))] <- fitted(fm1) # Overwrite non-missing with fitted values
tsp(Tt) <- tsp(Yt)
}
if (season == "none") {
Wt <- 0
St <- 0
bp.Wt <- NA
ci.Wt <- NA
nobp.Wt <- TRUE
} else {
### Change in seasonal component
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,
hpc = hpc)
nobp.Wt <- is.na(breakpoints(bp.Wt)[1])
}
else {
nobp.Wt <- TRUE
bp.Wt <- NA
}
if (nobp.Wt) {
## No seasonal change detected
sm0 <- lm(smod)
St <- ts(data=NA,start = ti[1], end = ti[length(ti)],frequency = f)
St[which(!is.na(Yt))] <- fitted(sm0) # Overwrite non-missing with fitted values
tsp(St) <- tsp(Yt)
Wt.bp <- 0
ci.Wt <- NA
} else {
if (season == "dummy")
sm1 <- lm(Wt[!is.na(Wt)] ~ -1 + D[!is.na(Wt),] %in% breakfactor(bp.Wt))
if (season == "harmonic")
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)
Wt.bp <- ci.Wt$confint[, 2]
# Define empty copy of original time series
St <- ts(data=NA,start = ti[1], end = ti[length(ti)],frequency = f)
St[which(!is.na(Yt))] <- fitted(sm1) # Overwrite non-missing with fitted values
tsp(St) <- tsp(Yt)
}
}
i <- i + 1
output[[i]] <- list(Tt = Tt, St = St, Nt = Yt - Tt - St, Vt = Vt, bp.Vt = bp.Vt, Vt.bp = Vt.bp, ci.Vt = ci.Vt,
Wt = Wt, bp.Wt = bp.Wt, Wt.bp = Wt.bp, ci.Wt = ci.Wt)
}
if (!nobp.Vt) { # probably only works well for dummy model!
Vt.nrbp <- length(bp.Vt$breakpoints)
co <- coef(fm1) # final fitted trend model
Mag <- matrix(NA, Vt.nrbp, 3)
for (r in 1:Vt.nrbp) {
if (r == 1) {
y1 <- co[1] + co[r + Vt.nrbp + 1] * ti[Vt.bp[r]]
} else {
y1 <- co[1] + co[r] + co[r + Vt.nrbp + 1] * ti[Vt.bp[r]]
}
y2 <- (co[1] + co[r + 1]) + co[r + Vt.nrbp + 2] * ti[Vt.bp[r] + 1]
Mag[r, 1] <- y1
Mag[r, 2] <- y2
Mag[r, 3] <- y2 - y1
}
index <- which.max(abs(Mag[, 3]))
m.x <- rep(Vt.bp[index], 2)
m.y <- c(Mag[index, 1], Mag[index, 2]) #Magnitude position
Magnitude <- Mag[index, 3] # Magnitude of biggest change
Time <- Vt.bp[index]
} else {
m.x <- NA
m.y <- NA
Magnitude <- 0 # if we do not detect a break then the magnitude is zero
Time <- NA # if we do not detect a break then we have no timing of the break
Mag <- 0
}
return(structure(list(Yt = Yt, output = output, nobp = list(Vt = nobp.Vt, Wt = nobp.Wt),
Magnitude = Magnitude, Mags = Mag, Time = Time,
jump = list(x = ti[m.x], y = m.y)), class = "bfast"))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.