/* ------------------    West 1994 (Wong, Jun95)  ------------------------ **
**  West (1994): "Another Heteroskedasticity Autocorrelation Consistent    **
**        Covariance Matrix Estimation," SSRI working paper #9434,         **
**        University of Wisconsin-Madison.                                 **
         **
** ----------------------------------------------------------------------- **

      Author: Ka-fu Wong                   email: kfwong@macc.wisc.edu
              Department of Economics,
              University of Wisconsin,
              Madison, WI 53706

    This code is written and submitted for public, non-commerical use.
    You are welcome to use and distribute them as long as you include proper
    attribution to the author.  Although this code has been tested
    extensively, we do not guarantee its performance.  Comments and
    suggestions are welcome.  It is last revised on June 12, 1995.
** ----------------------------------------------------------------------- **
The following procedure computes the covariance matrix nonparametrically
by ARMA regression of the residuals.  See West (1994) for details.  A
separate file EXAMPLE.PGM demonstrates the use of this procedure.

Format: {var,gam}=NEWMA(kkk,e,xx,mpar0)

Note that it uses OPTMUM procedure in gauss to compute the nonlinear
least square estimate of MA weights.  This procedure is present only in
the basic package of gauss 3.0 or older versions.

*
* INPUT:
*     kkk:= 2x1 vecotr = ARorder|MAorder.
*       e:= Tx1 residuals.
*      xx:= TxK matrix (usually the RHS variables of an OLS regression)
*   mpar0:= the initial parameters for the nonlinear estimation of the
*           arma model of the cross products of instruments and residuals.
*
* OUTPUT:
*   var:= covariance matrix.
*   gam:= same as the gamma for Bartlett and QS kernels in Newey and West(1994);
*         indicator for non-positive definiteness of the covariance matrix
*         for Truncated kernel for the specified bandwidth.  (When
*         this happens, the bandwidth is set to 1.)
*
-------------------------------------------------------------------------- */
/*
** Returns objective function to be minimized (see the GAUSS application
** OPTMUM  for details)
*/
library optmum;
clear e,tempseed;
proc (1) = objmin(aaa);
local eee,ra,ra_1,i,n;

   aaa=rev(aaa); ra=rows(aaa); ra_1=ra-1; n=rows(e);
   eee=zeros(ra,1);
   i=1; do until i>n;
      eee=eee|(e[i,.]-aaa'eee[i:i+ra_1,.]);
   i=i+1; endo;

   eee=eee[ra+1:rows(eee),.];

retp(sumc(eee^2));
endp;




/*
** Compute MA weights and MA residuals by MA regression
*/
clear expno,seedux,seedux0;
proc (2) = mawold(resid,maorder,mapar0);
local maweight,maresid,minvalue,mingrad,retcode,MA_1,i,n,ll,vcb,aic,sbc,
      mar,tempy,tempret;

   if rows(mapar0)/=maorder;
      "Experiment" expno ": MAorder does not match mapar0";
      maorder=rows(mapar0);
   endif;

   /* estimate MA weights by Nonlinear Least Square */
   {maweight,minvalue,mingrad,retcode}=optmum(&objmin,mapar0);

   if retcode /= 0;
      "Warning: nonlinear least square has not been completed normally!"
      "return code (check optmum procedure description) = " retcode;
      maweight=0; maresid=0;
   else;
      maweight=rev(maweight);
      maresid=zeros(maorder,1);  ma_1=maorder-1; n=rows(resid);
      i=1; do until i>n;
         maresid=maresid|(e[i,.]-maweight'maresid[i:i+ma_1,.]);
      i=i+1; endo;
      maresid=maresid[maorder+1:rows(maresid),.];
   endif;

retp((maweight)|1,maresid);
endp;





/*
** Compute the *new* covariance matrix using d(t).
** See West(1994) for details.
*/
proc (1) = newcov(xxx,arorder,maorder,maweight,maresid);
local tt,dum,dd,jhatnpw,nnn,zz,zz0;

     maweight=rev(maweight);
     nnn=rows(xxx);

     tt=maorder+1; dd={};
     do until tt>nnn;
        dd=dd|maresid[tt-1,.]*maweight'xxx[(tt-maorder):tt,.];
     tt=tt+1; endo;
        @ dd is (n-MAorder)x1 since ARorder<MAorder [see notes] @

     jhatnpw=dd'dd/rows(dd); @ Jhat with NO prewhitening @

retp(jhatnpw);
endp;




/*
** NEW(MA) Procedure
*/
proc (2) = newma(kkk,e,xx,mapar0);
local gam,jhatnpw,maweight,maresid,arorder,maorder;

   arorder=kkk[1,.]; maorder=kkk[2,.]; gam=0;

   {maweight,maresid}=mawold(e,maorder,mapar0);

     if maresid==0;
        gam=9; jhatnpw=0;
     else;
        jhatnpw=newcov(xx,arorder,maorder,maweight,maresid);
     endif;

retp(jhatnpw,gam);
endp;
