This is a long post so I hope you can bear with me, and please correct me where I'm wrong.
My goal is to produce a daily forecast based on 3 or 4 weeks of historical data.
The data is 15 minute data of the local load of one of a transformer lines. Im having trouble finding the model order of a seasonal ARIMA process. Consider the electricity demand time series:
Original Time Series http://i.share.pho.to/80d86574_l.png
When the first 3 weeks are taken as subset and differenced the follwing ACF/PACF plots are computed:
Subset http://i.share.pho.to/5c165aef_l.png
First difference http://i.share.pho.to/b7300cc2_l.png
Seasonal and first difference http://i.share.pho.to/570c5397_l.png
This looks like the series is kinda stationary. But the seasonality could also be weekly (see Seasonal difference week and second order differences [here]http://share.pho.to/3owoq, what do you think?)
So let's conclude that the model take shape of :
$$
ARIMA(p,1,q)(P,1,Q)_{96}
$$
In the last figure a distinct spike at lag 96 indicates a seasonal MA(1) component (maybe AR(1) could be aswell as there is a distinct spike in PACF aswell). The spikes at lag 1:4 indicate a MA(4) component which corresponds with exponential decay in the PACF with a bit of imagination. Thus the initial model manually selected could be:
$$
ARIMA(0,1,4)(0,1,1)_{96}
$$
with
Series: x
ARIMA(0,1,4)(0,1,1)[96]
Coefficients:
ma1 ma2 ma3 ma4 sma1
-0.2187 -0.2233 -0.0996 -0.0983 -0.9796
s.e. 0.0231 0.0234 0.0257 0.0251 0.0804
sigma^2 estimated as 364612: log likelihood=-15138.91
**AIC=30289.82 AICc=30289.87 BIC=30323.18**
The auto.arima function computes the following model (with stepwise and approximation on TRUE, else it takes to long to converge):
$$
ARIMA(1,1,1)(2,0,2)_{96}
$$
with
Series: x
ARIMA(1,1,1)(2,0,2)[96]
Coefficients:
ar1 ma1 sar1 sar2 sma1 sma2
0.7607 -1.0010 0.4834 0.4979 -0.3369 -0.4168
s.e. 0.0163 0.0001 0.0033 0.0116 0.0216 0.0255
sigma^2 estimated as 406766: log likelihood=-15872.02
**AIC=31744.99 AICc=31745.05 BIC=31784.25**
Which means no seasonal differencing is applied. Here
are the residuals of the both models. The Ljung Box statistic give a very small p value, indicating that there is still some autocorrelation present (? correct me if im wrong).
Forecasting
Thus to determine which is better, an out-sample accuracy test is then the best. So for both models an forecast is made 24 hours ahead which is compared with each other. The results are:
auto.arima http://i.share.pho.to/5d1dd934_l.png
manual model http://i.share.pho.to/7ca69c97_l.png
Auto:
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set -2.586653 606.3188 439.1367 -1.284165 7.599403 0.4914563 -0.01219792 NA
Test set -330.144797 896.6998 754.0080 -7.749675 13.268985 0.8438420 0.70219229 1.617834
Manual
ME RMSE MAE MPE MAPE MASE ACF1 Theil's U
Training set 2.456596e-03 589.1267 435.6571 -0.7815229 7.509774 0.4875621 -0.002034122 NA
Test set 2.878919e+02 919.7398 696.0593 3.4756363 10.317420 0.7789892 0.731013599 1.281764
Questions
As you can think of this is an analysis on the first three weeks of a dataset. I'm struggling in my mind with the following questions:
- How do I select the best ARIMA model (by trying all different orders and checking the best MASE/MAPE/MSE? where the selection of performance measurement can be a discussion in it's own..)
- If I generate a new model and forecast for every new day forecast (as in online forecasting), do I need to take the yearly trend into account and how? (as in such a small subset my guess would be that the trend is neglible)
- Would you expect that the model order stays the same throughout the dataset, i.e. when taking another subset will that give me the same model?
- What is a good way, within this method to cope with holidays? Or is ARIMAX with external holiday dummies needed for this?
- Do I need to use Fourier series approach to try models with
seasonality=672as discussed in Long seasonal periods
? - If so would this be like
fit<-Arima(timeseries,order=c(0,1,4), xreg=fourier(1:n,4,672)(where the function fourier is as defined in Hyndman's blog post) - Are initial P and Q components included with the fourier series?
Most theoretical knowledge obtained from FPP, great stuff!
Before advising on using exponential smoothing or (dynamic)linear regression this is also being worked on to compare.
Data
https://www.dropbox.com/sh/mzx61sskya5ze6x/Zq3A7Q6htH/trafo.txt
Code
data<-read.csv("file", sep=";")
load<-data[,3]
I removed the few zero values with week before values
stepback<-672
load[is.na(load)] <- 0 # Assumed no 0's in first 672 values!
idx <- which(load == 0)
idx <- idx[which(idx>stepback)]
load[idx] <- load[idx-stepback]
ED<-ts(load,start=0, end=c(760,96),frequency=96)
x<-window(ED,start=0, end=c(20,96))
It is also possible to post a reproducible example but this will make the post even longer, but possible if needed. So if there is anything I should provide please let me know.
Best Answer
Out of sample risk estimates are the gold standard for performance evaluation, and therefore for model selection. Ideally, you cross-validate so that your risk estimates are averaged over more data. FPP explains one cross-validation method for time series. See Tashman for a review of other methods:
Tashman, L. J. (2000). Out-of-sample tests of forecasting accuracy: an analysis and review. International Journal of Forecasting, 16(4), 437–450. doi:10.1016/S0169-2070(00)00065-0
Of course, cross-validation is time consuming and so people often resort to using in-sample criteria to select a model, such as AIC, which is how auto.arima selects the best model. This approach is perfectly valid, if perhaps not as optimal.
I'm not sure what you mean by yearly trend. Assuming you mean yearly seasonality, there's not really any way to take it into account with less than a year's worth of data.
I would expect that barring some change to how the data are generated, the most correct underlying model will be the same throughout the dataset. However, that's not the same as saying that the model selected by any procedure (such as the procedure used by auto.arima) will be the same if that procedure is applied to different subsets of the data. This is because the variability due to sampling will result in variability in the results of the model selection procedure.
External holiday dummies is the best approach.
You need to do something, because as mentioned in that article, the arima function in R does not support seasonal periods greater than 350. I've had reasonable success with the Fourier approach. Other options include forecasting after seasonal decomposition (also covered in FPP), and exponential smoothing models such as bats and tbats.
That looks correct. You should experiment with different numbers of terms. Note that there is now a
fourierfunction in the forecast package with a slightly different specification that I assume supersedes the one on Hyndman's blog. See the help file for syntax.I'm not sure what you're asking here. P and Q usually refer to the degrees of the AR and MA seasonal components. Using the fourier approach, there are no seasonal components and instead there are covariates for fourier terms related to season. It's no longer seasonal ARIMA, it's ARIMAX where the covariates approximate the season.