Article From: Department of Economics
Subject: GRAD2.PRC

Posted by: Information Provider
Phone Number: (202) 885-3770
E-Mail Address: econ@american.edu

Post Date: 23 Sep 1994
Expiration Date: 23 Oct 2004


               GRAD2.PRC (Ogaki)

/*
GRAD2 --   Procedure to compute the gradient matrix of a vector-valued
           function that has been defined in a procedure. Double-sided
           gradients are computed.  In  y=f(x) the assumption is that
           x is Kx1 and that y is Nx1.  The gradient matrix g will be
           nxk.

           g=grad2(&f,x,f0,dh);
           INPUTS
           f -- a vector valued function (f:Kx1 -> Nx1) defined in a proc.
                It is acceptable for f(x) to have been defined in terms of
                global arguments in addition to x, and thus f can return
                an Nx1 vector (e.g.:
                      proc f(x); y=exp(x*b); retp(y); endp;).
           x -- Kx1 vector
          f0 -- Nx1 vector containing value of f at x: f0=f(x)
          dh -- step length for forward gradient. If dh == 0,
                then step length will be computed internally.
          OUTPUT
           g -- NxK matrix containing the gradients of f with respect
                to the variable x at x.

           The procedure requires that f0 be passed. This makes the
           call identical to the call for GRAD1, and in addition
           the number of rows in the result are computed from this.

NOTE: GRAD2 will return a ROW for every row that is returned by f. For
instance, if f returns a 1x1 result, then GRAD2 will return a 1xK row
vector.  This violates the GAUSS convention that functions return column
vectors.  However, this allows the same function to be used regardless of
N, where N is the number of rows in the result returned by f.
............................................................................*/
proc 1=grad2(f,x,f0,dh);
local f:proc;
local n,k,grdd,ax,dax,xdh,ee,i,eei;

 n=rows(f0);
 k=rows(x);
 grdd=zeros(n,k);

 /* Computation of stepsize (dh) for gradient if dh == 0 */
if dh==0; ax=abs(x);
 if x /= 0; dax=x./ax; else; dax=1; endif;
 dh=(1e-5)*maxc(ax|(1e-2)*ones(rows(x),1)).*dax;
endif;

 xdh=x+dh; dh=xdh-x;        @ This increases precision slightly @
 ee=eye(k).*dh;

 i=1;
 do until i > k;

    eei=ee[.,i];
    grdd[.,i]=f(x+eei)-f(x-eei);

 i=i+1;
 endo;

  grdd=grdd./(2*dh');

retp(grdd);
endp;

************************************************************************
Department of Economics is entirely responsible for the information
provided above. Please direct any comments to Information Provider at:
          E-Mail Address: econ@american.edu
          Phone Number:   (202) 885-3770