# Non-Correlated errors from Generalized Least Square model (GLS)

As a financial institution, we often run into analysis of time series data. A lot of times we end up doing regression using time series variables. As this happens, we often encounter residuals with time series structure that violates basic assumption of independent errors in OLS regression. Recently we are building another model in which I believe we have regression with autocorrelated errors.The residuals from linear model have lm(object) has clearly a AR(1) structure, as evident from ACF and PACF. I took two different approaches, the first one was obviously to fit the model using Generalized least squares gls() in R. My expectation was that the residuals from gls(object) would be a white noise (independent errors). But the residuals from gls(object) still have the same ARIMA structure as in the ordinary regression. Unfortunately there is something wrong in what I am doing that I could not figure out. Hence I decided to manually adjust the regression coefficients from the linear model (OLS estimates). Surprisingly that seems to be working when I plotted the residuals from adjusted regression (the residuals are white noise). I really want to use gls() in nlme package so that coding will be lot simpler and easier. What would be the approach I should take here? Am I supposed to use REML? or is my expectation of non-correlated residuals (white noise) from gls() object wrong?

gls.bk_ai <- gls(PRNP_BK_actINV ~ PRM_BK_INV_ENDING + NPRM_BK_INV_ENDING,
correlation=corARMA(p=1), method='ML',  data  = fit.cap01A)

gls2.bk_ai  <- update(gls.bk_ai, correlation = corARMA(p=2))

gls3.bk_ai <- update(gls.bk_ai, correlation = corARMA(p=3))

gls0.bk_ai <- update(gls.bk_ai, correlation = NULL)

anova(gls.bk_ai, gls2.bk_ai, gls3.bk_ai, gls0.bk_ai)
## looking at the AIC value, gls model with AR(1) will be the best bet

acf2(residuals(gls.bk_ai)) # residuals are not white noise


Is there something wrong with what I am doing???????

The residuals from gls will indeed have the same autocorrelation structure, but that does not mean the coefficient estimates and their standard errors have not been adjusted appropriately. (There's obviously no requirement that $\Omega$ be diagonal, either.) This is because the residuals are defined as $e = Y - X\hat{\beta}^{\text{GLS}}$. If the covariance matrix of $e$ was equal to $\sigma^2\text{I}$, there would be no need to use GLS!

In short, you haven't done anything wrong, there's no need to adjust the residuals, and the routines are all working correctly.

predict.gls does take the structure of the covariance matrix into account when forming standard errors of the prediction vector. However, it doesn't have the convenient "predict a few observations ahead" feature of predict.Arima, which takes into account the relevant residuals at the end of the data series and the structure of the residuals when generating predicted values. arima has the ability to incorporate a matrix of predictors in the estimation, and if you're interested in prediction a few steps ahead, it may be a better choice.

EDIT: Prompted by a comment from Michael Chernick (+1), I'm adding an example comparing GLS with ARMAX (arima) results, showing that coefficient estimates, log likelihoods, etc. all come out the same, at least to four decimal places (a reasonable degree of accuracy given that two different algorithms are used):

# Generating data
eta <- rnorm(5000)
for (j in 2:5000) eta[j] <- eta[j] + 0.4*eta[j-1]
e <- eta[4001:5000]
x <- rnorm(1000)
y <- x + e

> summary(gls(y~x, correlation=corARMA(p=1), method='ML'))
Generalized least squares fit by maximum likelihood
Model: y ~ x
Data: NULL
AIC      BIC    logLik
2833.377 2853.008 -1412.688

Correlation Structure: AR(1)
Formula: ~1
Parameter estimate(s):
Phi
0.4229375

Coefficients:
Value  Std.Error  t-value p-value
(Intercept) -0.0375764 0.05448021 -0.68973  0.4905
x            0.9730496 0.03011741 32.30854  0.0000

Correlation:
(Intr)
x -0.022

Standardized residuals:
Min          Q1         Med          Q3         Max
-2.97562731 -0.65969048  0.01350339  0.70718362  3.32913451

Residual standard error: 1.096575
Degrees of freedom: 1000 total; 998 residual
>
> arima(y, order=c(1,0,0), xreg=x)

Call:
arima(x = y, order = c(1, 0, 0), xreg = x)

Coefficients:
ar1  intercept       x
0.4229    -0.0376  0.9730
s.e.  0.0287     0.0544  0.0301

sigma^2 estimated as 0.9874:  log likelihood = -1412.69,  aic = 2833.38


EDIT: Prompted by a comment from anand (OP), here's a comparison of predictions from gls and arima with the same basic data structure as above and some extraneous output lines removed:

df.est <- data.frame(list(y = y[1:995], x=x[1:995]))
df.pred <- data.frame(list(y=NA, x=x[996:1000]))

model.gls <- gls(y~x, correlation=corARMA(p=1), method='ML', data=df.est)
model.armax <- arima(df.est$y, order=c(1,0,0), xreg=df.est$x)

> predict(model.gls, newdata=df.pred)
 -0.3451556 -1.5085599  0.8999332  0.1125310  1.0966663

> predict(model.armax, n.ahead=5, newxreg=df.pred$x)$pred
 -0.79666213 -1.70825775  0.81159072  0.07344052  1.07935410


As we can see, the predicted values are different, although they are converging as we move farther into the future. This is because gls doesn't treat the data as a time series and take the specific value of the residual at observation 995 into account when forming predictions, but arima does. The effect of the residual at obs. 995 decreases as the forecast horizon increases, leading to the convergence of predicted values.

Consequently, for short-term predictions of time series data, arima will be better.

• Applying an arma structure to the residuals would I think be a little different from the result the gls gives regarding the regression parameters. – Michael Chernick Jul 21 '12 at 16:50
• Also the ARMA modeling approach fit the correlation structure in the residuals rather than specifying it in the covariance matrix. – Michael Chernick Jul 21 '12 at 16:57
• Bowman, thank you sir for such a concise and clear explanation. So, in short using residuals structure at the end of data series instead of Covariance matrix will be a better approach in prediction. What that means is predict.arima() will give you better prediction than predict.gls() correct? – Anand Jul 22 '12 at 15:30
• I've appended some prediction example stuff to the response, but the short version is yes, for time series data predict.arima() will give you better prediction than predict.gls(). – jbowman Jul 22 '12 at 16:12
• jBowman, may I asked two follow-ups to this rather old question (I was just logging in to ask the same as this OP)? 1) I assume the same holds for heteroskedasticity - i.e. the residuals from GLS still wont look homoscedastic? 2) Is the situation different for regression with arima errors (your arima example)? I have always read that to see if you have the model specified properly, the residuals should be white noise? – B_Miner Apr 14 '14 at 2:04

You want the normalized residuals. See ?residuals.lme.

#Reproducible code from ?corARMA
fm1Ovar.lme <- lme(follicles ~ sin(2*pi*Time) + cos(2*pi*Time),
data = Ovary, random = pdDiag(~sin(2*pi*Time)))
fm5Ovar.lme <- update(fm1Ovar.lme,
corr = corARMA(p = 1, q = 1))

#raw residuals divided by the corresponding standard errors
acf(residuals(fm5Ovar.lme),type="partial")

#standardized residuals pre-multiplied
#by the inverse square-root factor of the estimated error correlation matrix
acf(residuals(fm5Ovar.lme,type="normalized"),type="partial")