# How to calculate the regression variance for a GLS model?

I need to calculate the regression variance ($\sigma^2$) in order to estimate both the confidence intervals and the prediction intervals in a gls regression analysis. For the analysis, the covariance matrix ($V$) of the response variable ($y$) is known in advance, and so I use it directly as the weighting matrix (=$V^{-1}$) in the gls regression analysis.

The regression variance is a weighted sum of the residual error: $\sigma^2 = \frac{ (Y – X\beta)^T C^{-1} (Y – X\beta)}{n – p}$

My question/problem is how to determine the weighting matrix $C^{-1}$? $C$ cannot be set equal to $V$ since (according to the above equation) $C$ must be dimensionless while $V$ has the same units as $\sigma^2$.

Based on my reading of the literature and available texts, it seems that $C$ is the correlation matrix and is a scaled or normalized form of the covariance matrix $V$. i.e., $V = Var(\epsilon^2) = \sigma^2 C$. But my problem is that $\sigma^2$ is not yet known, and so I need another way find $C$ from $V$.

R functions such as gls() will compute the regression variance (if I knew how gls() does this, it would answer my question). However I cannot use gls() in this case since I am specifying a user-defined covariance (weighting) matrix, and gls() only accepts a limited set of specific correlation structures.

In fact a possible solution can be found in this earlier post where an equation for the SEE (or sigma2) for a GLS regression was cited :

GLS calc of SEE: sqrt( sum( ( residuals from linear model) ^ 2 * glsWeight ) ) / sum( glsWeight ) * length( glsWeight ) / residualDegreeFreedom )

However I am unable to ascertain the validity of this equation and cannot find its source reference.

You are a bit off I think because the formula for epsilon you seem to use relies on the neccesity to decompose your error variance into a constant term ($\sigma^2$) and a correlation term $\Omega$, which is the $C$ in your post. You actually have to use this $C$ as a weight for the regression. Let me explain.

It is assumed you know $\Omega$, otherwise you have to do FGLS instead. The formula you have is a result of a transformed model. First, let's look at our GLS estimators:
Under the assumption that $V[\epsilon|X]=\sigma^2\Omega$ you can transform the model and estimate the following $$\widehat{\beta}^{GLS} = [X'\Omega^{-1}X]^{-1}X'\Omega^{-1}y$$ This, generally, is a BLUE estimator under the assumptions.
Which has, under the above and usual assumptions this variance $$V[\widehat{\beta}^{GLS}] = \sigma^2[X'\Omega^{-1}X]^{-1}$$ Since you know $\Omega$, you can now calculate everything and estimate the $\sigma^2$. To understand how have a look at how the GLS estimator is reached.
You reached this result by applying a matrix $P$ to the model formula as such $$PY = PX\beta + P\epsilon$$ What is P? If you know omega, you can calculate the $P$ (Cholesky decomposition) with $$P'P = \Omega^{-1}$$ And then you have a consistent estimator in $$\widehat{\sigma}^2 = \frac{(P\hat{\epsilon}^{OLS})'P\hat{\epsilon}^{OLS}}{n-k} \\ = \frac{ \hat{\epsilon}'^{OLS} \Omega^{-1} \hat{\epsilon}^{OLS}}{n-k}$$ This is the formula you have listed, but it implies a possible decomposition $V[\epsilon|X]=\sigma^2\Omega$. Note how it uses the residuals of OLS. This is because of said existence of a transformed and non-transformed model.

You said you used $V[Y]=V$ as weight for GLS. Can you tell us the exact formula? Is it $$\widehat{\beta}^{GLS} = [X'V^{-1}X]^{-1}X'V^{-1}y$$ ?
In that case you used $V$ instead of $C$. You mentioned that you think $V[\epsilon]=V$ and $V=\sigma^2C$ (I assume this is what you mean).
You could go ahead and plug in $V$ instead of $C$. But of course this is not actually correct because $V$ is not equal to $C$ as you discovered.
This seems to be your problem.

If you know nothing about the composition of your error terms, you just can't get to that decomposition and the $C$ you need. That is the 'point' of GLS. It takes advantage of a known correlation structure $\Omega$ or $C$ to regain BLUE estimators.
So the equation you can not validate is the result of said transformed model. For this you also need to have your $\Omega$ or $C$ as you called it and by that the $P$ to transform the model.
If you do not have this, use FGLS instead.

• You could go ahead and plug in $V$ instead of $C$. But of course this is not actually correct because $V$ is not equal to $C$ as you discovered … In fact this is not true. In the equation for $\beta^GLS$ any scaling factor applied to $V$ (in the equation that you showed) will cancel. Thus $V$ or $C$ or any $V$ scaled by any constant will work and yield the same answer for $\beta$. In my case, I used $V$ to compute $\beta$, in the absence of knowing $C$ and $\sigma^2$. – treemake Apr 19 '13 at 12:05
• If you know nothing about the composition of your error terms, you just can't get to that decomposition and the C you need. That is the 'point' of GLS … For my particular problem I know $V$ but not $C$. This is because my empirical measurements of my response variable $y$ include estimates of their variances and covariances. (It is somewhat unusual to know $V$ a priori and not $C$, but nevertheless quite possible, if $y$ represents instrumental measurements or values taken from a complex survey design for example). – treemake Apr 19 '13 at 12:09
• And then you have a consistent estimator in … I would like to clarify a point that you made the first part of your response (before being able to evaluate it). The equation for $\hat{\sigma }^2$ is expressed in terms of $\hat{\epsilon}^{OLS}$. I do not believe this is correct. In fact $\hat{\epsilon}^{OLS}$ should be replaced by simply $\hat{\epsilon}$, where the latter are in fact the GLS residuals. This follows from examination of the preceding equations where $\epsilon$ represents $\epsilon^{GLS}$. (In fact $\epsilon^{OLS} = P\epsilon$) – treemake Apr 19 '13 at 12:10
• I would very be interested in your follow-up to my comments ... – treemake Apr 19 '13 at 12:14
• Post 1: It is true that the unbiasedness is not in question when using any weighting factor. It is not relevant to your question because you are looking for the covariances. The GLS formula in question results from the transformed model with the $P$ – IMA Apr 19 '13 at 12:29