Seasonal and autocorrelated regression residuals: difference raw data or residuals?

I have a multiple linear model for time series data for which the regression residuals are autocorrelated and display seasonal behavior. This seasonal behavior is induced deterministically by a cyclic variable written into the model. In order to calculate corrected standard errors for the regression coefficients, I intend to use generalized least squares with correction for the autocorrelation in the residuals.

But for the seasonality, I am not sure for which time series I should perform the differencing: If for raw original data (observations), or for the regression residuals. For the data, I could superimpose each cycle and take the mean, thus removing season effects, but this would alter the data, and also would generalize bad for other kinds of data. For the regression residuals, I dont know if I have to difference the residuals time series substracting for each period, or to handle it as some kind of seasonal ARMA process.

Any hints?

LAST EDIT: The problem behind this question might have been resolved already, and might have been product of a misconception.

Regression was being done over simulations of the model. These simulations did not contain any stochastic error factor, so the simulations were purely of a deterministic nature. Regression errors were just showing an unperfect fit to the data and the regression residuals were obviously following a deterministic pattern. This pattern was not a result of some neglected explanatory variable in the model, and thus there was no reason to model it (by means of an ARMA model or any kind). When adding white noise to the data -which should be a crucial step on any simulation of real data- regression residuals were mostly dominated by stochasticity, loosing the autocorrelated behavior.

In order to calculate corrected standard errors for the regression coefficients, I intend to use generalized least squares with correction for the autocorrelation in the residuals.

Note that generalized least squares (GLS) would affect not only the standard errors but also the point estimates. Anyhow, you could gain power by estimating regression with an explicitly specified error structure, e.g. regression with ARMA errors as can be done using functions stats::arima or forecast::auto.arima in R. There you use maximum likelihood estimation instead of GLS. See related blog posts by Francis X. Diebold "The HAC Emperor has no Clothes" and "The HAC Emperor has no Clothes: Part 2" where he encourages explicit error specification as a way to get better coefficient estimates and gain predictive power. Although he discusses the case of HAC there, I believe similar conclusions apply here, too.

But for the seasonality, I am not sure for which time series I should perform the differencing: If for raw original data (observations), or for the regression residuals.

Since the problem arises due to a cyclic regressor, you could remove the deterministic component of the cyclic variable before including it in the model, or alternatively you could include some seasonal terms (dummies or Fourier terms) in the model.

For the data, I could superimpose each cycle and take the mean, thus removing season effects, but this would alter the data, and also would generalize bad for other kinds of data.

I am a little confused here, but I will try addressing this nevertheless.
With regards to the regressor, you can adjust using a model, and so altering data is not really a problem because you keep track of how you did it and you can recreate the original variable if you need to.
Regarding generalization, if the cyclic behaviour is unique for this instance, keeping it untreated would not help. If, on the other hand, it is similar across this instance and the ones you want to generalize to, you would not lose by removing the deterministic component before running the regression but then using it to adjust the other cases similarly.

A technical note: If you are doing a regression with ARMA errors, then it is the error that gets differenced. If the errors is some SARIMA process, regular treatment of SARIMA models applies (roughly speaking, you do not have to worry that it is a regression error rather than raw data).

• I had never heard of kernel HAC and I am not sure that is the approach I am taking. Maybe I expressed it wrong in the question, but the approach I would like to take is a full model specification and parameterization of the residual time-series. For the estimation I would use GLS for estimation of the parameters with inclusion of SARIMA residuals (that is what I meant with "correction for autocorrelation"). The last paragraph answers my question. But is there a reason why I shouldnt remove seasonality from the data thus then only treating the residuals as an AR model? – hirschme Jul 29 '16 at 10:33
• After searching for a while I found this en.wikipedia.org/wiki/Newey%E2%80%93West_estimator , which I think is what you meant by kernel-HAC and yeah it is not the approach I am taking. I do want to parameterize the residual time series as an autoregressive model. The main problem is if to adress the seasonality in the raw data and then use an AR model, or if to leave it and then difference the residuals by/or using a SARIMA model (I guess the I = Integrated would come because of the differencing done before, and the Moving Average because I somehow have never seen a "SARI" model.) – hirschme Jul 29 '16 at 12:06
• Thanks for the updates. Let me come back to it later in the day. – Richard Hardy Jul 29 '16 at 12:16
• @RichardHardy Thank you for the anwer. I will update soon with an example. I was intrigued by a couple of points in the answer. First, I have never heard GLS changing the point estimates, just the SE. And the optimization would be max. Likelihood too, where the likelihood function would be adapted for the error model (and thus be other than a gaussian as with OLS). With a proper error specification it could even be transformed into a linear regression with gaussian errors and max.Like = OLS, then transforming the coefficients back (for ex. by Cochrane-orcutt estimation in the case of AR(1) ). – hirschme Jul 30 '16 at 16:09
• @hirschme, GLS not changing the point estimates was actually what was happening in STATA for a long time, not sure if they corrected the mistake. But you can see immediately that the estimates differ: $\beta_{OLS}=(X'X)^{-1}X'y \not \equiv \beta_{GLS}=(X'WX)^{-1}X'Wy$ for some weighting matrix $W$. – Richard Hardy Jul 30 '16 at 16:12