Automatic Differentiation

This example shows how GAUSS uses the symbolic procs created by Symbolic Tools to undertake automatic differentiation.  The GAUSS code for the likelihood function is given in the proc fct; the equivalent code for creating the gradients and Hessian is given in the string txt.  The code is effectively identical - the former deals with matrices, while the latter deals with symbols.  Almost all GAUSS operations and functions are permitted when dealing with symbols.  gradp and hessp are overloaded to permit symbol operation.  NLP, like most optimization techniques in GAUSS, permits the user to specify procs for evaluating the gradient or Hessian - in this case by defining the globals  _nlp_Gradproc  and   _nlp_Hessproc.  In this example, analytic estimation of the gradient and Hessian is 8.5 times faster than using forward differencing.

/* Symbolic example: nlp.e

       proc boxcox(y,indx,sigma, lamda);
       y is the censored variable
       indx is the vector of the index
       sigma is the parameter for the standard error of the residual
       lamda is the parameter for the boxcox process


   This is an estimation of a BoxCox process, with 5000 observations and
   10 parameters.  In this example, artificial data is generated
   and the parameters are estimated using the NLP optimization
   package (part of FANPAC).  The Newton-Raphson algorithm is used,
   with a step length of unity.

   NLP minimizes the function, so to maximize the llf, minimize -llf

   Pindyck and Rubinfeld, 4th ed, p 278
*/


library fanpac, symbolic;  
nlpset;
call symstate(reset);
call symdebug(off);

  rndseed 12345;
  grad_mode = 1;                    // 0 - numeric; 1 - analytic
  display_mode = 0;                 // 0 - off; 1 - on
  nobs = 5000;                      // number of observations
  
                                    // create the data

  xmat = ones(nobs,1) ~rndu(nobs,7);
 
  x1= xmat[.,1];   x2= xmat[.,2]; 
  x3= xmat[.,3];   x4= xmat[.,4];
  x5= xmat[.,5];   x6= xmat[.,6];
  x7= xmat[.,7];   x8= xmat[.,8];
  let bhat = 1 2 3 4 5 6 7 8;
  lmda = 2;
  y = (lmda*xmat*bhat + 1)^(1/lmda) + .01*rndn(nobs,1);

 // Likelihood functionfor Box-Cox process

 proc fct(avec);
   local llf1,  indx, bvec, sig,lam,cnst, llf;
   bvec = avec[1:8];
   sig = avec[9];
   lam = avec[10];
   indx = xmat*bvec;
   cnst = ln(2*pi);
   llf = -.5*( cnst + ln(sig^2) + (((y^lam -1)/lam) - indx)^2/sig^2 ) 
                  +  (lam-1)*ln(y);
   retp(-sumc(llf));
 endp;


 // Gauss symbolic representation of Box-Cox as string

txt = 
    "
     slist = symset(avec,b,10);
     xlist = symlist(8,x);
     sig = b9;
     lam = b10;
     indx = xlist * slist[1:8];
     cnst = ln(2*pi);
     llf = -sumc(-.5*( cnst + ln(sig^2) + (((y^lam -1)/lam) - indx)^2/sig^2 ) 
                     +  (lam-1)*ln(y));
     llfg = gradp(llf,slist);
     llfh = hessp(llf,slist);
";





proc fgrad; endp;
proc fhess; endp;
cls;
if grad_mode == 1 ;
    print /flush "Estimating BoxCox model using NLP, analytic grad and Hess";
    call symproc("fgrad","avec","llfg",txt);
    call symproc("fhess","avec","llfh",txt);
    _nlp_Gradproc = &fgrad;
    _nlp_Hessproc = &fhess;
else;
    print /flush "Estimating BoxCox model using NLP, numeric grad and Hess";
endif;

if display_mode == 1;
  _nlp_Diagnostic = 1;          // print out iterations
endif;
_nlp_GradMethod = 1;            // forward difference
_nlp_MaxIters = 1000;           // max iterations
_nlp_Switch = {.};              // not switch
_nlp_LineSearch = 1;            // step length = 1
_nlp_Algorithm = 3;             // Newton-Raphson

avec0 = 1*ones(10,1);           // starting values 
d1 = date;
{x, f, g, retcode} = nlp(&fct,avec0);
time1 = ethsec(d1,date)/100;
"\n ========================================
 Estimation Results for Box-Cox process
 ========================================";
"\n  LLF value:   " -f;
"\n  Iterations:  " _nlp_niter;
"\n  Coefficient vector " x;
"\n  Elapsed Time: " time1;

The output from this example is :

    grad_mode = 0;           (numeric gradients and hessian)

========================================
Estimation Results for Box-Cox process
========================================

   grad_mode = 0; (numeric)        	 grad_mode = 1; (analytic)    

LLF value: 15766.373 			LLF value: 15766.374 

Iterations: 42.000000 			Iterations: 40.000000 

Coefficient vector 			Coefficient vector 
1.0038481 				1.0033551 	
1.9995901 				1.9998821 
2.9934115 				2.9938492
3.9981160 				3.9987023
4.9974498 				4.9981792 
5.9960651 				5.9969432
6.9976282 				6.9986546
8.0011099 				8.0022799 
0.062763255 				0.062768393 
1.9998455 				1.9999265

Elapsed Time :16.922000 		Elapsed Time: 2.0630000 

Notes:

1. MAXLIK and CML require two input arguments for the gradient and Hessian functions - a Kx1 vector of parameter values, and an NxM matrix of data as input, and a NxK matrix of gradients as an output. This is supported by Symbolic Tools.

2. The gradient and Hessians can either be created on the fly before the estimation:

 call symproc("grdcode","bvec,dta","llfg",txt);
 __gdproc = &grdcode;
 {x,f,g,cov} = maxlik ... 

Alternatively one could run symproc and either retrieve the code from the string, or save it to file. The proc could then be incorporated into the optimization code in the normal way. 

str = symproc("grdcode","bvec","llfg",txt); 
print str; call keyw; // display in GAUSS 
fname = "c:\\temp\\grdcode.txt";
call putf(fname,str,1,strlen(str),0,0); // write code to file