How to decompose a time series with multiple seasonal components? I have a time series that contains double seasonal components and I would like to decompose the series into the following time series components (trend, seasonal component 1, seasonal component 2 and irregular component). As far as I know, the STL procedure for decomposing a series in R only allows one seasonal component, so I have tried decomposing the series twice. First, by setting the frequency to be the first seasonal component using the following code:
ser = ts(data, freq=48)
dec_1 = stl(ser, s.window="per")

Then, I decomposed the irregular component of the decomposed series (dec_1) by setting the frequency to be the second seasonal component, such that:
ser2 = ts(dec_1$time.series[,3], freq=336)
dec_2 = stl(ser2, s.window="per")

I'm not very confident with this approach. And I would like to know if there are any other ways to decompose a series that has multiple seasonalities. Also,I have noticed that the tbats() function in the R forecast package allows one to fit a model to a series with multiple seasonalities however, it doesn't say how to decompose a series with it.
 A: R's forecast package now has a function mstl() to handle multiple seasonal time series decomposition.
This page has got more details how to use it: https://pkg.robjhyndman.com/forecast/reference/mstl.html
A: R's forecast package bats() and tbats() functions can fit BATS and TBATS models to the data. The functions return lists with a class attribute either "bats" or "tbats". One of the elements on this list is a time series of state vectors, $x(t)$, for each time, $t$. 
See http://robjhyndman.com/papers/complex-seasonality/ for the formula's and Hyndman et al (2008) for a better description of ETS models. BATS and TBATS are an extension of ETS.
For example:
fit <- bats(myTimeseries)
fit$x

In this case, each row of x will be on fourier-like harmonic.
There are also plot.tbats() and plot.bats() functions to automatically decompose and view the components.
A: The facebook prophet package supports multiple seasonalities.
Yearly, weekly and daily seasonalities are built-in but custom seasonalities can be specified.
Here is a custom monthly seasonality:
df <- ...     # data to build model on or decompose
future <- ... # data to make forecasts on

m <- prophet(weekly.seasonality = FALSE)
m <- add_seasonality(m, name = 'monthly', period = 30.5, fourier.order = 5)
m <- fit.prophet(m, df)
forecast <- predict(m, future)
prophet_plot_components(m, forecast)

If predict() is called without passing in a data frame, then it will decompose the time series used to build the model.
