R help - Dec 2010 - GLS with corAR(1) correlation structure residual/standard error calculation

I am using the gls function to fit a two-stage least squares model with
first order autoregressive error terms. Since there is no automated
adjustment for the use of two-stage least squares in this package, I am
trying to manually replicate standard errors of the coefficient estimates in
order to adjust for a first stage OLS estimate of endogenous variables.
However, thus far I have been unable to replicate the residuals or standard
errors produced by this function. My understanding is outlined below, but
using this approach does not yield the reported results. Is anyone familiar
with the inner workings of this function and can either explain the
calculation of the standard errors or provide code that explains the inner
workings of this function.

Thanks!

Example of the model I am running:
model1<- gls(Y~ X1I + X2 + X3 + X4, data=Dat1, correlation = corAR1(),
method = "ML")

My understanding of model errors:
Y = b_0 + X1 b_1+ ...Xk b_k + Z
Z_t =phi Z_{t-1) + e_t

The residuals reported by GLS are the Z's, while the white noise terms are
the e's. I cannot replicate the reported residuals using this approach. I
also do not know how Z_0 should be calculated, i.e. what does the first step
of this recursive procedure look like?
>From the residuals, I also cannot replicate the reported standard errors. I
am using se(b_j) = sqrt(sigma^2/sum(x_i-x_mean)^2) where sigma =sqrt(SSR/df)

Any help on this or explanation of how GLS works would be much appreciated.

Any clarification would be much appreciated.

	[[alternative HTML version deleted]]