STLF function in the FORECAST package

Question

I am trying to forecast a yearly time series on a weekly bases (52 weeks a year and I have 164 weeks data). As the frequency is larger than 24, R advices me to use "stlf" rather than "ets" to avoid seasonality being ignored. The "stlf" function works perfectly well and I got the following:

> WR.ets<-stlf(WeeklyReferral,method="ets")
> summary(WR.ets)

Forecast method: STL +  ETS(A,A,N)

Model Information:
ETS(A,A,N) 

Call:
 ets(y = x.sa, model = etsmodel) 

  Smoothing parameters:
    alpha = 0.0262 
    beta  = 1e-04 

  Initial states:
    l = 93.1548 
    b = 0.1159 

  sigma:  12.6201

     AIC     AICc      BIC 
1675.954 1676.205 1688.353 

Error measures:
                     ME     RMSE      MAE       MPE     MAPE      MASE
Training set -0.1869514 12.62011 9.790321 -2.589141 11.12905 0.5990874

Forecasts:
         Point Forecast     Lo 80     Hi 80     Lo 95     Hi 95
2013.423       95.39869  79.22537 111.57201  70.66373 120.13364
2013.442       95.03434  78.85538 111.21330  70.29075 119.77793
...............................................................

The point forecast gives the mean of the predicted value. However, what I want is the actual forecast value rather than the mean. Thus I am trying to understand how it works and break down the steps. I use "stl" decomposition firstly on the time series

temp<-stl(WeeklyReferral,s.window="periodic", robust=TRUE)
> temp
 Call:
 stl(x = WeeklyReferral, s.window = "periodic", robust = TRUE)

Components
Time Series:
Start = c(2010, 15) 
End = c(2013, 22) 
Frequency = 52 
            seasonal     trend   remainder
2010.269   7.1597729  82.33453  -0.4943046
2010.288  -1.4283001  82.69446   5.7338358
..........................................
2013.404   8.0046803 117.74388  -0.7485615

Then I use "trend+remainder" as the new time series to forecast for 3 months (12 periods). I use the last state vector obtained by "stlf" function as the initial state vector in my following formulas. And add the seasonal values at the same week last year back to the forecasted values as the "stlf" function shows the model is ETS(A,A,N).

y<-c(rep(NA,13))
l<-c(rep(NA,13))
b<-c(rep(NA,13))
e<-c(rep(NA,12))
alpha<-0.0262
beta<-0.0001

y[1]<-117.74388-0.7485615
l[1]<-109.66913
b[1]<-0.11284923

for (j in 1:1000){
  for(i in 2:13){
e[i-1]=rnorm(sd=12.6201,n=1)
b[i]<-b[i-1]+beta*e[i-1]
l[i]<-l[i-1]+b[i-1]+alpha*e[i-1]
y[i]<-l[i-1]+b[i-1]+e[i-1]+temp$time.series[i+164-52,1]
}}

Am I right?

I tried to use "ets" function on the new decomposed time series and it gave different parameters (alpha, beta, l,b, sigma) and it didn't give any forecasted values.

Any opinions are appreciated.

What do you mean "I want the actual forecast value rather than the mean"? A point forecast, by definition, is a summary of the probability distribution of a future random variable, and the mean is the most common way of summarizing it. — Rob Hyndman, Jul 28 '14 at 11:37
Because I need to use the forecast values (y) in the calculation of other values (say I want to calculate x). So what I want is to find the distribution of x at each forecasted time period by calculating it 1000 times. If I use the point forecast then I will only have the mean value of x. Thus I want to have y, which I am forecasting in the above model, varing in that 1000 times. — MichelleX, Jul 28 '14 at 12:15
Do you mean you want the whole distribution rather than just point forecasts? Your use of `x` and `y` makes little sense to me. — Rob Hyndman, Jul 28 '14 at 12:23
@Rob Hyndman. I am forecasting the referral number in a hospital so I can calculate the number of patients in the queue. There are cancellations etc., which have their own distribution. So I want to calculate the number of patients in the queue 1000 times to find the possible distribution at each future time point. Thus I want the forecast to give a predicted referral number each time I calculate rather than use the same mean value. — MichelleX, Jul 28 '14 at 12:33
So I am not sure if I need the whole distribution, but need the future random variable not just the mean. — MichelleX, Jul 28 '14 at 12:54

score 4 · Answer 1 · answered Jul 28 '14 at 22:43

As far as I can tell from the comments above, you actually want to simulate future sample paths from the model rather than obtain point forecasts or interval forecasts. The following code will do it.

# STL decomposition
temp <- stl(WeeklyReferral, s.window="periodic", robust=TRUE)
# Seasonally adjusted data
sa <- seasadj(temp)
seascomp <- tail(temp$time.series,52)[,1]
# ETS model
fit <- ets(sa, "ZZN")
# Simulations from ETS model with re-seasonalization
sim <- matrix(0, nrow=52, ncol=1000)
for(i in 1:1000)
  sim[,i] <- simulate(fit, nsim=52) + seascomp

The matrix sim contains 1000 future sample paths each of length 52.

One more question, Professor. When I apply "stlf" function, it shows that "Forecast method: STL + ETS(A,N,N)". However, when I apply STL decomposition manually and apply ETS, the fitted model I got from ETS is "ETS(A,A,N)". Any possible reason for this? And also how do I know whether the seasonal components need to be added or multiplied to the forecast value. In other words, when I apply the STL decomposition, how to test if the model is additive or multiplicative? — MichelleX, Aug 05 '14 at 12:08

STLF function in the FORECAST package

1 Answers1