function results = ols_ts(y,X,info);

%OLS_TS is my version of OLS, including a scaling of data to avoid
%instabilities. It also uses my version of LeSage's ols (with tinv)
%and my version of prt_reg (prt_reg_ts).
%
% Written by: T.E.SMITH  01/20/03  (modified: 01/28/08) 
%
%INPUTS: y = n-vector dependent variable data
%        X = nxk vector of explanatory variables 
%            (constant is included automatically)
%    info  = (optional) structure of additional options
%
%             info.names = names of variables, for example
%                          strvcat('y','var1',var2').
%             info.noint = 1 if no intercept is to be used
%             info.dw    = 1 if Durban-Watson test is to be performed
%
%
%OUTPUTS: a structure
%        results.meth  = 'ols'
%        results.beta  = bhat     (nvar x 1)
%        results.tstat = t-stats  (nvar x 1)
%        results.yhat  = yhat     (nobs x 1)
%        results.resid = residuals (nobs x 1)
%        results.sige  = e'*e/(n-k)   scalar
%        results.rsqr  = R-squared     scalar
%        results.rbar  = adjusted R-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


[n k] = size(X);

names_flag = 0;

int_flag = 0;

dw_flag = 0;

if nargin == 3 % parse user input options
    
    fields = fieldnames(info);
    
    nf = length(fields);
    
    for i=1:nf
        
        if strcmp(fields{i},'names') 
            
            vnames = info.names;
            
            names_flag = 1;
            
        elseif strcmp(fields{i},'noint')            
            
            int_flag = 1;
            
        elseif strcmp(fields{i},'dw')
            
            dw_flag = 1;           
        
        end;
        
   end;
   
end;

if int_flag == 0;  %Make intercept
    
	X = [ones(n,1),X];
    
end

%Now rescale data to avoid possible instabilities

S = inv(diag(max(abs(X)))) ;

x = X*S ;

%Now run LeSage's ols using (y,x):

[nobs nvar] = size(x); 

nobs2 = size(y,1);

 if (nobs ~= nobs2); 
     
     error('x and y must have same # obs in ols'); 
     
 end

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

[q r] = qr(x,0);

xpxi = (r'*r)\eye(nvar);

results.beta = xpxi*(x'*y);

bb = results.beta;

b = S*bb ; %This factors out the rescaling effect as seen below:

%EXPLANATION: Z = X*S => bb = inv(Z'*Z)*Z'*y
%                              = inv[S'*(X)'*(X)*S]*S'*(X)'*y
%                              = inv(S)*inv[(X)'*(X)]*inv(S')*S'*(X)*y
%                              = inv(S)*{inv[(X)'(X)]*(X)'*y}
%                              = inv(S)*b
%             => b = S*bb.

results.beta = b;

results.yhat = X*results.beta;
results.resid = y - results.yhat;
sigu = results.resid'*results.resid;
results.sige = sigu/(nobs-nvar);
tmp = (results.sige)*(diag(xpxi));
sigb= S*sqrt(tmp); %must rescale sigb the same as beta.
tcrit=-tinv(.025,nobs-nvar);
results.bint=[results.beta-tcrit.*sigb, results.beta+tcrit.*sigb];
results.tstat = results.beta./sigb;
ym = y - mean(y);
rsqr1 = sigu;
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)/sigu; % durbin-watson

%LeSage ols done.

% results.yhat = X*b;
% 
% e = y - X*b ;
% 
% results.resids = e;
% 
% results.sige = e'*e/(n-k-1);

%Construct Variable Names

if names_flag == 0;      
       
    if int_flag == 1;
        
        vnames = 'Y';
        
    else
        
        vnames = strvcat('Y','constant');
        
    end
    
    nvar = k; 
    
    for i=1:nvar

        tmp = ['Variable ',num2str(i)];

        vnames = strvcat(vnames,tmp);

    end;
    
else;

    if int_flag == 0;

        names = vnames(2:k+1,:);

        vnames = strvcat(vnames(1,:),'constant',names);

    end
    
end

%Log Likelihood

e = results.resid;

var = e'*e/n;

L = -(n/2)*(log(2*pi*var) + 1);

%Compute corrected AIC (Burnham and Anderson, MODEL SELECTION .. pp.12,66

K = k + 2; %Number of parameters estimated

AIC = 2*(-L + K);

AIC_c = AIC + 2*K*(K+1)/(n-K-1);

results.AIC = AIC;

results.AIC_c = AIC_c;

%prt(results,vnames);

%Finally, run appropriate version of LeSages's prt:

prt_reg_ts(results,vnames);

if dw_flag == 1
    
    num = 0;
    
    for t = 2:n
        
        num = num + e(t)*e(t-1);
        
    end
    
    rho_hat = num/(e'*e);
    
    [pval,dw] = dwtest(e,X);
    
    disp('DURBAN-WATSON TEST:')    
    disp(['DW = ',num2str(dw)]);
    disp(['rho_hat = ',num2str(rho_hat)]);
    disp(['P-Value = ',num2str(pval)]);
    disp(' ');
    
end