Dual seasonal cycles in ts object

Question

I want to strip out seasonality from a ts. This particular ts is daily, and has both yearly and weekly seasonal cycles (frequency 365 and 7).

In order to remove both, I have tried conducting stl() on the ts with frequency set to 365, before extracting trend and remainders, and setting the frequency of the new ts to 7, and repeat.

This doesn't seem to be working very well and I am wondering whether it's my approach, or something inherent to the ts which is causing me problems. Can anyone critique my methodology, and perhaps recommend an alternate approach?

score 7 · Answer 1 · answered Jan 08 '14 at 21:58

7

There is a very easy way to do it using a TBATS model implemented in the forecast package. Here is an example assuming your data are stored as x:

library(forecast)
x2 <- msts(x, seasonal.periods=c(7,365))
fit <- tbats(x2)
x.sa <- seasadj(fit)

Details of the model are described in De Livera, Hyndman and Snyder (JASA, 2011).

answered Jan 08 '14 at 21:58

Rob Hyndman

30,301
7
73
85

Thanks so much for your response, Rob. Your forecast package is the best, and it looks as though I've only scratched the surface. The seasadj(fit) line in your answer returns an error which I am struggling to interpret, perhaps you could help? "Error in `[.default`(comp, , "season") : subscript out of bounds" – Jonathan Mulligan Jan 09 '14 at 13:00
Can you send me a minimal example to replicate the error? This functionality is relatively new and not well-tested. It worked on the example I tried. Send bug reports to https://github.com/robjhyndman/forecast/issues?state=open – Rob Hyndman Jan 10 '14 at 00:30

score 2 · Answer 2 · answered Jan 08 '14 at 11:40

An approach that can handle not only seasonal components (cyclically reoccurring events) but also trends (slow shifts in the norm) admirably is stl(), specifically as implemented by Rob J Hyndman.

The decomp function Hyndman gives there (reproduced below) is very helpful for checking for seasonality and then decomposing a time series into seasonal (if one exists), trend, and residual components.

decomp <- function(x,transform=TRUE)
{
  #decomposes time series into seasonal and trend components
  #from http://robjhyndman.com/researchtips/tscharacteristics/
  require(forecast)
  # Transform series
  if(transform & min(x,na.rm=TRUE) >= 0)
  {
    lambda <- BoxCox.lambda(na.contiguous(x))
    x <- BoxCox(x,lambda)
  }
  else
  {
    lambda <- NULL
    transform <- FALSE
  }
  # Seasonal data
  if(frequency(x)>1)
  {
    x.stl <- stl(x,s.window="periodic",na.action=na.contiguous)
    trend <- x.stl$time.series[,2]
    season <- x.stl$time.series[,1]
    remainder <- x - trend - season
  }
  else #Nonseasonal data
  {
    require(mgcv)
    tt <- 1:length(x)
    trend <- rep(NA,length(x))
    trend[!is.na(x)] <- fitted(gam(x ~ s(tt)))
    season <- NULL
    remainder <- x - trend
  }
  return(list(x=x,trend=trend,season=season,remainder=remainder,
    transform=transform,lambda=lambda))
}

As you can see it uses stl() (which uses loess) if there is seasonality and penalized regression splines if there is no seasonality.

Prasanna Nandakumar · Answer 3 · 2014-01-08T11:37:02.077

Check if this is useful:
Start and End Values depends on your Data - Change the Frequency values accordingly

splot <- ts(Data1, start=c(2010, 2), end=c(2013, 9), frequency=12)

additive trend, seasonal, and irregular components can be decomposed using the stl() Function

fit <- stl(splot, s.window="period")
monthplot(splot) 
library(forecast)
vi <-seasonplot(splot)

vi should give seperate values for a seasonal indices

Also Check the below one:

splot.stl <- stl(splot,s.window="periodic",na.action=na.contiguous)
    trend <- splot.stl$time.series[,2]
    season <- splot.stl$time.series[,1]
    remainder <- splot - trend - season

Dual seasonal cycles in ts object

3 Answers3