Is there a remedy for removing autocorrelations from residuals of seasonally fitted ARIMA model?

I fitted a number of SARIMA models using R and chose the ARIMA(0,0,0)(3,1,0)[12] as the best fitted model to the univariate data with 180 points (periodicity=12). This model is chosen as the best model according to the criteria of lowest MAPE among other fitted 624 models.

The residuals of the model violates the assumption of independently distributed residuals (and same for the 2nd best, 3rd best model etc.). Actually the residuals are also non-normally distributed; however the model is fitted with the method of conditional sum of squares in order to bypass the violation of normality assumption.

In the data, the most of the values are close to zero and this does not allow any data transformation.

The data represent the evolution of coefficents of a 11th degree polynomial equation (in total 15 equations representing different years of electricity load duration curves). The purpose is to forecast the coefficients of e.g. the 16th equation and so the corresponding load duration curve.

Can anybody sugggest/provide any solutions to this case?

x=c(1.887090e+04, -6.023007e+00,  1.193635e-02, -1.455856e-05,  1.064251e-08, -4.953592e-12,  1.517229e-15, -3.090332e-19,
4.137144e-23, -3.491891e-27,  1.682794e-31, -3.527046e-36,  1.904962e+04, -7.394189e+00,  1.600849e-02, -2.077511e-05,
1.585519e-08,-7.587987e-12,    2.363570e-15, -4.859251e-19,  6.534816e-23, -5.525202e-27,  2.663420e-31, -5.580438e-36,
2.009098e+04, -1.061082e+01,  2.319182e-02, -2.917768e-05,  2.171827e-08, -1.019917e-11,  3.133564e-15, -6.379905e-19,
8.520995e-23, -7.168462e-27,  3.442102e-31, -7.188143e-36,  2.067028e+04, -8.034999e+00,  1.761326e-02, -2.240562e-05,
1.680919e-08, -7.961614e-12,  2.469832e-15, -5.081494e-19,  6.861040e-23, -5.835236e-27,  2.831898e-31, -5.974519e-36,
2.233604e+04, -1.033148e+01,  2.287039e-02, -2.952031e-05,  2.255568e-08, -1.086351e-11,  3.419260e-15, -7.123005e-19,
9.720229e-23, -8.341734e-27,  4.079166e-31, -8.660882e-36,  2.392045e+04, -8.246481e+00,  1.585412e-02, -2.056180e-05,
1.636424e-08, -8.253437e-12,  2.710813e-15, -5.858824e-19,  8.245204e-23, -7.258003e-27,  3.624039e-31, -7.827743e-36,
2.636514e+04, -9.886355e+00,  1.951992e-02, -2.504930e-05,  1.963158e-08, -9.789139e-12,  3.190186e-15, -6.856046e-19,
9.606813e-23, -8.427664e-27,  4.196799e-31, -9.046539e-36,  2.866210e+04, -8.866902e+00,  1.734494e-02, -2.387617e-05,
1.957175e-08, -9.993900e-12,  3.300201e-15, -7.152619e-19,  1.008517e-22, -8.892694e-27,  4.448060e-31, -9.626143e-36,
3.002254e+04, -1.007403e+01,  2.151203e-02, -2.984675e-05,  2.427803e-08, -1.226036e-11,  3.997630e-15, -8.550747e-19,
1.190499e-22, -1.037815e-26,  5.140218e-31, -1.103334e-35,  2.929311e+04, -1.123255e+01,  2.282206e-02, -2.968240e-05,
2.323868e-08, -1.146069e-11,  3.677709e-15, -7.777557e-19,  1.073806e-22, -9.301478e-27,  4.584147e-31, -9.800725e-36,
3.306894e+04, -1.396117e+01,  2.326777e-02, -2.724425e-05,  2.023428e-08, -9.690231e-12,  3.055811e-15, -6.392630e-19,
8.763020e-23, -7.552202e-27,  3.707622e-31, -7.901994e-36,  3.491666e+04, -1.315883e+01,  2.554492e-02, -3.194439e-05,
2.437661e-08, -1.184053e-11,  3.762542e-15, -7.896499e-19,  1.082565e-22, -9.310722e-27,  4.554895e-31, -9.664092e-36,
3.775600e+04, -2.101521e+01,  4.695457e-02, -6.000206e-05,  4.510264e-08, -2.134088e-11,  6.600784e-15, -1.352465e-18,
1.817468e-22, -1.538166e-26,  7.429410e-31, -1.560507e-35,  3.699341e+04, -1.019327e+01,  1.761360e-02, -2.428662e-05,
2.084200e-08, -1.112473e-11,  3.796505e-15, -8.415154e-19,  1.204392e-22, -1.072641e-26,  5.402195e-31, -1.174885e-35,
4.009280e+04, -1.887174e+01,  3.441926e-02, -4.161190e-05,  3.152055e-08, -1.535050e-11,  4.911316e-15, -1.040003e-18,
1.440215e-22, -1.251900e-26,  6.190925e-31, -1.327693e-35)

fit=arima(x, order = c(0, 0, 0),seasonal = list(order = c(3, 1, 0), period =12),method=c("CSS"))

par(mfrow=c(1,2))
x1<-acf(fit$residuals,180,ylab="Sample ACF",main ="",xaxt="n") axis(1, at=seq(0, 15, by=2), labels = TRUE) abline(v=(seq(0,15,1)), col="black", lty="dotted") x2<-pacf(fit$residuals,180,ylab="Sample PACF",main ="",xaxt="n")
axis(1, at=seq(0, 15, by=2), labels = TRUE)
abline(v=(seq(0,15,by=1)), col="black", lty="dotted")

• Sometimes a bit of autocorrelation can be tolerated if the alternative is an overfitted model. Model selection by AIC, BIC or MAPE would not necessarily give the model with the least autocorrelation, but there is a reason for that. Given the limited amount of data, trying to explain all observable patterns (such as autocorrelation) may lead to overfitting. – Richard Hardy May 10 '15 at 20:15
• @RichardHardy, thx for your reply. The sample acf and pacf indicate significant spikes at least for three consecutive seasonal lags (suggesting not a good fit ). If an ARIMA(0,0,0)X(6,1,6) [12] model is fittted to the data by "ML Method", then there are not any significant spikes in sample acf & pacf of the residuals. This is an overfitted model; however the sample acf & pacf of the only once differenced data suggests this model. Although it is overfitted, this seems to be the possible model for this data. I think in this case, overfitting should not be a big problem. what is your opinion? – Dirk May 11 '15 at 12:30
• It depends on what you want to use the model for. Asymptotically, the model selected by AIC should give the best forecasts in terms of mean squared error. Meanwhile, BIC is a consistent model selector in a certain sense; you could use it if your interest was explanation rather than prediction. You could also try time series cross validation as an alternative to AIC or BIC. In general, I would be cautious about choosing a model with a higher AIC value over another model with a lower AIC value, even if the former model produces nicer residuals. – Richard Hardy May 11 '15 at 12:39
• @RichardHardy, I compared the model ARIMA(0,0,0)X(3,1,0) [12] with ARIMA(0,0,0)X(6,1,6) [12] and saw that the first model has a log likelihood = -1261.85, aic = 2531.7 whereas the second has log likelihood = -1206.36, aic = 2438.73. The second model with higher orders seems as a better fit for the data; although the model is overparameterized. I think the data make this case as an exceptional case. – Dirk May 11 '15 at 13:51
• AIC of the second model is lower than that of the first model. Lower is better when it comes to AIC. Hence, AIC points to using the second model. Note also that the likelihood and the AIC are not comparable across models with different dependent variables. In SARIMA models, this may happen, for example, if you discard $p_1$ observations when calculating the likelihood of an AR($p_1$) model but $p_2$ observations when calculating the likelihood of an AR($p_2$) model. The dependent variable differs by $| p_1-p_2 |$ observations. You should carefully adjust for that. – Richard Hardy May 11 '15 at 14:21

This does not answer the general question "Is there a remedy for removing autocorrelations from residuals of seasonally fitted ARIMA model?" but relates to the specific question of this data.

I used the auto.arima() function in the "forecast" package in R, and the residuals from ARIMA(0,0,0)(2,1,0)[12] seem to come out clean. Here's the code:

plot(x)

xts <- ts(x,start=1,frequency=12) #convert to a time series

plot(xts)

library(fpp)  #load forecasting package (with files related to Hyndman's textbook)

mod1 <- auto.arima(xts)

mod1

acf(mod1$residual) pacf(mod1$residual)


• Also, I'm not sure MAPE is the best criterion for this data, since I would think you would be more concerned with the spike than the percent error when the data are very low. – zbicyclist May 12 '15 at 3:44
• thx for your reply. The plotted correlogram of yours show only one seasonal lag but not all the seasonal lags. If you increase the lags, u will see that the residuals are correlated at most of the seasonal lags. – Dirk May 14 '15 at 12:15
• par(mfrow=c(1,2)) x1<-acf(mod1$residual,180,ylab="Sample ACF",main ="",xaxt="n") axis(1, at=seq(0, 15, by=2), labels = TRUE) abline(v=(seq(0,15,1)), col="black", lty="dotted") x2<-pacf(mod1$residual,180,ylab="Sample PACF",main ="",xaxt="n") axis(1, at=seq(0, 15, by=2), labels = TRUE) abline(v=(seq(0,15,by=1)), col="black", lty="dotted") – Dirk May 14 '15 at 12:20

The degree of the SAR and SMA terms should be set at least to 6th order in order to get residuals with no correlations. This is due to the available data. And for model selection AIC, AICc and BIC cannot be used since the residuals are not normally distributed. That's why only MAPE should be considered, in my opinion.

fit=arima(x, order = c(0, 0, 0),seasonal = list(order = c(6, 1, 6), period =12),method=c("ML"))
par(mfrow=c(1,2))
x1<-acf(fit$residuals,180,ylab="Sample ACF",main ="",xaxt="n") axis(1, at=seq(0, 180, by=12), labels = TRUE) abline(v=(seq(0,180,12)), col="black", lty="dotted") x2<-pacf(fit$residuals,180,ylab="Sample PACF",main ="",xaxt="n")
axis(1, at=seq(0, 180, by=12), labels = TRUE)
abline(v=(seq(0,180,by=12)), col="black", lty="dotted")