bfast {bfast}

Break Detection in the Seasonal and Trend Component of a Univariate Time Series

Description

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.

Usage

bfast(
  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",
  ...
)

Arguments

Yt

univariate time series to be analyzed. This should be an object of class "ts" with a frequency greater than one.

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).

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.

max.iter

maximum amount of iterations allowed for estimation of breakpoints in seasonal and trend component.

breaks

integer specifying the maximal number of breaks to be calculated. By default the maximal number allowed by h is used.

hpc

A character specifying the high performance computing support. Default is "none", can be set to "foreach". Install the "foreach" package for hpc support.

level

numeric; threshold value for the sctest.efp test; if a length 2 vector is passed, the first value is used for the trend, the second for the seasonality

decomp

"stlplus" or "stl": the function to use for decomposition. stl can handle sparse time series (1 < frequency < 4), stlplus can handle NA values in the time series.

type

character, indicating the type argument to efp

...

additional arguments passed to stlplus::stlplus(), if its usage has been enabled, otherwise ignored.

Details

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.

Value

An object of the class "bfast" is a list with the following elements:

Yt

equals the Yt used as input.

output

is a list with the following elements (for each iteration):

Tt	the fitted trend component
St	the fitted seasonal component
Nt	the noise or remainder component
Vt	equals the deseasonalized data Yt - St for each iteration
bp.Vt	output of the breakpoints function for the trend model. Note that the output breakpoints are index numbers of na.omit(as.numeric(Vt)).
ci.Vt	output of the breakpoints confint function for the trend model
Wt	equals the detrended data Yt - Tt for each iteration
bp.Wt	output of the breakpoints function for the seasonal model. Note that the output breakpoints are index numbers of na.omit(as.numeric(Wt)).
ci.Wt	output of the breakpoints confint function for the seasonal model

nobp

is a list with the following elements:

nobp.Vt	logical, TRUE if there are no breakpoints detected
nobp.Wt	logical, TRUE if there are no breakpoints detected

Magnitude

magnitude of the biggest change detected in the trend component

Time

timing of the biggest change detected in the trend component

Author(s)

Jan Verbesselt

References

Verbesselt J, Hyndman R, Newnham G, Culvenor D (2010). “Detecting trend and seasonal changes in satellite image time series.” Remote Sensing of Environment, 114(1), 106–115. ISSN 0034-4257, doi: 10.1016/j.rse.2009.08.014, https://doi.org/10.1016/j.rse.2009.08.014.

Verbesselt J, Hyndman R, Zeileis A, Culvenor D (2010). “Phenological change detection while accounting for abrupt and gradual trends in satellite image time series.” Remote Sensing of Environment, 114(12), 2970–2980. ISSN 0034-4257, doi: 10.1016/j.rse.2010.08.003, https://doi.org/10.1016/j.rse.2010.08.003.

See Also

plot.bfast for plotting of bfast() results.
breakpoints for more examples and background information about estimation of breakpoints in time series.

Examples

## Simulated Data
plot(simts) # stl object containing simulated NDVI time series
datats <- ts(rowSums(simts$time.series))
# sum of all the components (season,abrupt,remainder)
tsp(datats) <- tsp(simts$time.series) # assign correct time series attributes
plot(datats)

fit <- bfast(datats, h = 0.15, season = "dummy", max.iter = 1)
plot(fit, sim = simts)
fit
# prints out whether breakpoints are detected
# in the seasonal and trend component
 
## Real data
## The data should be a regular ts() object without NA's
## See Fig. 8 b in reference
plot(harvest, ylab = "NDVI")
# MODIS 16-day cleaned and interpolated NDVI time series

(rdist <- 10/length(harvest))
# ratio of distance between breaks (time steps) and length of the time series

fit <- bfast(harvest, h = rdist, season = "harmonic", max.iter = 1, breaks = 2)
plot(fit)
## plot anova and slope of the trend identified trend segments
plot(fit, ANOVA = TRUE)
## plot the trend component and identify the break with
## the largest magnitude of change
plot(fit, type = "trend", largest = TRUE)

## plot all the different available plots
plot(fit, type = "all")

## output
niter <- length(fit$output) # nr of iterations
out <- fit$output[[niter]]
# output of results of the final fitted seasonal and trend models and
## #nr of breakpoints in both.

## running bfast on yearly data
t <- ts(as.numeric(harvest), frequency = 1, start = 2006)
fit <- bfast(t, h = 0.23, season = "none", max.iter = 1)
plot(fit)
fit

## handling missing values with stlplus
(NDVIa <- as.ts(zoo::zoo(som$NDVI.a, som$Time)))
fit <- bfast(NDVIa, season = "harmonic", max.iter = 1, decomp = "stlplus")
plot(fit)
fit