Article From: Department of Economics
Subject: MINQUAD.SET

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


               MINQUAD.SET (Heaton & Ogaki)

@ MINQUAD.SET@    @ This file is for Version 2 of GAUSS @

@------------------------------------------------------------------------------
 Written by: John Heaton and Masao Ogaki
             Department of Economics, University of Chicago

 Last Revision:  05/26/92
------------------------------------------------------------------------------@
/* This program defines the proc MINQUAD and the proc OFN and GRAD which
MINAQUAD uses.  MINQUAD computes minimum of a specified quadratic form.
We thank George Tauchen for suggesting the algorithm used here.
   This program uses procedures GRAD2.PRC whcih was in MAXPROC2.ARC and
CONVTEST.PRC, RELVALS.PRC in MAXPROC3.ARC on MODULE7 of version 1.49B  of
GAUSS.  The program package contains these procedures.

   For approximation of hessian, DFP update algorithm or outer products of
derivatives is used depending on a value of global variable "hflag"
(See below).

   After each iteration, the keyboard buffer is tried.
If the buffer contains '1 or 2' "hflag" is set accordingly (See below).

   For each Step Size tried, the keyboard buffer is tried.

If the buffer contains 'c' or 'C', the
   search is stopped just as convergence is attained.
                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            SO WHEN YOU THINK CONVERGENCE IS ATTAINED, HIT 'C' ONCE.
                    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   User is supposed to define global variables and procedures as going through
Step 1 through 3 according to his\her specific problems.  User can replace
numerical derivatives of GRAD2 by his\her own analytical derivatives by
changing "gradname" (See below).

Min Q(b)=u(b)'Wu(b).
 b

where b: kp by 1 vector
      u(b): L by 1 vector  (proc ufn defined by the user)
      W: L by L matrix     (global W0)
*/
@ *********************************************************************** @
@============ Step 1:  Define Global Variables ===========================@
@ *********************************************************************** @
#include "grad2.prc";
#include "convtest.prc";
#include "relvals.prc";
clear gufn,ufnv,vof,relx;  @ No change should be made for this line @

output file=mm.out reset; @ Specify the output file @

   @ Specify the following globals @
w0=eye(1);                             @ Weighting Matrix @
bmin=2;                                @ Initial Values of parametes @

gradname=&GRAD2;       @ Specify the name of proc that returns du(b)/db'
                         Set to &GRAD1 or & GRAD2 for numerical derivatives.
                         User can write his/her own analytical derivatives,
                         and specify the proc's name here.@

hflag=2;               @ Specify the inverse hessian method to use.
                         If hflag = 1, the method is DFP update.
                         If hflag = 2, the method used is outer products of
                         derivatives.
                          Our experience shows that outer products method work
                         better when parameters are near convergence.
                         Otherwise DFP method works better in some examples.@
dfpflag=1;       @ if dfpflag /= 1, initial h=dfpflag @
sstol=1e-100;   @ When 'Step Size<sstol the search is stopped
                 just as convergence is attained.@
    @ See MAXMUM.DOC on MODULE9 of Version 1.49B of GAUSS for
     the following globals @
gradtol=1e-5;
btol=1e-5;
typf=1;
typb=1;
@ ********************************************************************* @
@ ======== Step 2: Define the proc that prints out results ============ @
@ ********************************************************************* @
proc (0)=prntrslt(x,s);

  @ The following message will be printed when converged.  Modify as you like@
"

========================= Minimization Results =======================";
datestr(0);
timestr(0);
?
"Data Set is
Step Size: " s "
Value of the objective function: " vof "
Minimiser is"  x';
? "
=======================================================================";
endp;
@ ********************************************************************* @
@------------- Step 3: Define the proc that returns u(b) -------------- @
@ ********************************************************************* @
proc (1)=ufn(b);
retp(b);
endp;


/* --------------------------------------------------------------------------
Usually the user does not have to change the code below.
--------------------------------------------------------------------------- */

@-------------------- The initial value of the parameters for the
                       routine                                                @

proc startval;
retp(bmin);
endp;

@-------------------- Proc that defines the function to be
                          minimized.                                      ----@


proc ofn(b);

local fval;

ufnv=ufn(b);

fval=ufnv'w0*ufnv;

retp(fval);

endp;

@------------------- Proc that defines the derivatives ------@


proc grad(&gradn,b);

local gradn:proc;
local gqfn;

gufn=gradn(&ufn,b,ufnv,0);

gqfn=gufn'*w0*ufnv;

retp(gqfn');

endp;

@---------------- The proc that defines the inverse of the
                          hessian                                         ----@


proc hessiv(coff);

local hess1, norm1, norm2,kp;

kp=cols(gufn);

hess1=gufn'*w0*gufn;

trap 1;

hess1=invpd(hess1);

if scalerr(hess1); hess1=eye(kp)/maxc(abs(vof)|typf); endif;

trap 0;

@norm2=maxc(maxc(abs(hess1)));

if norm2 > 100;

norm1=minc(minc(abs(hess1)));

hess1=eye(kp)*norm1;   endif; @

retp(hess1);

endp;

@------------ The main procedure -------------------------- @
clear crparv;  @ store current parameter value @
proc minquad;

local x0,g,converge,h,iter,d,s,dx,x,vof1,relgrad,ky,g0jch,vjch,dgjch,
      h1;


@----------------------  Initializations -------------------------------------@

     x0 = startval;                                 @initial parameter value  @
     vof= ofn(x0);                                  @initial function value   @
? "
initial function value is " vof;
     g = grad(gradname,x0);                        @initial gradient value @
     converge=0;                                    @used for convergence test@


@Initial hessian calculation                                    @

    If hflag == 1;
      if dfpflag==1;
        h=eye(rows(x0))/maxc(abs(vof|typf));
      else;
        h=dfpflag;
      endif;
    Else;
        h=hessiv(x0);
    Endif;

     iter = 1;

do until converge;

     d = -h*g';                                     @direction for step       @

     @Choose step direction                                                   @

     s = 1;

     Do while s > sstol;

     dx=s*d;
     x = x0 + dx;

     vof1 = ofn(x);

         if vof1 < vof;
              vof=vof1;
              x0=x;
              goto l1;
              endif;

     s = s*0.5;

     Endo;
" !!!!! WARNING !!!!!!
In MINQUAD Attempted Step Size is smaller than 'ssftol'. Step Size=" s;
 ? "Check the output to see if convergence is attained
!!!!!!!!!!!!!!!!!!";
      goto labc;

l1:

@Calculate relative gradients using gauss routine                             @

     relgrad=relvals(x,dx,g,vof,typb,typf,"relx");

@Convergence test using Gauss routine                                         @

     converge=convtest(relx,relgrad,btol,gradtol);

     if converge; goto labc; endif;

@Update the inverse hessian and the gradient                               @

g0jch=g;

g=grad(gradname,x0);

If hflag == 1;
   dgjch=g-g0jch;

   vjch=h*dgjch';

       @DFP update for the inverse hessian                        @

   h=h + dx*dx'/(dx'dgjch') - vjch*vjch'/(dgjch*vjch);

Else;
    h=hessiv(x0);
Endif;

@Printing results of iteration                           @

"------------- Iteration number " iter " ---------------
Step Size " s "
Value of objective function " vof "
Current parameter values: " x' "
Current relative gradients: " relgrad;
If hflag == 1;
"DFP update used for inverse hessian: DFP=1, Outer Product=2";
Else;
"Outer product used for Hessian: DFP=1, Outer Prodct=2";
Endif;
iter=iter+1;
save crparv=x;
save crhess=h;

    ky=key;
@Now check whether user wants to change the inverse hessian routine  @
    If ky == 49 or ky == 33;  @ 1 or ! @
        hflag=1;
     @   h=eye(rows(x0))/maxc(abs(vof|typf));@
    Elseif ky == 50 or ky == 64;  @ 2 or @
        hflag=2;
    Elseif ky == 67 or ky == 99; @ C or c @
        ? "!!!!!!!!
           MINQUAD interupted by the user";
        goto labc;
    Endif;

endo;
labc:
ky=key;ky=key;ky=key;   @ Clear the buffer @
call prntrslt(x,s);
retp(x);
endp;

@--------------------------End of program------------------------------------@

************************************************************************
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