/*  LS.PRO  PROGRAM TO CALCULATE LEAST-SQUARES COEFFICIENTS,
    REGULAR STANDARD ERRORS AND VARHAC STANDARD ERRORS   (MARCH 1997)

    BOTH STANDARD ERRORS CORRECT FOR DEGREES OF FREEDOM

    The file with data used in this example is added at the end

    Implementing this program would require the following lines
    in gauss

    > run ols.src
    > run vhgauss3.src
    > run ls.pro

    Output is written to the file exam.out

    This program assumes that the user also runs the GAUSS OLS procedure.
    This is not necessary, however, if one just wants to get
    parameter estimates and standard errors.

    PARAMETER SECTION:

*/

ldim  =   5;
nobs  =  1000;
imax   = 4;
ilag   = 1;
imodel = 2;
iprint = 1;

/*
ldim:   the number of independent variables (including the constant if there
        is one)
nobs:   number of observations

the next four options deal with the varhac procedure

imax:   maximum lag order considered
ilag:   if equal to 1, then all elements enter with the same lag
        if equal to 2, then the own lag can enter with a different
        lag
        if equal to 3, then only the own lag enters
imodel: if equal to 1, then AIC is used
        if equal to 2, then SCHWARTZ is used
        if equal to 3, then fixed lag order imax is used
iprint: if equal to 1, then if the covariance matrix and the chosen lag
        order will be printed on the screen


        DATA SECTION:

        The program assumes that the dependent variable is
        in the matrix y with dimensions (nobs x 1) and that the
        independent variables are in the (nobs x ldim) matrix x.

        On my web page I have included two example data files.  In both
        cases, the generating process is equal to

        y = 0.1 + 0.2*x1 + 0.3*x3 + 0.4*x4 + 0.5*x5 + u

        In exam2.dat, the independent variables and the error term are
        i.i.d. random variables and the standard formula for the
        standard errors inv(x'x)*sigma^2 is correct

        In exam1.dat, x1 and u are serially correlated, which means that
        the standard formula is not correct for the first two standard
        errors.

        Note that even when VARHAC chooses zero lag orders, the VARHAC
        standard errors still corrects for heteroskedasticity and are,
        therefore, in general, not equal to the uncorrected standard errors.



*/

load datin[nobs,ldim] = exam2.dat;

x2 = datin[.,2:ldim];
y  = datin[.,1];
x1 = ones(nobs,1);
x  = x1~x2;

/*      The user of this program could use the GAUSS OLS procedure
        to check the answers of this program and to get some more
        standard statistics.  However, obtaining parameter estimates
        and standard errors does not require doing this.
*/

output file = exam.out reset;

print "RESULTS FROM THE GAUSS OLS PROCEDURE:";
print " ";

call ols(0,y,x);

output off;

/*      BEGIN OF THE PROGRAM,  THE USER DOES NOT HAVE TO CHANGE
        ANYTHING BELOW

*/

beta = inv(x'*x)*x'*y;
u  = y - x*beta;
r  = ones(nobs,ldim);

jdim = 1;
do while jdim <=ldim;
r[.,jdim] = u.*x[.,jdim];
jdim = jdim + 1;
endo;

datin = r;


it1  = 1;
it2  = nobs;
k1   = 1;
k2   = ldim;
imean = 0;
ccc = varhac(datin,it1,it2,k1,k2,imax,ilag,imodel,imean,iprint);

format /M1 /RD 10,4;

xpx = inv(x'*x);
cov = xpx*ccc*xpx;

cova = xpx*(u'*u)/(nobs-ldim);

se  = ones(ldim,1);
sea = ones(ldim,1);

jj = 1;
do while jj <=ldim;
se[jj,1]  = (cov[jj,jj]*nobs*nobs/(nobs-ldim))^0.5;
sea[jj,1] = (cova[jj,jj])^0.5;
jj = jj+1;
endo;

output file = exam.out on;

print " ";
print " ";
print " RESULTS FROM THE LS.PRO PROGRAM:";
print " ";
print "   LS COEF       SE      VARHAC SE";
results = beta~sea~se;
results;

end;