proc(1) = varhac(datin,it1,it2,k1,k2,imax,ilag,imodel,imean,iprint);


/*                      VHGAUSS3.SRC  


                        MARCH 1997

This program is identical to the old vhgauss1.src but allows imax
to be equal to 0

                        VHGAUSS1.SRC

VARHAC ESTIMATOR TO ESTIMATE THE COVARIANCE MATRIX
(SPECTRAL DENSITY AT FREQUENCY ZERO)


This GAUSS procedure is more flexible in choosing lag orders than the
RATS procedure and the FORTRAN subroutine.  That is, in the
regression of the i-th element, the i-th element and the other elements
can enter with a different lag.  The other elements still enter with the
same lag.  This means that if the dimension of the vector of residuals
is equal to 5 and the maximum lag order is equal to 4, then 17 different
models are considered (for each of the 5 residuals).

The input is a data matrix "datin" where each row corresponds with a
different observation


it1:        first observation to be used
it2:        last  observation to be used
k1:         first column to be used
k2:         last  column to be used
imax:       maximum lag order considered  (if imax = 0, then no
            correction for serial correlation will be made)
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 a fixed lag order equal to imax is used
imean:      if equal to 1, then the mean will be subtracted from each series
iprint:     if equal to 1, then the covariance matrix and the chosen lag
            order will be printed on the screen

*/


/*              INITIALIZE PARAMETERS AND MATRICES                   */


local nt,kdim,dat,ave,minpar,minres,minorder,ddd,aic,l,iorder,ex2,ex1,
      dep,cov,bbb,aicnew,ccc,exo1,exo2,iorder1,iorder2,ss1,ss2,tt1,tt2,
      temp1,temp2,npar,iii,t,
      vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat;


nt    = it2 - it1 + 1;
kdim  = k2  - k1  + 1;

dat = datin[it1:it2,k1:k2];

if imean == 1;
ave = meanc(dat);
dat = dat - ave';
endif;


ddd   = dat[imax+1:nt,1:kdim];


if imax == 0;
goto label0;
endif;


minpar   =  zeros(kdim,kdim*imax);
minres   =  ddd';
minorder =  zeros(kdim,2);



aic = ( (nt-imax-1.)/(nt-imax) ) * stdc(ddd)^2.;
aic = ln(aic);


/*     DETERMINE HOW TO ESTIMATE THE PARAMETRIC MODEL    */




if imodel == 3;
goto label2;
endif;

if ilag == 1;
goto label3;
endif;

if ilag == 3;
goto label4;
endif;


/*
        ILAG=2, THUS THE OWN ELEMENT CAN ENTER WITH A DIFFERENT LAG
        FOR EACH L WE FIND THE MODEL THAT MINIMIZES THE MODEL SELECTION
        CRITERION

*/


l = 1;

do while l <= kdim;


/*      step 1: construct dependent and independent variables       */


exo1  = dat[1:nt,l];

    if l == 1;
    exo2  = dat[1:nt,2:kdim];
    endif;

    if l == kdim;
    exo2  = dat[1:nt,1:kdim-1];
    endif;

    if l > 1;
    if l < kdim;
    temp1 = dat[1:nt,1:l-1];
    temp2 = dat[1:nt,l+1:kdim];
    exo2  = temp1~temp2;
    endif;
    endif;


dep    = dat[imax+1:nt,l:l];

iorder1 = 0;
ss1    = exo1[imax+1-1:nt-1,1];


do while iorder1 <= imax;

    if iorder1 >= 2;
    tt1    = exo1[imax+1-iorder1:nt-iorder1,1];
    ss1    = ss1~tt1;
    endif;

iorder2 = 0;
ss2    = exo2[imax+1-1:nt-1,1:kdim-1];



do while iorder2 <= imax;



    if iorder1+iorder2 == 0;
    goto labelx;
    endif;




    if iorder2 >= 2;
    tt2    = exo2[imax+1-iorder2:nt-iorder2,1:kdim-1];
    ss2    = ss2~tt2;
    endif;


    if iorder2 == 0;
    ex1 = ss1;
    endif;

    if iorder1 == 0;
    ex1 = ss2;
    endif;

    if iorder1 > 0;
    if iorder2 > 0;
    ex1 = ss1~ss2;
    endif;
    endif;

/*      step 2: do the regression       */

screen off;
__con   = 0;
_olsres = 1;
{vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat}
        = ols(0,dep,ex1);
screen on;


/*      step 3: calculate the model selection criterion         */

npar = cols(ex1);
aicnew = ( (nt-imax-1.)/(nt-imax) ) * stdc(resid)^2.;

    if imodel == 1;
    aicnew = ln(aicnew) +          2.*npar/(nt-imax);
    else;
    aicnew = ln(aicnew) + ln(nt-imax)*npar/(nt-imax);
    endif;


/*      step 4: if there is improvement store minpar, minres, and minorder */

    if aicnew < aic[l];

    aic[l] = aicnew;

    t = 1;
    do while t <= kdim*imax;
    minpar[l,t] = 0.0;
    t = t + 1;
    endo;

    iorder = 1;
    do while iorder <= iorder1;
    minpar[l,l+(iorder-1)*kdim] = b[iorder,1];
    iorder = iorder + 1;
    endo;

    iorder = 1;
    do while iorder <= iorder2;


    if l == 1;
    minpar[l,2+(iorder-1)*kdim:kdim+(iorder-1)*kdim] =
    b[iorder1+1+(iorder-1)*(kdim-1):iorder1+kdim-1+(iorder-1)*(kdim-1),1]';
    endif;

    if l == kdim;
    minpar[l,1+(iorder-1)*kdim:kdim-1+(iorder-1)*kdim] =
    b[iorder1+1+(iorder-1)*(kdim-1):iorder1+kdim-1+(iorder-1)*(kdim-1),1]';
    endif;

    if l > 1;
    if l < kdim;
    minpar[l,1+(iorder-1)*kdim:(l-1)+(iorder-1)*kdim]  =
    b[iorder1+1+(iorder-1)*(kdim-1):iorder1+(l-1)+(iorder-1)*(kdim-1),1]';
    minpar[l,l+1+(iorder-1)*kdim:kdim+(iorder-1)*kdim] =
    b[iorder1+l+(iorder-1)*(kdim-1):iorder1+kdim-1+(iorder-1)*(kdim-1),1]';
    endif;
    endif;

    iorder = iorder + 1;
    endo;

    minres[l,1:nt-imax] = resid';
    minorder[l,1] = iorder1;
    minorder[l,2] = iorder2;

    endif;

labelx:

iorder2 = iorder2 + 1;
endo;
iorder1 = iorder1 + 1;
endo;

print "l = "  l;
l = l + 1;
endo;

goto label1;

/*
        ILAG=3, THUS ONLY THE OWN ELEMENT ENTERS
        FOR EACH L WE FIND THE MODEL THAT MINIMIZES THE MODEL SELECTION
        CRITERION

*/

label4:

l = 1;

do while l <= kdim;


/*      step 1: construct dependent and independent variables       */


exo1  = dat[1:nt,l];
dep    = dat[imax+1:nt,l:l];

iorder1 = 1;
ss1    = exo1[imax+1-1:nt-1,1];


do while iorder1 <= imax;

    if iorder1 >= 2;
    tt1    = exo1[imax+1-iorder1:nt-iorder1,1];
    ss1    = ss1~tt1;
    endif;

ex1 = ss1;

/*      step 2: do the regression       */

screen off;
__con   = 0;
_olsres = 1;
{vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat}
        = ols(0,dep,ex1);
screen on;


/*      step 3: calculate the model selection criterion         */

npar = cols(ex1);
aicnew = ( (nt-imax-1.)/(nt-imax) ) * stdc(resid)^2.;

    if imodel == 1;
    aicnew = ln(aicnew) +          2.*npar/(nt-imax);
    else;
    aicnew = ln(aicnew) + ln(nt-imax)*npar/(nt-imax);
    endif;


/*      step 4: if there is improvement store minpar, minres, and minorder */

    if aicnew < aic[l];

    aic[l] = aicnew;

    t = 1;
    do while t <= kdim*imax;
    minpar[l,t] = 0.0;
    t = t + 1;
    endo;

    iorder = 1;
    do while iorder <= iorder1;
    minpar[l,l+(iorder-1)*kdim] = b[iorder,1];
    iorder = iorder + 1;
    endo;

    minres[l,1:nt-imax] = resid';
    minorder[l,1] = iorder1;


    endif;

iorder1 = iorder1 + 1;
endo;

print "l"  l;
l = l + 1;
endo;

goto label1;


/*      IMODEL=3, THUS A FIXED LAG FOR ALL ELEMENTS     */

label2:

l = 1;
do while l <= kdim;

iorder = 1;
dep    = dat[imax+1:nt,l:l];
ex1    = dat[imax+1-1:nt-1,1:kdim];

iorder = 2;
do while iorder <= imax;

ex2    = dat[imax+1-iorder:nt-iorder,1:kdim];
ex1    = ex1~ex2;

iorder = iorder + 1;
endo;

screen off;
__con   = 0;
_olsres = 1;
{vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat}
        = ols(0,dep,ex1);
screen on;


    minpar[l,1:kdim*imax] = b';
    minres[l,1:nt-imax] = resid';
    minorder[l,1] = imax;
    minorder[l,2] = imax;

l = l + 1;
endo;

goto label1;


/*      ILAG = 1, THUS ALL VARIABLES ENTER WITH THE SAME LAG     */

label3:

l = 1;
do while l <= kdim;

iorder = 1;
dep    = dat[imax+1:nt,l:l];
ex1    = dat[imax+1-1:nt-1,1:kdim];

do while iorder <= imax;

if iorder >= 2;
ex2    = dat[imax+1-iorder:nt-iorder,1:kdim];
ex1    = ex1~ex2;
endif;

screen off;
__con   = 0;
_olsres = 1;
{vnam,m,b,stb,vc,stderr,sigma,cx,rsq,resid,dwstat}
        = ols(0,dep,ex1);
screen on;

aicnew = ( (nt-imax-1.)/(nt-imax) ) * stdc(resid)^2.;
if imodel == 1;
aicnew = ln(aicnew) +          2.*kdim*iorder/(nt-imax);
else;
aicnew = ln(aicnew) + ln(nt-imax)*kdim*iorder/(nt-imax);
endif;

    if aicnew < aic[l];

    aic[l] = aicnew;
    minpar[l,1:kdim*iorder] = b';
    minres[l,1:nt-imax] = resid';
    minorder[l,1] = iorder;
    endif;


iorder = iorder + 1;
endo;


l = l + 1;
endo;



/*      THE TIME-SERIES REPRESENTATION HAS NOW BEEN ESTIMATED, AND
        WE CAN CALCULATE THE VARHAC ESTIMATOR
*/


label1:

cov = (minres*minres')/(nt-imax);

bbb = eye(kdim);

format /M1 /RD 10,6;
iorder = 1;
do while iorder <= imax;
bbb = bbb - minpar[1:kdim,(iorder-1)*kdim+1:iorder*kdim];
iorder = iorder + 1;
endo;

ccc = inv(bbb)*cov*inv(bbb');

goto labely;


/*      IF NO CORRECTION FOR SERIAL CORRELATION IS REQUIRED, THEN
        THE VARIANCE IS EQUAL TO
*/

label0:

ccc = (ddd'*ddd)/nt;

labely:

/*                      WE ARE DONE             */


format /M1 /RD 10,4;
if iprint == 1;
print " ";
print " ";
print " VARHAC VARIANCE-COVARIANCE MATRIX";

print ccc;
print " ";
print " ";

if imax   == 0;
goto label00;
endif;

if imodel == 1;
print "LAG ORDERS CHOSEN WITH AIC";
print " ";
if ilag <=2;
minorder[.,1:ilag]';
endif;
if ilag == 3;
minorder[.,1]';
endif;
endif;

if imodel == 2;
print "LAG ORDERS CHOSEN WITH SCHWARTZ INFORMATION CRITERION";
print " ";
if ilag <=2;
minorder[.,1:ilag]';
endif;
if ilag == 3;
minorder[.,1]';
endif;
endif;

if imodel == 3;
print "FIXED LAG ORDERS";
print " ";
minorder[.,1]';
endif;

label00:

endif;

retp(ccc);
endp;







