function results=gls(y,X,omega)
%PURPOSE: generalized least-squares regression (by default with AR(1))
%---------------------------------------------------
% USAGE: results = gls(y,X,omega)
% where: y = dependent variable vector    (nobs x 1)
%        x = independent variables matrix (nobs x nvar)
%        omega = optional - positive definite symmetric matrix (nobs x nobs)
%           by default - omega for AR(1) process
%---------------------------------------------------
% RETURNS: a structure
%         results.meth  = 'gls'
%         results.beta  = bhat     (nvar x 1)
%         results.tstat = t-stats  (nvar x 1)
%         results.bstd  = std deviations for bhat (nvar x 1)
%         results.yhat  = yhat     (nobs x 1)
%         results.resid = residuals (nobs x 1)
%         results.sige  = e'*e/(n-k)   scalar
%         results.rsqr  = rsquared     scalar
%         results.rbar  = rbar-squared scalar
%         results.dw    = Durbin-Watson Statistic
%         results.nobs  = nobs
%         results.nvar  = nvars
%         results.y     = y data vector (nobs x 1)
%         results.bint  = (nvar x2 ) vector with 95% confidence intervals on beta


nIn = nargin;

if nIn < 2
  error('Not enough input arguments.');
elseif nIn > 3
  error('Too many input arguments.');
end

[nobs,nvar] = size(X);

if nIn == 2
   res_pom = ols(y,X);
   u = res_pom.resid(2:end);
   u_1 = res_pom.resid(1:end-1);
   res_u = ols(u,u_1);
   Rho = res_u.beta;
   for i = 1:nobs
        omega(i,1)=Rho^(i-1);
   end
   for j = 2:nobs
        omega(:,j)=omega(:,j-1)./Rho;
   end
        omega = tril(omega)+triu(omega',1); 
        omega = omega/(1-Rho^2);
end

P = chol(inv(omega));

results.meth = 'gls';
results.y = y;
results.nobs = nobs;
results.nvar = nvar;

xpxi = inv(X'*(inv(omega))*X);
results.beta = xpxi*(X'*inv(omega)*y);

results.yhat = X*results.beta;
results.pyhat = P*X*results.beta;

results.resid = P*y - results.pyhat;
%sigu = results.resid'*inv(omega)*results.resid;
sigu = results.resid'*results.resid;
results.sige = sigu/(nobs-nvar);
results.bstd = sqrt(diag(results.sige*xpxi));
tcrit=-tdis_inv(.025,nobs);
results.bint=[results.beta-tcrit.*results.bstd, results.beta+tcrit.*results.bstd];
results.tstat = results.beta./results.bstd;
ym = y - mean(y);
rsqr1 = (y-results.yhat)'*(y-results.yhat);
rsqr2 = ym'*ym;
results.rsqr = 1.0 - rsqr1/rsqr2; % r-squared
rsqr1 = rsqr1/(nobs-nvar);
rsqr2 = rsqr2/(nobs-1.0);
if rsqr2 ~= 0
results.rbar = 1 - (rsqr1/rsqr2); % rbar-squared
else
    results.rbar = results.rsqr;
end;
ediff = results.resid(2:nobs) - results.resid(1:nobs-1);
results.dw = (ediff'*ediff)/(results.resid'*results.resid); % durbin-watson
