Can anyone explain what is happening in the stl function of R? I am recently working with seasonal-trend decomposition.
Yet I am not that familiar with the approach that R is using.
Can anyone one kindly explain the mechanism of the stl function?
OP questions copied from comments:


*

*In step 2, is $C^{k+1}_v$ again a complete time series (joining all the cycle subseries accordingly)?

*In step 3, I dont really understand how the filter looks like. What does "the filter consists of a moving average ... followed by another moving average ... followed by another moving average, followed by a loess smoothing" means? Does it means smoothing the series 4 times, each done by the corresponding technique?

 A: This answer is in response to the specific questions around steps 2 and 3.
Regarding step 2, yes, $C^{k+1}_v$ is again a complete time series. In fact it extends the original time series (i.e. the one being decomposed) by pre- and post-padding, since the low-pass filter in step 3 needs that padding to avoid truncating the original series.
Regarding step 3, yes, the filter smooths the series four times during any given run of the inner loop.
If you're feeling adventurous...
You can look at the source code itself. From the R shell, just type
$ stl

and you will see that R's stl function is really just a wrapper around a Fortran subroutine:
z <- .Fortran(C_stl, x, n, as.integer(period), as.integer(s.window),
    as.integer(t.window), as.integer(l.window), s.degree,
    t.degree, l.degree, nsjump = as.integer(s.jump), ntjump = as.integer(t.jump),
    nljump = as.integer(l.jump), ni = as.integer(inner),
    no = as.integer(outer), weights = double(n), seasonal = double(n),
    trend = double(n), double((n + 2 * period) * 5))

The STL docs say that the Fortran function lives at "netlib":


*

*stl: http://www.netlib.org/a/stl

*loess: http://www.netlib.org/a/loess
In stl you'll find the fts subroutine, which has the moving average calls:
subroutine fts(x,n,np,trend,work)

integer n, np
real x(n), trend(n), work(n)

call ma(x,n,np,trend)
call ma(trend,n-np+1,np,work)
call ma(work,n-2*np+2,3,trend)
return
end

That way you can get as down-in-the-weeds as you'd like.
