R devel - Apr 2011 - Proposed modification to decompose() and plot.decomposed.ts()

The decompose() function truncates the seasonal component
unnecessarily. I've modified the function to fix this problem, and
also added the original data to the object returned (to enable better
plotting).

I've also modified the plot.decomposed.ts() function so that it plots
the original data in the top panel rather than the reconstructed data.
The difference between the two is that the reconstructed data has
missing values at the two ends whereas the original data does not.

Please find the revised code below. I hope this can find its way into
the stats package in a future release.

Regards,
Rob Hyndman


# Replacement for decompose function in the stats package
# Only changes are (1) the seasonal component is complete instead of truncated
# and (2) the original data is included in the returned object.

decompose <- function (x, type = c("additive",
"multiplicative"),
filter = NULL)
{
? ?type <- match.arg(type)
? ?l <- length(x)
? ?f <- frequency(x)
? ?if (f <= 1 || length(na.omit(x)) < 2 * f)
? ? ? ?stop("time series has no or less than 2 periods")
? ?if (is.null(filter))
? ? ? ?filter <- if (!f%%2)
? ? ? ? ? ?c(0.5, rep(1, f - 1), 0.5)/f
? ? ? ?else rep(1, f)/f
? ?trend <- filter(x, filter)
? ?season <- if (type == "additive")
? ? ? ?x - trend
? ?else x/trend
? ?season <- na.omit(c(as.numeric(window(season, start(x) +
? ? ? ?c(1, 0), end(x))), as.numeric(window(season, start(x),
? ? ? ?start(x) + c(0, f)))))
? ?periods <- l%/%f
? ?index <- c(0, cumsum(rep(f, periods - 2)))
? ?figure <- numeric(f)
? ?for (i in 1L:f) figure[i] <- mean(season[index + i])
? ?figure <- if (type == "additive")
? ? ? ?figure - mean(figure)
? ?else figure/mean(figure)
? ?# Following lines are only changes
? ?seasonal <- ts(rep(figure, periods+1)[1:l], start = start(x), frequency =
f)
? ?structure(list(x=x, seasonal = seasonal, trend = trend, random if (type =? ?
? ?"additive") x - seasonal - trend else x/seasonal/trend,
? ? ? ?figure = figure, type = type), class = "decomposed.ts")
}

# Changed definition of observed to use passed data rather than construct on fly

plot.decomposed.ts <- function (x, ...)
{
? ?plot(cbind(observed = x$x, trend = x$trend, seasonal = x$seasonal,
random = x$random),
? ? ? ?main = paste("Decomposition of", x$type, "time
series"), ...)
}



-- 
-------------------------------------------------------------
Rob J Hyndman
Professor of Statistics, Monash University
www.robjhyndman.com