/*************     Newey-West and West Example (Wong, Jun95) ********/

/*
Example file for West (1994) and Newey and West (1994) code
archived by Ka-fu Wong <KFWONG@macc.wisc.edu>.
This example assumes you have the archived Newey-West 1994
and West 1994 source code by Ka-fu Wong and have saved
it as nw1994.src and west1994.src.
*/

NEW; FORMAT /M1/RD 12,5; OUTWIDTH 250;  tstart=DATE;
outfile="example.out"; @<-- Filename for output file. @
OUTPUT FILE = ^outfile RESET;/*          DATESTR(0); TIMESTR(0);*/
"
 Source Program: example.pgm
 Created: 6  January 1994;  Revised:  13  June, 1995    by Ka-fu Wong
 This is an example file to compute the covariance matrix using
 different procedures discussed in Newey and West (1994) and West (1994).
 TR, QS and BT are discussed in Newey and West (1994).  NEWMA is discussed
 in West(1994).
 QS and BT use prewhitening and SVD while TR and NEWMA do not.
 Truncated kernel uses m=population MAorder.


TR   : m = ? if resulting cov p.d.; 0 otherwise.
QS   : n = fit AR(1) to h series
BT   : n = int[4(T/100)^(2/9)]
NEWMA: new estimator (ARorder=0, MAorder=?) using NLS as  21 May, 1994
       ?= true MAorder";?;


#include "nw1994.src"; /* procedures to compute TR, QS and BT covariance */
#include "west1994.src"; /* procedures to compute NEWMA covariance */
                         /* it calls "library optmum" of gauss 3.0 */

@ -----------    Set the initial values   ---------- @
n=128;          @ # of sample size @
nd=0;   @ # of obs to be discarded @

ww={0,1,1,1,1}; @ give equal weight to all variables except constant term @
npar=5;         @ # of regressors in equation @
usigx=1;        @ std dev of x @
usigu=1;        @ std dev of u @
rhoxx=0.5;
psiuu=0.5;
psixx={}; rhouu={};
sigexx=0.8660254038;
sigeuu=0.8944271910;
_opmiter=100;
seedx=34287; seedx0=49375; seedu=134287; seedu0=149375;



    OUTPUT FILE = ^outfile ON; FORMAT /M1/RD 12,5;
    " x's AR parameter = " rhoxx' " MA parameter = " psixx';
    " u's AR parameter = " rhouu' " MA parameter = " psiuu';
    "With homoskedastic disturbances ";
    "random number seeds of the first experiment(seedx0,seedx,seedu0,seedu):";
     seedx0~seedx~seedu0~seedu;

    /** Generate data **/
    {x,seedx,seedx0}=ARMA1(0,rhoxx,psixx,seedx,seedx0,0,sigexx,usigx,4,
                             n+1,nd);
    x=ONES(ROWS(x),1)~x;
    {u,seedu,seedu0}=ARMA1(0,rhouu,0,seedu,seedu0,0,sigeuu,usigu,1,
                             n+1,nd); /* psiuu set to zero to get u iid */
    y=u;


    /** OLS **/
    {b,e,xxinv}=OLSRESID(x,y,1);
    h=x.*e;  @ cross products of instruments and residuals @

    /* Prewhiten h before computing covariances, using VAR(1)
    ** Note: if prewhitening is not needed, do not use this procedure
    **       and set vstar=h; ahat=0;
    */
     {Ahat,Vstar}=PWhite1(-0.97|0.97,h);


  /*** Compute the covariance matrix ***/
     /* trucated kernel */
       {trcov,trgam}=tru(rows(psiuu),h,0); @ using the population order of ma @
?;?;"truncated kernel";
"cov=" trcov;
"gam=" trgam;

     /* quadratic spectral */
       {qscov,qsgam}=qsbar(vstar,ahat,ww);
?;?;"quadratic spectral kernel";
"cov=" qscov;
"gam=" qsgam;

     /* bartlett kernel */
       {btcov,btgam}=btbar(trunc(4*((n/100)^(2/9))),vstar,ahat,ww);
?;?;"bartlett kernel";
"cov=" btcov;
"gam=" btgam;

     /* newma */
       {macov,magam}=newma(0|rows(psiuu),e,x,psiuu);
                                              @ psiuu is used as the initial
                                                value of parameters in the
                                                nonlinear least squares @
?;?;"newma";
"cov=" macov;
"gam=" magam;

OUTPUT FILE=^outfile ON;
tend=DATE; "time elapsed: " ETSTR(ETHSEC(tstart,tend));
END;

/* special procedures used in this program */
/* --------------------------------------------------------------------------
SOURCE PROGRAM: ARMA1.PRC
Purpose: Generate k indepedent ARMA(p,q) time series.
Format:
   xt=ARMA1(const,arpar,mapar,seedx,seedx0,x0,usige,usigx,kseries,n,nd)

*
* INPUT:
*    const:= constant term in the series.
*    arpar:= px1 AR parameters
*    mapar:= qx1 MA parameters
*    seedx:= scalar random number seed
*   seedx0:= scalar random number seed for the initial observation(s);
*            set to zero if want to impose x0.
*        x0:= initial observation(s).
*     usige:= scalar, standard deviation of innovation e (unconditional)
*     usigx:= scalar, standard deviation of x (unconditional).
*   kseries:= scalar, number of independent series generated
*         n:= scalar, number of observations in each series
*        nd:= number of observations to be discarded
*
* OUTPUT:
*       xt:= nxk, k indepedent series of n observations each
*    seedx:= update random seed
*   seedx0:= update random seed
*
* Remarks:
*   Procedure originally written for West(1994) "new covariance estimate"
*         x(t) = const + ARpar(1)*x(t-1) + ... + ARpar(p)*x(t-p) + u(t)
*  where  u(t) = e(t) + MApar(1)*e(t-1) + ... + MApar(q)*e(t-q)
--------------------------------------------------------------------------- */
PROC (3) = ARMA1(const,arpar,mapar,seedx,seedx0,x0,usige,usigx,kseries,n,nd);
    LOCAL p,rarpar,q,rmapar,xt,t,et,ut;
        p=ROWS(arpar);
        rarpar=REV(arpar)';
        q=ROWS(mapar);
            IF q>0; IF mapar[1,.]==0; mapar={}; q=0; ENDIF; ENDIF;
        rmapar=REV(1|mapar)';

        et=usige*RNDNS(n+nd+q,kseries,seedx);

       /** compute the MA(q) disturbances **/
        IF q==0;
           ut=et;
        ELSE;
           t=1; ut={};
           DO UNTIL t>n+nd;
              ut=ut|(rmapar*et[t:t+q,.]);
           t=t+1;
           ENDO;
        ENDIF;

       /** form xt AR(p) using MA(q) disturbances **/
        IF p==0;
            xt=const+ut;
        ELSE;
           t=1; IF seedx0==0; xt=x0.*ONES(p,kseries);
                ELSE;         xt=usigx*RNDNS(p,kseries,seedx0);
                ENDIF;
           DO UNTIL t>n+nd;
               xt=xt|(const+rarpar*xt[t:t+p-1,.]+ut[t,.]);
           t=t+1;
           ENDO;
        ENDIF;

 RETP(xt[nd+p+1:n+nd+p,.],seedx,seedx0);
 ENDP;




@ ------------------------------------------------------------------------
SOURCE FILE: OLSRESID.PRC
Purpose: returns OLS estimates, residuals and x'x inverse.

Format:  {bhat,resid,xxinv}=OLSRESID(x,y,ret_ctrl)

*
* Input:
*          x:= NxK regressors (RHS variables) matrix.
*          y:= Nx1 dependent variable.
*   ret_ctrl:= 1x1.  If 1, will return residuals; if 0, will
*              return a dummy (scalar).
*
* Output:
*    bhat:= Kx1 OLS parameter estimates.
*   resid:= Nx1 OLS residuals.
*   xxinv:= KxK, x'x inverse.
*
* Remarks:
*   See IVRESID.PRC for estimation with Instrumental Variable method.
*
-------------------------------------------------------------------------- @
PROC (3) = OLSresid(xxx,yyy,ret_ctrl);
  LOCAL qqq,bbb,resid;
    qqq=INVPD(MOMENT(xxx,0));
    bbb=qqq*xxx'yyy;    @ regression coefficient estimates @

    IF ret_ctrl==1;
       resid=yyy-xxx*bbb;    @ vector of residuals from OLS @
    ELSE; resid=0;
    ENDIF;

RETP(bbb,resid,qqq);
ENDP;

/* ----------------------------- End of Program --------------------------- */
/**********************   EXAMPLE PROGRAM OUTPUT  **************

 Source Program: example.pgm
 Created: 6  January 1994;  Revised:  13  June, 1995    by Ka-fu Wong
 This is an example file to compute the covariance matrix using
 different procedures discussed in Newey and West (1994) and West (1994).
 TR, QS and BT are discussed in Newey and West (1994).  NEWMA is discussed
 in West(1994).
 QS and BT use prewhitening and SVD while TR and NEWMA do not.
 Truncated kernel uses m=population MAorder.


TR   : m = ? if resulting cov p.d.; 0 otherwise.
QS   : n = fit AR(1) to h series
BT   : n = int[4(T/100)^(2/9)]
NEWMA: new estimator (ARorder=0, MAorder=?) using NLS as  21 May, 1994
       ?= true MAorder

 x's AR parameter =      0.50000  MA parameter =           {}
 u's AR parameter =           {}  MA parameter =      0.50000
With homoskedastic disturbances
random number seeds of the first experiment(seedx0,seedx,seedu0,seedu):
 49375.00000  34287.00000 149375.00000 134287.00000


truncated kernel
cov=
     0.95438     -0.41686      0.07470     -0.06254     -0.18666
    -0.41686      0.94041      0.18613     -0.46633      0.28144
     0.07470      0.18613      0.60453     -0.03835      0.08275
    -0.06254     -0.46633     -0.03835      1.04736     -0.01107
    -0.18666      0.28144      0.08275     -0.01107      0.61920
gam=     0.00000


quadratic spectral kernel
cov=
     1.03891     -0.54426      0.08008      0.03322     -0.21664
    -0.54426      1.04702      0.17935     -0.47799      0.33660
     0.08008      0.17935      0.63429     -0.03427      0.08307
     0.03322     -0.47799     -0.03427      0.96518     -0.11792
    -0.21664      0.33660      0.08307     -0.11792      0.62421
gam=     0.15072


bartlett kernel
cov=
     0.93152     -0.43235      0.15876      0.09594     -0.16949
    -0.43235      0.93381      0.10756     -0.46191      0.34435
     0.15876      0.10756      0.71687      0.07818      0.05868
     0.09594     -0.46191      0.07818      0.92138     -0.04333
    -0.16949      0.34435      0.05868     -0.04333      0.60480
gam=     0.77725


newma
cov=
  0.95502  -0.34420   0.08216  -0.16216  -0.07301
 -0.34420   0.72317   0.09200  -0.18421   0.09694
  0.08216   0.09200   0.50684  -0.06886   0.08467
 -0.16216  -0.18421  -0.06886   1.25836  -0.11732
 -0.07301   0.09694   0.08467  -0.11732   0.74248
gam=  0.00000
time elapsed: 2.04 seconds

************* END OF EXAMPLE PROGRAM OUTPUT ***************/

