# NOTE: THIS IS CARRIED OVER FROM stlboot.R
# CLEARLY FLAG ANY CHANGES!!!!!!!!!!!!!!!!!!!!!!!
#____________________________________
## stlboot {TStools} R Documentation
#____________________________________
## STL/Loess bootstrapping
##
## Description
##
## STL/Loess bootstrapping by Bergmeir, Hyndman and Benitez, 2014.
##
## Usage
## stlboot(ts,k=1,test.season=c(TRUE,FALSE),outplot=c(FALSE,TRUE))
##
##
## Arguments
##
## ts: Time series to bootstrap. Must be ts object.
## k: Number of bootstraps.
## test.season: If TRUE then test for presence of seasonality.
## If FALSE then all non-yearly are decomposed using STL, assuming seasonality.
## outplot: If TRUE provide a plot of the bootstrapped series.
##
##
## Value
##
## ts.recon: Array of bootstrapped series. Each column is a time series.
##
## Author(s): Nikolaos Kourentzes
##
## References
##
## Bergmeir C., Hyndman R. J., Benitez J. M., Bagging Exponential Smoothing Methods using STL Decomposition and Box-Cox Transformation. 2014, Working Paper. http://robjhyndman.com/working-papers/bagging-ets/
##
## See Also
## stl, loess.
##
## Examples
##
## stlboot(referrals,k=20,outplot=TRUE)
##
#__________________________________________
stlboot <- function (ts, k = 1, test.season = c(TRUE, FALSE), outplot = c(FALSE,TRUE)){
test.season <- test.season[1]
outplot <- outplot[1]
lambda <- forecast::BoxCox.lambda(ts, method = "guerrero", lower = 0,
upper = 1)
ts.bc <- forecast::BoxCox(ts, lambda)
m <- frequency(ts)
if (m > 1) {
season.exist <- TRUE
if (test.season == TRUE) {
#mf: seasplot is from tsutils, which is a dependency of forecast.
season.exist <- seasplot(ts.bc, outplot = 0)$season.exist
}
} else {
season.exist <- FALSE
}
if (season.exist == TRUE) {
ts.decomp <- stl(ts.bc, s.window = "periodic", s.degree = 0,
t.degree = 1, l.degree = 1)
remainder <- ts.decomp$time.series[, "remainder"]
} else {
x <- 1:length(ts.bc)
ts.decomp <- loess(ts.bc ~ x, data.frame(x = x, y = ts.bc),
degree = 1, span = 6/length(ts.bc))
ts.decomp.trend <- predict(ts.decomp, data.frame(x = x))
remainder <- ts.bc - ts.decomp.trend
}
n <- length(remainder)
l <- m * 2
boot <- function(x) {
sapply(x, function(x) {
remainder[(x - l + 1):x]
})
}
boot.sample <- array(NA, c(k, n))
for (i in 1:k) {
endpoint <- sample(l:n, size = (floor(n/l) + 2))
ts.bt <- as.vector(sapply(endpoint, function(x) {
remainder[(x - l + 1):x]
}))
ts.bt <- ts.bt[-(1:(sample(l, size = 1) - 1))]
ts.bt <- ts.bt[1:n]
boot.sample[i, ] <- ts.bt
}
if (season.exist == TRUE) {
ts.bc.recon <- boot.sample + t(replicate(k, as.vector(ts.decomp$time.series[,
"trend"]))) + t(replicate(k, as.vector(ts.decomp$time.series[,
"seasonal"])))
ts.recon <- forecast::InvBoxCox(ts.bc.recon, lambda)
} else {
ts.bc.recon <- boot.sample + t(replicate(k, ts.decomp.trend))
ts.recon <- forecast::InvBoxCox(ts.bc.recon, lambda)
}
rownames(ts.recon) <- paste0("Boot", 1:k)
if (outplot == TRUE) {
plot(1:n, ts, type = "l", xlab = "Period", ylab = "",
lwd = 2)
for (i in 1:k) {
lines(1:n, ts.recon[i, ], col = "red")
}
lines(1:n, ts, col = "black", lwd = 2)
}
ts.recon <- t(ts.recon)
ts.recon <- ts(ts.recon, frequency = m, start = start(ts))
return(ts.recon)
}
#' @title Seasonal plots with simplistic trend/season tests
#'
#' @description Construct seasonal plots of various styles for a given time
#' series. The series can automatically detrended as needed.
#'
#' @param y input time series. Can be \code{ts} object.
#' @param m seasonal period. If \code{y} is a \code{ts} object then the default
#' is its frequency.
#' @param s starting period in the season. If \code{y} is a \code{ts} object
#' then this is picked up from \code{y}.
#' @param trend if \code{TRUE}, then presence of trend is assumed and removed.
#' If \code{FALSE} no trend is assumed. Use \code{NULL} to identify
#' automatically.
#' @param colour single colour override for plots.
#' @param alpha significance level for statistical tests.
#' @param outplot type of seasonal plot \itemize{ \item{0}{: none.} \item{1}{:
#' seasonal diagram.} \item{2}{: seasonal boxplots.} \item{3}{: seasonal
#' subseries.} \item{4}{: seasonal distribution.} \item{5}{: seasonal
#' density.} }
#' @param decomposition type of seasonal decomposition. This can be
#' \code{"multiplicative"}, \code{"additive"} or \code{"auto"}. If \code{y}
#' contains non-positive values then this is forced to \code{"additive"}.
#' @param cma input precalculated level/trend for the analysis. This overrides
#' \code{trend=NULL}.
#' @param labels external labels for the seasonal periods. Use \code{NULL} for
#' none. If \code{length(labels) < m}, then this input is ignored.
#' @param ... additional arguments passed to plotting functions. For example,
#' use \code{main=""} to replace the title.
#'
#' @return An object of class \code{seasexpl} containing: \itemize{
#' \item{\code{season}}{: matrix of (detrended) seasonal elements.}
#' \item{\code{season.exist}}{: \code{TRUE}/\code{FALSE} results of
#' seasonality test.} \item{\code{season.pval}}{: p-value of seasonality test
#' (Friedman test).} \item{\code{trend}}{: CMA estimate (using
#' \code{cmav}) or \code{NULL} if \code{trend=FALSE}.}
#' \item{\code{trend.exist}}{: \code{TRUE}/\code{FALSE} results of trend
#' test.} \item{\code{trend.pval}}{: p-value of trend test (Cox-Stuart).}
#' \item{\code{decomposition}}{: type of decomposition used.} }
#'
#' @author Nikolaos Kourentzes, \email{nikolaos@kourentzes.com}.
#'
#' @keywords ts
#'
#' @examples
#'
#' @keywords internal
#'
#' @example
#' \dontrun{
#' seasplot(referrals,outplot=1)
#' }
seasplot <- function(y,m=NULL,s=NULL,trend=NULL,colour=NULL,alpha=0.05,
outplot=c(1,0,2,3,4,5),decomposition=c("multiplicative","additive","auto"),
cma=NULL,labels=NULL,...)
{
# Defaults
decomposition <- match.arg(decomposition,c("multiplicative","additive","auto"))
outplot <- outplot[1]
# Get m (seasonality)
if (is.null(m)){
if (any(class(y) == "ts")){
m <- frequency(y)
} else {
stop("Seasonality not defined (y not ts object).")
}
}
# Get starting period in seasonality if available
if (is.null(s)){
if (any(class(y) == "ts")){
s <- start(y)
s <- s[2]
# Temporal aggregation can mess-up s, so override if needed
if (is.na(s)){s<-1}
} else {
s <- 1
}
}
# Make sure that labels argument is fine
if (!is.null(labels)){
if (length(labels) < m){
labels <- NULL
warning("Incorrect number of labels. It must be equal to m. Defaulting to labels=NULL.")
} else {
labels <- labels[1:m]
}
}
if (is.null(labels)){
labels <- paste(1:m)
}
n <- length(y)
if (min(y)<=0){
decomposition <- "additive"
}
# Override trend if cma is given
if (!is.null(cma)){
trend <- NULL
if (n != length(cma)){
stop("Length of series and cma do not match.")
}
}
# Calculate CMA
if ((is.null(trend) || trend == TRUE) && (is.null(cma))){
cma <- cmav(y=y,ma=m,fill=TRUE,outplot=FALSE)
}
# Test for changes in the CMA (trend)
if (is.null(trend)){
trend.pval <- coxstuart(cma)$p.value
trend.exist <- trend.pval <= alpha/2
trend <- trend.exist
} else {
trend.exist <- NULL
trend.pval <- NULL
}
if (trend == TRUE){
is.mult <- mseastest(y,m)$is.multiplicative
if (decomposition == "auto"){
if (is.mult == TRUE){
decomposition <- "multiplicative"
} else {
decomposition <- "additive"
}
}
if (decomposition == "multiplicative"){
ynt <- y/cma
} else {
ynt <- y - cma
}
title.trend <- "(Detrended)"
} else {
if (decomposition == "auto"){
decomposition <- "none"
}
ynt <- y
title.trend <- ""
cma <- NULL
}
ymin <- min(ynt)
ymax <- max(ynt)
ymid <- median(ynt)
# Fill with NA start and end of season
k <- m - (n %% m)
ks <- s-1
ke <- m - ((n+ks) %% m)
if (ke == m){ke <- 0}
if (ks == m){ks <- 0}
ynt <- c(rep(NA,times=ks),as.vector(ynt),rep(NA,times=ke))
ns <- length(ynt)/m
ynt <- matrix(ynt,nrow=ns,ncol=m,byrow=TRUE)
colnames(ynt) <- labels
rownames(ynt) <- paste("s",1:ns,sep="")
# Check seasonality with Friedman
if (m>1 && (length(y)/m)>=2){
# Check if ynt is a constant, if it is do not get it through friedman, as it fails
if (diff(range(colSums(ynt,na.rm=TRUE))) <= .Machine$double.eps ^ 0.5){
season.pval <- 1
} else {
season.pval <- friedman.test(ynt)$p.value
}
season.exist <- season.pval <= alpha
if (season.exist==TRUE){
title.season <- "Seasonal"
} else {
title.season <- "Nonseasonal"
}
} else {
season.pval <- NULL
season.exist <- NULL
title.season <- "Nonseasonal"
}
# Produce plots
if (outplot != 0){
yminmax <- c(ymin - 0.1*(ymax-ymin),ymax + 0.1*(ymax-ymin))
yminmax <- c(-1,1)*max(abs(ymid-yminmax))+ymid
if (is.null(season.pval)){
plottitle <- paste(title.trend, "\n", title.season,sep="")
} else {
plottitle <- paste(title.trend, "\n", title.season,
" (p-val: ",round(season.pval,3),")",sep="")
}
# Allow user to override plot defaults
args <- list(...)
if (!("main" %in% names(args))){
addtitle <- TRUE
} else {
addtitle <- FALSE
}
if (!("xlab" %in% names(args))){
args$xlab <- "Period"
}
if (!("ylab" %in% names(args))){
args$ylab <- ""
}
if (!("yaxs" %in% names(args))){
args$yaxs <- "i"
}
if (!("ylim" %in% names(args))){
args$ylim <- yminmax
}
# Remaining defaults
args$x <- args$y <- NA
}
if (outplot == 1){
# Conventional seasonal diagramme
if (is.null(colour)){
cmp <- colorRampPalette(RColorBrewer::brewer.pal(9,"YlGnBu")[4:8])(ns)
} else {
cmp <- rep(colour,times=ns)
}
# Additional plot options
if (!("xlim" %in% names(args))){
args$xlim <- c(1,m)
}
if (!("xaxs" %in% names(args))){
args$xaxs <- "i"
}
if (addtitle){
args$main <- paste0("Seasonal plot ",main=plottitle)
}
args$xaxt <- "n"
# Produce plot
do.call(plot,args)
for (i in 1:ns){
lines(ynt[i,],type="l",col=cmp[i])
}
lines(c(0,m+1),c(ymid,ymid),col="black",lty=2)
legend("topleft",c("Oldest","Newest"),col=c(cmp[1],cmp[ns]),lty=1,bty="n",lwd=2,cex=0.7)
axis(1,at=1:m,labels=labels)
}
if (outplot == 2){
# Seasonal boxplots
if (is.null(colour)){
cmp <- RColorBrewer::brewer.pal(3,"Set3")[1]
} else {
cmp <- colour
}
# Additional plot options
if (!("xlim" %in% names(args))){
args$xlim <- c(1,m)
}
if (addtitle){
args$main <- paste0("Seasonal boxplot ",main=plottitle)
}
args$x <- ynt
args$col <- cmp
# Produce plot
do.call(boxplot,args)
lines(c(0,m+1),c(ymid,ymid),col="black",lty=2)
}
if (outplot == 3){
# Subseries plots
if (is.null(colour)){
cmp <- RColorBrewer::brewer.pal(3,"Set1")
} else {
cmp <- colour
}
# Additional plot options
if (!("xlim" %in% names(args))){
args$xlim <- c(1,m*ns)
}
if (addtitle){
args$main <- paste0("Seasonal subseries ",main=plottitle)
}
if (!("xaxs" %in% names(args))){
args$xaxs <- "i"
}
args$xaxt <- "n"
# Produce plot
do.call(plot,args)
lines(c(1,ns),median(ynt[,1],na.rm=TRUE)*c(1,1),col=cmp[1],lwd=2)
for (i in 1:m){
lines((1+(i-1)*ns):(ns+(i-1)*ns),ynt[,i],type="o",col=cmp[2],pch=20,cex=0.75)
lines(c((1+(i-1)*ns),(ns+(i-1)*ns)),median(ynt[,i],na.rm=TRUE)*c(1,1),col=cmp[1],lwd=2)
if (i < m){
lines(c(1,1)*i*ns+0.5,yminmax,col="gray")
}
}
lines(c(0,m*ns+1),c(ymid,ymid),col="black",lty=2)
axis(1,at=seq(0.5+(ns/2),0.5+m*ns-ns/2,ns),labels=labels)
}
if (outplot == 4){
# Seasonal distribution
if (is.null(colour)){
cmp <- "royalblue"
}
qntl <- matrix(NA,nrow=9,ncol=m)
for (i in 1:m){
qntl[,i] <- quantile(ynt[!is.na(ynt[,i]),i], c(0, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 1))
}
# Additional plot options
if (!("xlim" %in% names(args))){
args$xlim <- c(1,m)
}
if (addtitle){
args$main <- paste0("Seasonal distribution ",main=plottitle)
}
if (!("xaxs" %in% names(args))){
args$xaxs <- "i"
}
args$xaxt <- "n"
# Produce plot
do.call(plot,args)
polygon(c(1:m,rev(1:m)),c(qntl[7,],rev(qntl[1,])),col=gray(0.8),border=NA)
polygon(c(1:m,rev(1:m)),c(qntl[6,],rev(qntl[2,])),col="lightblue",border=NA)
polygon(c(1:m,rev(1:m)),c(qntl[5,],rev(qntl[3,])),col="skyblue",border=NA)
lines(1:m,qntl[4,],col=cmp,lwd=2)
lines(c(0,m*ns+1),c(ymid,ymid),col="black",lty=2)
legend("topleft",c("Median","25%-75%","10%-90%","MinMax"),col=c(cmp,"skyblue","lightblue",gray(0.8)),lty=1,bty="n",lwd=2,cex=0.7)
axis(1,at=1:m,labels=labels)
box()
}
if (outplot == 5){
dnst <- matrix(NA,nrow=m,ncol=512)
for (i in 1:m){
tmp <- density(ynt[!is.na(ynt[,i]),i], bw = "SJ", n = 512, from = yminmax[1], to = yminmax[2])
dnst[i,] <- tmp$y
if (i == 1){
llc <- tmp$x
}
}
cmp <- c(rgb(0,0,0,0), colorRampPalette((RColorBrewer::brewer.pal(9,"Blues")[3:9]))(100))
# Additional plot options
if (addtitle){
args$main <- paste0("Seasonal density ",main=plottitle)
}
args$xaxt <- "n"
args$col <- cmp
args$x <- 1:m
args$y <- llc
args$z <- dnst
# Produce plot
do.call(image,args)
lines(colMeans(ynt,na.rm=TRUE),type="o",lty=1,bg=RColorBrewer::brewer.pal(3,"Set1")[1],pch=21,cex=1.1,lwd=2)
lines(c(0,m*ns+1),c(ymid,ymid),col="black",lty=2)
box()
axis(1,at=1:m,labels=labels)
}
out <- list(season=ynt,season.exist=season.exist,season.pval=season.pval,
trend=cma,trend.exist=trend.exist,trend.pval=trend.pval,decomposition=decomposition)
out <- structure(out,class="seasexpl")
return(out)
}
#' @export
#' @method summary seasexpl
summary.seasexpl <- function(object,...){
print(object)
}
#' @export
#' @method print seasexpl
print.seasexpl <- function(x,...){
cat("Results of statistical testing\n")
if (!is.null(x$trend.exist)){
cat(paste0("Evidence of trend: ",x$trend.exist, " (pval: ",round(x$trend.pval,3),")\n"))
} else {
cat("Presence of trend not tested.\n")
}
cat(paste0("Evidence of seasonality: ",x$season.exist, " (pval: ",round(x$season.pval,3),")\n"))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.