@ GMM.EXP @
@------------------------------------------------------------------------------
Prepared by Lars P. Hansen, John C. Heaton, and Masao Ogaki
    Financial support from the National Science Foundation
         Grant numbers: SES-8512371 and SES-9213930

 Last Revision:  04/07/93
------------------------------------------------------------------------------@
/*
This program has been used and seems to be free of errors. However,
we (Lars P. Hansen, John C. Heaton, and Masao Ogaki) do not assume
responsibility for any remaining errors.  In no event shall we be liable to
for any damages whatsoever arising out of the use of or inability to use
this program.

This program is for the generalized method of moments (GMM) procedure.
See "Generalized Instrumental Variable Estimation of Nonlinear Rational
Expectations Models," by L.P. Hansen and K.J. Singleton (Econometrica
1982, Vol.50, 1269-1286).  This example file is for Hansen and Singleton
model, but uses different data than theirs.

   @@ To run this example file, type
       RUN MINQUAD.SET [F4][F2]
       RUN GMM.SET [F4][F2]
       RUN GMM.EXP [F4][F2]

      The programs MINQUAD.SET and GMM.SET define necessary procedures.
      To use these programs for the user's problem, the user can
      modify this example file for the user's problem.
      Usually the user does not have to modify MINQUAD.SET and GMM.SET
      because all the paramters for these two programs are controlled
      by this example file.

   @@ The user should go through Step 1-4 and modify this example file for
      the user's own problem.

Model:   E[h(t,b0)]=0:  nmr moments restrictions, h(t,b0) is nmr by 1

Algorithm:
   This program calculates and iterates the following:
     ititial W0 given.

  1. g(b)=C*[h(1,b)+,...,+h(T,b)]
     b=argmin{g(t,b)'W0*g(t,b)}

  2. Rzw(j)=(1/T)[h(1,b)h(1,b)'+,..,+h(T,b)h(T-j,b)']
                       j=0,1,..,mas
     W1=Rzw(0)+{Rzw(1)+Rzw(1)'}+,..,+{Rzw(mas+1)+Rzw(mas+1)'}
     W0=inv(W1)
     (W1 can be calculated by other methods depending on the value of
      calwflag.  See below for more explanations.)

  3. go back to 1.

   Here C is a scaling multiplier, and taken to be 'const' or 'const2'
   depending on whether or not W0 is I.

OUTPUT:
        chi: chi-square test
              statistic for the overidentifying restrictioins
    The following values in GAUSS matrix files:
        a. Current parameter values in MINQUAD.SET search: crparv.fmt
        b. Current inverse hessian matrix in MINQUAD.SET search: crhess.fmt
        c. GMM estimates which uses W0=I: bgmi.fmt
        d. GMM estimates of the last GMM iteration: bgm.fmt
        e. W0 calculated from the last GMM iteration results: w0.fmt
        f. Covariance matrix for bgm: varb.fmt
*/
@ ======================== User Definition Area ============================ @
@ ********************************************************************** @
@ ----------- Step 1: Prepare Output File ---------------------------@
@ ********************************************************************** @
output file=gmmq.out reset; @ Specify the name of the output file. @

? "
************************** GMM Results *********************************";

    @ Prepare the following message which will be printed
      at the top of the output file @
?"
GMMQ.OUT";
? "Hansen and Singleton Model (Single return, VWR) (Econometrica '82)";
datestr(0);
timestr(0);
? "Parameter Names  'beta,gamma'";

@ See GMMQ.EXP for an expamle with multiple returns.@

@ *********************************************************************** @
@------------ Step 2:  Define Global Variables  ------------------------- @
@ *********************************************************************** @

@----------------------------@
   @ Section A. User must specify variables in this section. @
@----------------------------@

tend=334;   @ # of observations=T; scalar @

bgm=1|1;     @ bgm=b(1)|..|b(kgm); kgm by 1 vector
                      Initial values of the coefficients @

nmr=3;           @ scalar: the number of moment restrictions@

mas=0;           @ The order of MA of the disturbance @

@---------------------------------@
     @ Section B.  User can leave the variables in this section as they are,
                   and go to STEP 3.

                   Only advanced users should modify this section.@

@---------------------------------@

gradname=&GRAD2; @ Specify the name of proc that calculates
                    the gradient dgT(b)/db. @

const=1/sqrt(tend);    @ Scaling Multiplier when W0=eye(L) @
const2=1/sqrt(tend);   @ Scaing Multiplier when W0 is not eye(L) @
                       @ These scaling multipliers should be set so that
                         the value of function (vof) in nonlinear search of
                         MINQUAD.SET is near 1.  If vof is too close to 0,
                         the search will not work properly.

                         const2=1/sqrt(tend) will generally be good. @
w0flag=0;    @ scalar;
                If w0flag=0, W0=I is used as the initial weighing matrix W0.

                If w0flag=1, initial bgm is used to calculate initial W0.

                If w0flag=2, W0 in the memory is used as initial W0.

                If w0flag=3, W0 and bgm in the memory are used to give
                       the first GMM result.                       @

maxitegm=5; @ Sets maximum # of iteration over weighting matrix, W0.  Set
                 w0flag=0 and maxitegm=2 to execute usual 2-Stage GMM. @

zero=1E-2;    @ Iteration over W0 continues until the maximum
                 difference of the current and the previous W0 in
                 absolute value becomes less than 'zero', or
                 the # of iteration exceeds maxitegm. @

calwflag=0;  @ This variable is used to choose the method to calculate
               the distance matrics, W0.

               If calwflag=0, Durbin's method will be used when W0
                   is singular.

               If calwflag=1, the QS kernel estimator (nonprewhitened
        or prewhitened) will be used when W0 is singular.

               If calwflag=3, Durbin's method will be used.

               If calwflag=4, the QS kernel estimator (nonprewhitened
                or prewhitened) will be used.
                    @

ordard=mas+3; @ Order of AR representation for
                 Durbin's method @

   @The following four variables control the QS kernel estimator.@

st=0;   @ A scalar to control the bandwidth parameter for the QS kernel.
         if st=0, then an automatic bandwidth estimator is used.
             if st/=0, then st is used as the bandwidth parameter.@

wav=ones(nmr,1); @nmr by 1 vector; weights given to the
      a-th element of z(t)w(t) for the automatic bandwidth estimator
          (used only if st=0). @

maxd=0;  @  scalar:
             if maxd=0, then nonprewhitened HAC with the QS kernel is used.
         if maxd>0, then the elements of DeltaLS with the absolute
                value greater than maxd is replaced by maxd.  See Andrews
                and Monahan's (1990) footnote 4. @

bst=10e+5;  @ scalar:
     When automatic bandwidth parameter is calculated to be bigger than
     bst, bst is used.  @

    @ See MINQUAD.SET for the following globals @
hflag=2;
dfpflag=1;
sstol=1e-25;
    @ See MAXMUM.DOC on MODULE9 of GAUSS for the following globals @
gradtol=1e-5;
btol=1e-5;
typf=1;
typb=1;

@ *************************************************************** @
@ ----------- Step 3: READING IN DATA --------------------------- @
@ *************************************************************** @

load x[335,3]=gmmq.dat;
c=x[1:tend+1,1];     @ c(t)/c(t-1) @
re=x[1:tend+1,2];    @ vwr @
rf=x[1:tend+1,3];    @ T-Bill rate @
clear x;

@ ******************************************************************** @
/* -------------------------------------------------------------------
          Step 4: Define the proc hu(b) that returns tend by L matrix
                             |[z(1)w(1)]'      |
                             |    |            |
                             |[z(tend)w(tend)]'|

                     where [z(t)w(t)]'=[w1(t)z1(t)',...,wnw(t)znw(t)']
----------------------------------------------------------------------- */
@ ********************************************************************** @

@ The proc must have the name "hu", and should take only one argument. @

@ Note that the proc "hu" will be called many times in GMM procedure.
  If some calculations can be done outside the proc, that will speed up
  the program. For example, here vectors of instrumental variables are
  constructed outside the proc. @

  @ Constructing instrumental variables @

  zp=ones(tend,1)~c[1:tend,.]~re[1:tend,.];

  @ Now define the proc "hu" @

proc hu(b);
 local ps,w1,w2,zw1,zw2;

   ps=b[1].*(c[2:tend+1,.]^(-b[2]));
   w1=ps.*re[2:tend+1,.]-1;
   zw1=zp.*w1;

 retp(zw1);
endp;

@=========== The User does not have to change the code below. ============== @

@--------------------- Print results of estimation---------------------@
proc (0)=prntrslt(x,s);

? "
--------------- Minimization Results ----------------------
Step size: " s "
Value of the objective function: " vof "
Minimiser is"  x' "
-----------------------------------------------------------";
endp;

@ ------------------------------------------------------------------------- @
@ ************************************************************************ @
@ THE PROGRAM STARTS @

chi=gmm(gradname,tend,nmr,mas,rows(bgm),zero,maxitegm,st,wav,maxd,bst);
output off;

@ ------------------------ END OF PROGRAM --------------------------------- @