Como decompor uma série temporal com vários componentes sazonais?

Tenho uma série temporal que contém componentes sazonais duplos e gostaria de decompor a série nos seguintes componentes da série temporal (tendência, componente sazonal 1, componente sazonal 2 e componente irregular). Até onde eu sei, o procedimento STL para decompor uma série em R permite apenas um componente sazonal, então tentei decompor a série duas vezes. Primeiro, definindo a frequência como o primeiro componente sazonal usando o seguinte código:

ser = ts(data, freq=48)
dec_1 = stl(ser, s.window="per")

Em seguida, decompus o componente irregular da série decomposta ( dec_1) definindo a frequência como o segundo componente sazonal, de modo que:

ser2 = ts(dec_1$time.series[,3], freq=336)
dec_2 = stl(ser2, s.window="per")

Não estou muito confiante com essa abordagem. E eu gostaria de saber se existem outras maneiras de decompor uma série com várias sazonalidades. Além disso, notei que a tbats()função no pacote de previsão R permite ajustar um modelo a uma série com várias sazonalidades, no entanto, não diz como decompor uma série com ele.

forecastbats()tbats()$x(t)$$t$

Veja http://robjhyndman.com/papers/complex-seasonality/ para obter as fórmulas e Hyndman et al (2008) para uma melhor descrição dos modelos de ETS. BATS e TBATS são uma extensão do ETS.

Por exemplo:

fit <- bats(myTimeseries)
fit$x

Nesse caso, cada linha de xestará em harmônica do tipo fourier.

Há também plot.tbats()e plot.bats()funções para automaticamente se decompor e visualizar os componentes.