PROCEDURE VARHAC start end 

***********************************************************************
***
***	   MAY 17, 1996
***
***    DOS VERSION
***
***	   This program calculates the spectral density at frequency zero
***	   of the input series.  The program assumes that the input series
***	   are residuals, that is, no mean is taken out. 
***	   The spectrum is calculated by fitting a VAR to the model, and a
***	   model selection criterion is used to find the best model.
***	   The minimum lag considered is one.  The GAUSS and FORTRAN
***	   routines allow for a zero lag lenght. 
***
***
***    To use interactively (i.e., menues):
*** 
***    @VARHAC start end
***    # the dependent variables
***
***    NOTE THAT YOU ANSWER QUESTIONS BY TYPING THE ITEM NUMBER!!!
***
***
***    To use noninteractively:
***
***    @VARHAC(nointeractive,det= ,crit= ,maxlag= ) start end
***    # the dependent variables
***
***    where:
***        
***    det  =  "none" or "constant"    --  deterministic component in the VAR
***    crit =  "AIC", "SIC" or "Kbar"  --  model selection criteria
***        (Kbar means use Kbar as the lag in each eq. in the VAR)
***    maxlag =  maxlag to consider for the model selection criteria
***
***
***    For interactive and non-interactive use NOPRINT
***	   will suppress printed output
***
***
***    The Spectral density estimator is saved to the symmetric matrix SDM  
***
***	   Note that it does not make much difference whether you include
***	   a constant in the regressions since the mean of the full 
***    series is (supposed to be) zero.  By lagging variables, this
***    may no longer be exactly true for the series used in the regression
***    in which case including a constant makes a small difference
***
***
***    EXAMPLES
***
***    source(noecho) VARHAC.SRC
***
***    @varhac(nointeractive,det=constant,crit=aic,maxlag=4) 59:1 93:4 
***    #u1 u2 u3 u4 u5
***
***    or
***
***    @varhac(nointeractive,noprint,det=none,crit=aic,maxlag=4) 59:1 93:4
***    #u1 u2 u3 u4 u5
***
***
***
***    !!!  ALL VARIABLES IN THE VARHAC PROCEDURE ARE LOCAL.  HOWEVER,
***    THE VARHAC PROCEDURE ESTIMATES A VAR SYSTEM, AND THE PARAMETERS
***    ESTIMATED BY THIS SYSTEM ARE NOT.
***
***
***
***		REFERENCES
***
*** 	DenHaan, Wouter J. and Andrew T. Levin, Inferences from parametric
***		and non-parametric covariance matrix estimation procedures, 
***		NBER technical working paper #195.
***
*** 	DenHaan, Wouter J. and Andrew T. Levin, A practioner's guide to
***		robust covariance matrix estimation, NBER technical working paper #197.


*********************************************************************

DECLARE SYMM VCOVMAT SDM

TYPE INTEGER  start end

OPTION CHOICE    DET  1  NONE CONSTANT  
OPTION INTEGER   MAXLAG 4
OPTION CHOICE    CRIT 1  AIC SIC Kbar
OPTION SWITCH    INTERACTIVE 1
OPTION SWITCH    PRINT 1

LOCAL INDEX    varlist levsx ccmat  laglevsx 
LOCAL INTEGER  lag vars counter startl endl nprint num deter critmax
LOCAL INTEGER  Kbar Kmax ii jj kk i j k nvars msc 
LOCAL RECT     lev laglev f11 F FP 
LOCAL VECT[series]  DEPVAR RESIDVEC RESIDS
LOCAL RECT     cmat  AIC SIC  U  Y  YLAGS 
LOCAL SYMM     VCOV VCORR SDMAT
LOCAL VECT     bb bb2 STDEV RSSMAT 
LOCAL SERIES   AICVEC SICVEC 
LOCAL VECT     AICMINENT SICMINENT AICMIN SICMIN
LOCAL VECT[rect]    outmatrix
LOCAL VECT[labels]  vlabel
LOCAL VECT[rect]    VMAMAT COEFFMAT
LOCAL VECT[INDEX]   laglevsx2
LOCAL VECT[VECT]    COEFFVEC 
  
ENTER(varying) varlist
COMPUTE vars = %rows(varlist)

INQUIRE(regressorlist)  startl>>start  endl>>end
#varlist

DIM vlabel(vars)
DO ii = 1,vars
   COMPUTE vlabel(ii) = %label([series] varlist(ii)) 
END DO 

IF PRINT.eq.1 {
 DISP; DISP  
 DISP '*******************************************************'
 DISP '* The VARHAC Estimator of the Spectral Density Matrix *'
 DISP '*     by Norman Morin                                 *'
 DISP '*******************************************************' 
 DISP 
}

DIMENSION depvar(vars)
DIMENSION residvec(vars) 
DIMENSION coeffvec(vars)
DIMENSION levsx(0) 
DIMENSION laglevsx(0)
DIMENSION ccmat(0)

IF INTERACTIVE.eq.0 {
   COMPUTE Kbar = maxlag
   COMPUTE deter = det
   COMPUT MSC = crit
 }

ELSE { 
   MENU 'Deterministic components'
      CHOICE 'None'
          COMPUTE deter = 1
      CHOICE 'Constant'
          COMPUTE deter = 2
   END MENU
 
   DISPLAY(store=qstring) $
     'Input the maximum lag to consider for AIC/Schwarz' 
   QUERY(prompt=qstring,status=check) Kbar

   MENU 'Model Selection Criterion'
      CHOICE 'AIC'
          COMPUTE MSC = 1
      CHOICE 'Schwarz'
          COMPUTE MSC = 2
      CHOICE 'The maximum lag previously input'
          COMPUTE MSC = 3
   END MENU
   }

DO i = 1,vars
    ENTER(varying) levsx
    # levsx varlist(i) 
    ENTER(varying) laglevsx 
    # laglevsx varlist(i){1 to Kbar}
END DO

MAKE(trans) lev startl endl
# levsx
 MAKE(trans) laglev startl endl
# laglevsx

DO ii = 1,vars
   SET depvar(ii) startl endl = lev(ii,t-startl+1)
   LABLES DEPVAR(ii)
   #VLABEL(ii) 
END DO  
 

IF MSC.lt.3 { *********************************************
IF (deter.eq.1)
    {
    DIM aic(vars,Kbar)
    COMPUTE aic = %CONST(0.) 
    DIM sic(vars,Kbar)
    COMPUTE sic = %CONST(0.) 
    DO lagnum = 1,Kbar
        DIM RESIDVEC(vars)
        
        SYSTEM 1 to vars
        LAGS 1 to lagnum  
        VARIABLES depvar
        END(SYSTEM)
        ESTIMATE(noftest,noprint) startl+Kbar endl RESIDVEC(1)
     
        DO jj = 1,vars
           set rr1 startl+Kbar endl = RESIDVEC(jj)**2
           acc rr1 startl+Kbar endl rr2
           COMPUTE aic(jj,lagnum) = $
               log(rr2(endl)/%nobs) + 2.*%nreg/%nobs
           COMPUTE sic(jj,lagnum) = $
               log(rr2(endl)/%nobs) + log(%nobs)*(1.*%nreg/%nobs) 
        END DO

    END DO 
    }
ELSE 
   {
    DIM aic(vars,Kbar)
    COMPUTE aic = %CONST(0.) 
    DIM sic(vars,Kbar)
    COMPUTE sic = %CONST(0.) 
    DO lagnum = 1,Kbar

        DIM RESIDVEC(vars)
        SYSTEM 1 to vars
        LAGS 1 to lagnum  
        VARIABLES depvar
        DETERMINISTIC CONSTANT
        END(SYSTEM)
        ESTIMATE(noftest,noprint) startl+Kbar endl RESIDVEC(1)
     
        DO jj = 1,vars
           set rr1 startl+Kbar endl = RESIDVEC(jj)**2
           acc rr1 startl+Kbar endl rr2 
           COMPUTE aic(jj,lagnum) = $
                log(rr2(endl)/%nobs) + 2.*%nreg/%nobs 
           COMPUTE sic(jj,lagnum) = $
                log(rr2(endl)/%nobs) + log(%nobs)*(1.*%nreg/%nobs) 
        END DO

    END DO 
    }
   
SET AICVEC 1 Kbar = 0.; SET SICVEC 1 Kbar = 0.
DIM AICMIN(vars) SICMIN(vars) AICMINENT(vars) SICMINENT(vars)
 
DO jj = 1,vars
   DO ii = 1,Kbar
      COMPUTE aicvec(ii) = aic(jj,ii)
      COMPUTE sicvec(ii) = sic(jj,ii)
   END DO
   EXTREMUM(noprint) aicvec 1 Kbar
   COMPUTE aicminent(jj) = %minent
   COMPUTE aicmin(jj) = %minimum 
   EXTREMUM(noprint) sicvec  1 Kbar
   COMPUTE sicminent(jj) = %minent
   COMPUTE sicmin(jj) = %minimum
END DO  

 IF MSC.eq.1 {
    COMPUTE critmax = FIX(AICMINENT(1))
    DO jj = 1,vars
       IF critmax.lt.(FIX(AICMINENT(jj)))
          COMPUTE critmax = FIX(AICMINENT(jj))
    END DO
    }
 ELSE {
    COMPUTE critmax = FIX(SICMINENT(1))
    DO jj = 1,vars
       IF critmax.lt.(FIX(SICMINENT(jj)))
          COMPUTE critmax = FIX(SICMINENT(jj))
    END DO
    }

IF PRINT.eq.1 {
 DISP '  The maximum lag length considered is' Kbar
 DISP  
 WRITE(noskip,format='(7X8A9)') VLABEL 
 DISP ' '
 DISP 'SIC'
 WRITE(noskip,format='(8F9.0)') SICMINENT 
 DISP 'AIC'
 WRITE(noskip,format='(8F9.0)') AICMINENT  
 DISP 
 IF MSC.eq.1
    DISP '  The AIC is being used as the model selection criterion'
 ELSE
    DISP '  The SIC is being used as the model selection criterion'
 DISP
 }
} ****************************

ELSE {
 IF PRINT.eq.1 {
   DISP '   Using' Kbar 'lags in each equation in the VAR'
   DISP
  }
   DIM AICMIN(vars) SICMIN(vars) AICMINENT(vars) SICMINENT(vars)
   COMPUTE AICMINENT = %const(kbar)
   COMPUTE SICMINENT = %const(kbar) 
   COMPUTE critmax = kbar
  }
 
*****************************************************************
 
 DIM laglevsx2(vars) RESIDS(vars) COEFFVEC(vars)
 DO jj = 1,vars
    DIM COEFFVEC(jj)(FIX(AICMINENT(jj)))
    DIM LAGLEVSX2(jj)(0)
    DO i = 1,vars
        IF (MSC.eq.1)
           ENTER(varying) laglevsx2(jj) 
           # laglevsx2(jj) varlist(i){1 to FIX(AICMINENT(jj))}
        ELSE
           ENTER(varying) laglevsx2(jj) 
           # laglevsx2(jj) varlist(i){1 to FIX(SICMINENT(jj))}
    END DO
    IF (DETER.eq.1) {
      LINREG(noprint) depvar(jj) startl+Kbar endl RESIDS(jj)
      #laglevsx2(jj)  
      COMPUTE COEFFVEC(jj) = %beta
      }
    ELSE {
      LINREG(noprint) depvar(jj) startl+Kbar endl RESIDS(jj)
      #laglevsx2(jj) constant
      COMPUTE COEFFVEC(jj) = %beta
      }
 END DO
  
 DIM U(vars,endl-startl-Kbar+1)  Y(vars,endl-startl-Kbar+1)
 DO jj = 1,vars
    DO ii = startl+Kbar,endl
       COMPUTE U(jj,ii-(startl+Kbar)+1) = RESIDS(jj)(ii)
       COMPUTE Y(jj,ii-(startl+Kbar)+1) = DEPVAR(jj)(ii)
    END DO
 END DO
 
 COMPUTE nin =  1./(endl-(startl+Kbar)+1)
 COMPUTE VCOVMAT = nin*(U*tr(U))
  

DIM COEFFMAT(critmax)
DO jj = 1,critmax
   DIM COEFFMAT(jj)(vars,vars)
   COMPUTE COEFFMAT(jj) = %CONST(0.)
END DO


IF MSC.eq.1 {
  DO ii = 1,vars                   ; ********  EQ #
     DO kk = 1,FIX(AICMINENT(ii))  ; ********  # LAGS
       DO jj = 1,vars             ; ********  # of VARIABLE in iith EQ
          COMPUTE COEFFMAT(kk)(ii,jj) = $
             coeffvec(ii)(kk+(jj-1)*FIX(AICMINENT(ii)) )
       END DO
    END DO 
  END DO   
  }
ELSE {
  DO ii = 1,vars                   ; ********  EQ #
     DO kk = 1,FIX(SICMINENT(ii))  ; ********  # LAGS
       DO jj = 1,vars              ; ********  # of VARIABLE in iith EQ
          COMPUTE COEFFMAT(kk)(ii,jj) = $
             coeffvec(ii)(kk+(jj-1)*FIX(SICMINENT(ii)) )
       END DO
    END DO 
  END DO   
  }

COMPUTE SUMMAT = %IDENTITY(vars)
COMPUTE SUMMAT2 = %IDENTITY(vars)
   
DO jj = 1,critmax
   COMPUTE SUMMAT = SUMMAT - COEFFMAT(jj) 
END DO
COMPUTE SDM = INV(SUMMAT)*VCOVMAT*TR(INV(SUMMAT))  

IF PRINT.eq.1 {
 DISP ' '
 DISP 'The Estimated Spectral Density Matrix at the Zero Frequency'
 DISP ' '
 WRITE(noskip,format='(7X8A14)') VLABEL
 WRITE(noskip,format='(8g14.5)') SDM 
} 

END 
