STLF function in the FORECAST package
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.
Answer
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.