PREVContentsINDEX

Chapter 3 Macro Economics Modeling with AD Model Builder

An active field in macroeconomic modeling is the area of ``regime switching''. This is discussed in greater generality in Hamilton (1994, chapter 22) Hamilton, James D. 1994. Time Series Analysis, Princeton, N.J.: Princeton University Press.. The code for the following example is based on the domain switching model taken from Hamilton (1989) A new approach to the economic analysis of nonstationary time series and the business cycle, Econometrica, 57(2):357-384.. This example is not ideal for exploiting AD Model Builder's greatest advantage, the ability to estimate parameters in models with a large number of independent variables. However it does illustrate the efficacy of the use of higher (up to seven dimensional) arrays in AD Model Builder.

Section 1 Anaylsis of economic data from Hamilton's 1989 paper

For this model The observed quantities are the where

and the state variables satisfy the fourth order autoregressive relationship

where the are independent, normally distributed random variables with mean and standard deviation . These equations correspond to Hamilton's equations 4.3. The state variable is the realized value of a Markov process, , whose evolution is described below. This coefficient takes on the value when the system is in state . In the current example there are two states so that takes on one of the two values or . We can solve for the values of conditioned on the unknown value of the state at time . Let be defined by

Let be a quintuplet of state values for the states at time Define , the realized values of the random variables by

Notice that we due to the lags we can only begin to calculate values for the in time period . It is assumed that the states transitions are given by a Markov process with transition matrix Hamilton seems to index his matrices with the column index first in some cases. We use the row index first. Thus Hamilton's may correspond to our . If we are in state at time the probability of being in state at time is .

If we consider the quintuple of the last 5 states to be the states of a new markov process then we can define the transition matrix for this process by

and

If is the probability of being in state at period the probability of being in state at time period is given by

In particular if

is the probability of being in the state before observing and is the probability of being in the state after observing then

Let be the conditional probability (or probability density) for given . Then, ignoring a constant term which is irrelevant for the calculations,

Define by

Then can be calculated from the relationship

The log-likelihood function for the parameters can be calculated from the . It is equal to

The sums needed for the calculations in can be saved from the calculations for ).

Section 2 The code for Hamilton's fourth order autoregressive model

The complete AD Model Builder template (TPL) code is in the file ham4.tpl The C++ (CPP) code produced from this is in the file ham4.cpp Here is the TPL code split up with comments.

DATA_SECTION
  init_number a1init   // read in the initial value of a1 with the data
  init_int nperiods1   // the number of observations
  int nperiods  // nperiods-1 after differencing
 !! nperiods=nperiods1-1;
  init_vector yraw(1,nperiods1)  //read in the observations
  vector y(1,nperiods)   // the differenced observations
 !! y=100.*(--log(yraw(2,nperiods1)) - log(yraw(1,nperiods))); 
  int order 
  int op1  
 !! order=4; //order of the autoregressive process
 !! op1=order+1;
  int nstates  // the number of states (expansion and contraction)
 !! nstates=2;
The DATA_SECTION contains constant quantities or ``data''. This is in contrast to quantities which depend on parameters being estimated which go into the PARAMETER_SECTION. All quantities in the PARAMETER_SECTION with the init_ prefix are initial data which must be read in from somewhere. By default they are read in from the file ROOT.dat (DAT file) where ROOT is the root part of the name of the program being run (in this case ham4.exe), so ham4.dat.

The first quantity is a number, a1init which will be used for initializing the value of a1 in the program. This is a simple way to try different initial values for a1 simply by modifying the input data file. Such procedures are often valuable to ensure that the correct global value of the objective function has been found. The second quantity nperiods1 is the number of data points in the file. Notice that as soon as a quantity has been defined it is available to use for defining other quantities. The quantity nperiod does not have an init_ before it so it will not be read in an must be calulated in terms of other quantities at some point. Since we want it now it is calculated immediately.


 !! nperiods=nperiods1-1;
The !! are used to insert any valid C++ code into the DATA_SECTION or PARAMETER_SECTION (see LOCAL_CALCS). This code will be executed verbatim (after the !! have been stripped off of course) at the appropriate time. The init_vector yraw is defined and give a size with indices going from 1 to nperiods1. The nperiods1 data points will be read into yraw from the DAT file. The data are immediately transformed and the resulting nperiods data point are put into y.

PARAMETER_SECTION
  init_vector f(1,order,1)  // coefficients for the atuoregressive
                            // process
  init_bounded_matrix Pcoff(0,nstates-1,0,nstates-1,.01,.99,2)  
        // determines the transition matrix for the markov process
  init_number a0(5)  // equation 4.3 in Hamilton (1989)
  init_bounded_number a1(0.0,10.0,4);  
 !! if (a0==0.0) a1=a1init;  // set initial value for a1 as specified
                     // in the top of the file nham4.dat
  init_bounded_number smult(0.01,1,3)  // used in computing sigma
  matrix z(1,nperiods,0,1)  // computed via equation 4.3 in 
                          // Hamilton (1989)
  matrix qbefore(op1,nperiods,0,1);  // prob. of being in state before
  matrix qafter(op1,nperiods,0,1); // and after observing y(t)
  number sigma // variance of epsilon(t) in equation 4.3
  number var  // square of sigma
  sdreport_matrix P(0,nstates-1,0,nstates-1);
  number ff1;
  vector qb1(0,1); 
  matrix qb2(0,1,0,1); 
  3darray qb3(0,1,0,1,0,1);
  4darray qb4(0,1,0,1,0,1,0,1);
  6darray qb(op1,nperiods,0,1,0,1,0,1,0,1,0,1); 
  6darray qa(op1,nperiods,0,1,0,1,0,1,0,1,0,1);
  6darray eps(op1,nperiods,0,1,0,1,0,1,0,1,0,1);
  6darray eps2(op1,nperiods,0,1,0,1,0,1,0,1,0,1);
  6darray prob(op1,nperiods,0,1,0,1,0,1,0,1,0,1);
  objective_function_value ff;
The PARAMETER_SECTION describes the parameters of the model, that is, the quantities to be estimated. Quantities which which have the prefix init_ are akin to the independent variables from which the log-likelihood function (or more generally any objective function) can be calculated. Other objects are dependent variables which must be calculated from the independent variables. The default behaviour of AD Model Builder is to read in initial parameter values for the parameters from a PAR file if it finds one. Otherwise they are given default values consistent with their type. The quantity f is a vector of four coefficents for the autoregressive process. Pcoff is a matrix which is used to parameterize The transition matrix P for the Markov process. Its values are restricted to lie between and . smult is a number used to parameterize sigma and var (which is the variance) as a multiple of the mean squared residuals. This reparameterization undimensionalizes the calculation and is a good technique to employ for nonlinear modeling in general. The transition matrix P is defined to be of type sdreport_matrix so that the standard deviation estimates for its members will be included in the standard deviation report contained in the STD file. To date AD Model Builder suports up to seven dimensional arrays. For historical reasons one and two dimensional arrays are referred to as vector and matrix. This becomes a bit difficult for higer dimensional arrays so they are simply referred to as 3darray, 4darray,... , 7darray.

PROCEDURE_SECTION
  P=Pcoff;
  dvar_vector ssum=colsum(P);  // form a vector whose elements are the
                           // sums of the columns of P
  ff+=norm2(log(ssum)); // this is a penalty so that the hessian will
                        // not be singular and the coefficients of P 
                        // will be well defined
  // normalize the transition matrix P so its columns sum to 1
  int j;
  for (j=0;j<=nstates-1;j++)
  {
    for (int i=0;i<=nstates-1;i++)
    {
      P(i,j)/=ssum(j);
    }
  }  

// get z into a useful format dvar_matrix ztrans(0,1,1,nperiods); ztrans(0)=y-a0; ztrans(1)=y-a0-a1; z=trans(ztrans); int t,i,k,l,m,n;

qb1(0)=(1.0-P(1,1))/(2.0-P(0,0)-P(1,1)); // unconditional distribution qb1(1)=1.0-qb1(0);

// for periods 2 through 4 there are no observations to condition // the state distributions on so we use the unconditional distributions // obtained by multiplying by the transition matrix P. for (i=0;i<=1;i++) { for (j=0;j<=1;j++) qb2(i,j)=P(i,j)*qb1(j); }

for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) qb3(i,j,k)=P(i,j)*qb2(j,k); } }

for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) qb4(i,j,k,l)=P(i,j)*qb3(j,k,l); } } }

// qb(5) is the probabilibility of being in one of 32 // states (32=2x2x2x2x2) in periods 5,4,3,2,1 before observing // y(5) for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) qb(op1,i,j,k,l,m)=P(i,j)*qb4(j,k,l,m); } } } } // now calculate the realized values for epsilon for all // possible combinations of states for (t=op1;t<=nperiods;t++) { for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) { eps(t,i,j,k,l,m)=z(t,i)-phi(z(t-1,j), z(t-2,k),z(t-3,l),z(t-4,m),f); eps2(t,i,j,k,l,m)=square(eps(t,i,j,k,l,m)); } } } } } } // calculate the mean squared "residuals" for use in // "undimensionalized" parameterization of sigma dvariable eps2sum=sum(eps2); var=smult*eps2sum/(32.0*(nperiods-4)); sigma=sqrt(var);

for (t=op1;t<=nperiods;t++) { for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) prob(t,i,j,k)=exp(eps2(t,i,j,k)/(-2.*var))/sigma; } } }

for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) qa(op1,i,j,k,l,m)= qb(op1,i,j,k,l,m)* prob(op1,i,j,k,l,m); } } } } ff1=0.0; qbefore(op1,0)=sum(qb(op1,0)); qbefore(op1,1)=sum(qb(op1,1)); qafter(op1,0)=sum(qa(op1,0)); qafter(op1,1)=sum(qa(op1,1)); dvariable sumqa=sum(qafter(op1)); qa(op1)/=sumqa; qafter(op1,0)/=sumqa; qafter(op1,1)/=sumqa; ff1-=log(1.e-50+sumqa); for (t=op1+1;t<=nperiods;t++) { // notice that the t loop includes 2 for (i=0;i<=1;i++) { // i,j,k,l,m blocks for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) { qb(t,i,j,k,l,m).initialize(); // here is where having 6 dimensional arrays makes the // formula for moving the state distributions form period // t-1 to period t easy to program and understand. // Throw away n and accumulate its two values into next // time period after multiplying by transition matrix P for (n=0;n<=1;n++) qb(t,i,j,k,l,m)+=P(i,j)*qa(t-1,j,k,l,m,n); } } } } } for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) qa(t,i,j,k,l,m)=qb(t,i,j,k,l,m)* prob(t,i,j,k,l,m); } } } } qbefore(t,0)=sum(qb(t,0)); qbefore(t,1)=sum(qb(t,1)); qafter(t,0)=sum(qa(t,0)); qafter(t,1)=sum(qa(t,1)); dvariable sumqa=sum(qafter(t)); qa(t)/=sumqa; qafter(t,0)/=sumqa; qafter(t,1)/=sumqa; ff1-=log(1.e-50+sumqa); // add small constant to avoid log(0) } ff+=ff1; //ff1 is minus the log-likelihood ff+=.1*norm2(f); // add small penalty to stabilize estimation

The PROCEDURE SECTION is where the calculation of the objective function are carried out. First the transition matrix P is calculated from the Pcoff. The function colsum forms a vector whose elements are the column sums of the matrix. This is used to normalize P so that its columns sum to . A penalty is added to the objective function for the colum sums so that the hessian matrix with respect to the independent variables will not be singular. This does not affect the ``statistical'' properties of the parameters of interest. The matrix z is calculated using a transformed matrix because AD Model Builder deals with vector rows better than columns. The probability distribution for the states in period 1, qb1 is set equal to the uncondtional distribution for a Markov process in terms of its transition matrix, P, as discussed in Hamilton (1994). The transition matrix is used to compute the probability distribution of the states in periods , , , and finally . For the last quintuplet this is the probability distribution before observing y(5). The quantities eps in the code correspond to the possible realized values of the random variable . The quantities qa and qb correspond to and in the documentation. The sum function is defined for arrays of any dimension and simply forms the sum of all the components. In AD Model Builder if xx is an n dimensional array then x(i) is an n-1 dimensional array. So the statement

    qbefore(t,0)=sum(qb(t,0));
takes the sum of the probabilities for the sixteen quintuples of states at time period t through t-4 for which the state at time period t is . These are used in the REPORT_SECTION to write out a report of the estimated state probabilities at time period t before and after observing y(t).

REPORT_SECTION
  dvar_matrix out(1,2,op1,nperiods);
  dvar_matrix out1(1,1,op1,nperiods);
  out(1)=trans(qbefore)(1);
  out(2)=trans(qafter)(1);
  {
    ofstream ofs("qbefore.rep");
    out1(1)=trans(qbefore)(0);
    ofs << trans(out1)<< endl;
  }  
  {
    ofstream ofs("qafter.rep");
    out1(1)=trans(qafter)(0);
    ofs << trans(out1) << endl;
  }  
  report << "#qbefore    qafter" <<  endl;
  report << setfixed << setprecision(3) << setw(7) << trans(out) << endl;
The REPORT_SECTION is used to report any result in a manner not already carried out by the models default behaviour. The probabilites of being in state before and after observing y(t) are printed into the files qbefore.rep and qafter.rep. These vectors were stored in files so that they could be easily imported into graphing programs. The results are very similar to figure 1 in Hamilton (1989) as one might hope.

RUNTIME_SECTION
  maximum_function_evaluations 20000 
  convergence_criteria 1.e-6 
maximum_function_evaluations The maximum_functio_evaluations 20000 will simply let the program run a long time by setting the maximum number of function evaluations in the function minimizer equal to 20,000 (nowhere near this many are actually needed.) The convergence_criteria 1.e-6 was needed becuase the default value of 1.e-4 caused the program to exit from the miniization before convergence had been achieved.

TOP_OF_MAIN_SECTION
  arrmblsize=500000;
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(200000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(2100000);
The TOP_OF_MAIN_SECTION is for including code which will be included at the top of the main() function in the C++ program. Any desired legal code may be included. There are a number of common statements which are used to control aspects of AD Model Builder's performance. The statement arrmblsize=500000; reserves 500,000 bytes of memory for variable objects. If it is not large enough a message will be printed out at run time. See the index for references to more discussion of this matter. The statements gradient_structure::set_GRADSTACK_BUFFER_SIZE(200000); and gradient_structure::set_CMPDIF_BUFFER_SIZE(2100000); set the amount of memory that AD Model Builder reserves for variable objects. Setting these is a matter of tuning for optimum performance. If you have a lot of memory available making them larger may improve performance. However models will run without including these statements as long as there is enough memory for AD Model Builder's temporary files.

GLOBALS_SECTION  
  #include 

dvariable phi(const dvariable& a1,const dvariable& a2,const dvariable& a3, const dvariable& a4,const dvar_vector& f) { return a1*f(1)+a2*f(2)+a3*f(3)+a4*f(4); }

The GLOBALS_SECTION is used to include statements at the top of the file containing the CPP program. This is generally where global declarations are made in C++, hence its name. However it may be used for any legal statments such as including header files for the users data structures etc. In this case it has been used to define the function phi which is used to simplify the code for the model's calculations. The header file admodel.hpp is included to define the AUTODIF structures used in the definition of the function. This header is automatically included near the top of the file, but this would be too late as GLOBALS_SECTION material is included first.

Section 3 Results of the analysis

The parameter estimates for the initial parameters are written into a file HAM4.PAR. This is an ASCII file wich can be easily read. (The results are also stored in a binary file HAM4.BAR which can be used to restart the model with more accurate parameters estimates.)

# Objective function value = 60.8934
# f:
 0.0139989 -0.0569580 -0.246292 -0.212250
# Pcoff:
 0.754133 0.0955834
 0.245118 0.900333
# a0:
-0.357964
# a1:
1.52138
# smult:
0.281342
The estimates are almost identical to those reported in Hamilton (1989) Our method for parameterizing the intial state probability distribution qb1 is slightly different from Hamilton's which would explain the small discrepancy. The first line reports the value of the log-likelihood function. This value can be used in hypothesis (likelihood-ratio) tests. the file ham5.para for the fifth order autoregressive model fit to the data in Hamilton (1989) is shown below. there is one more parameter in this model. Twice the difference in the log-likelihood functions is . For one extra parameter the significance level is , the improvement in fit is not significant.


# Objective function value = 59.6039
# f:
 -0.0474771 -0.113829 -0.241966 -0.225535 -0.192585
# Pcoff:
 0.779245 0.0951739
 0.219775 0.900719
# a0:
-0.271318
# a1:
1.46301
# smult:
0.259541

The plot of qa and qb demonstrates the extra information about the probability distribution of the current state contained in in the current value of y(t).

The standard deviation and correlation report for the model are in the file ham4.cor reproduced below.

Section 4 Extending Hamilton's model to a fifth order autoregressive process

Hamilton (1989, page 372) remarks that investigating higher order autoregressive processes might be a fruitful area of research. The form of the model is. The first extension of the model is a fifth order autoregressive process.

and the state variables satisfy the fourth order autoregressive relationship

which extend equations 3.1 and 3.2. The TPL file ham5.tpl for the fifth order autoregressive model is reproduced here. By employing higher dimensional arrays the cionversion of the TPL file from a fourth order autoregressive process to a fifth order one is straightforward. An experienced AD Model Builder user can carry out the modifications in under 1 hour. Places where modifications were made were tagged with the comment //!!5.

DATA_SECTION
  init_number a1init   // read in the initial value of a1 with the data
  init_int nperiods1   // the number of observations
  int nperiods  // nperiods-1 after differencing
 !! nperiods=nperiods1-1;
  init_vector yraw(1,nperiods1)  //read in the observations
  vector y(1,nperiods)   // the differenced observations
 !! y=100.*(--log(yraw(2,nperiods1)) - log(yraw(1,nperiods))); 
  int order 
  int op1  
 !! order=5; // !!5 order of the autoregressive process
 !! op1=order+1;
  int nstates  // the number of states (expansion and contraction)
 !! nstates=2;
PARAMETER_SECTION
  init_vector f(1,order,1)  // coefficients for the atuoregressive
                            // process
  init_bounded_matrix Pcoff(0,nstates-1,0,nstates-1,.01,.99,2)  
        // determines the transition matrix for the markov process
  init_number a0(5)  // equation 4.3 in Hamilton (1989)
  init_bounded_number a1(0.0,10.0,4);  
 !! if (a0==0.0) a1=a1init;  // set initial value for a1 as specified
                     // in the top of the file nham4.dat
  init_bounded_number smult(0.01,1,3)  // used in computing sigma
  matrix z(1,nperiods,0,1)  // computed via equation 4.3 in 
                          // Hamilton (1989)
  matrix qbefore(op1,nperiods,0,1);  // prob. of being in state before
  matrix qafter(op1,nperiods,0,1); // and after observing y(t)
  number sigma // variance of epsilon(t) in equation 4.3
  number var  // square of sigma
  sdreport_matrix P(0,nstates-1,0,nstates-1);
  number ff1;
  vector qb1(0,1); 
  matrix qb2(0,1,0,1); 
  3darray qb3(0,1,0,1,0,1);
  4darray qb4(0,1,0,1,0,1,0,1);
  5darray qb5(0,1,0,1,0,1,0,1,0,1); // !!5
  7darray qb(op1,nperiods,0,1,0,1,0,1,0,1,0,1,0,1); 
  7darray qa(op1,nperiods,0,1,0,1,0,1,0,1,0,1,0,1);
  7darray eps(op1,nperiods,0,1,0,1,0,1,0,1,0,1,0,1);
  7darray eps2(op1,nperiods,0,1,0,1,0,1,0,1,0,1,0,1);
  7darray prob(op1,nperiods,0,1,0,1,0,1,0,1,0,1,0,1);
  objective_function_value ff;
PROCEDURE_SECTION
  P=Pcoff;
  dvar_vector ssum=colsum(P);  // forma a vector whose elements are the
                           // sums of the columns of P
  ff+=norm2(log(ssum)); // this is a penalty so that the hessian will
                        // not be singular and the coefficients of P 
                        // will be well defined
  // normalize the transition matrix P so its columns sum to 1
  int j;
  for (j=0;j<=nstates-1;j++)
  {
    for (int i=0;i<=nstates-1;i++)
    {
      P(i,j)/=ssum(j);
    }
  }  

dvar_matrix ztrans(0,1,1,nperiods); ztrans(0)=y-a0; ztrans(1)=y-a0-a1; z=trans(ztrans); int t,i,k,l,m,n,p;

qb1(0)=(1.0-P(1,1))/(2.0-P(0,0)-P(1,1)); // unconditional distribution qb1(1)=1.0-qb1(0);

// for periods 2 through 4 there are no observations to condition // the state distributions on so we use the unconditional distributions // obtained by multiplying by the transition matrix P. for (i=0;i<=1;i++) { for (j=0;j<=1;j++) qb2(i,j)=P(i,j)*qb1(j); }

for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) qb3(i,j,k)=P(i,j)*qb2(j,k); } }

for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) qb4(i,j,k,l)=P(i,j)*qb3(j,k,l); } } } // !!5 for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) qb5(i,j,k,l,m)=P(i,j)*qb4(j,k,l,m); } } } } // qb(6) is the probabilibility of being in one of 64 // states (64=2x2x2x2x2x2) in periods 5,4,3,2,1 before observing // y(6) for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) { // !!5 for (n=0;n<=1;n++) qb(op1,i,j,k,l,m,n)=P(i,j)*qb5(j,k,l,m,n); } } } } } // now calculate the realized values for epsilon for all // possible combinations of states for (t=op1;t<=nperiods;t++) { for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) { for (n=0;n<=1;n++) { // !!5 eps(t,i,j,k,l,m,n)=z(t,i)-phi(z(t-1,j), z(t-2,k),z(t-3,l),z(t-4,m),z(t-5,n),f); eps2(t,i,j,k,l,m,n)=square(eps(t,i,j,k,l,m,n)); } } } } } } } // calculate the mean squared "residuals" for use in // "undimensionalized" parameterization of sigma dvariable eps2sum=sum(eps2); var=smult*eps2sum/(64.0*(nperiods-4)); //!!5 sigma=sqrt(var);

for (t=op1;t<=nperiods;t++) { for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) //!!5 prob(t,i,j,k,l)=exp(eps2(t,i,j,k,l)/(-2.*var))/sigma; } } } }

for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) { for (n=0;n<=1;n++) qa(op1,i,j,k,l,m,n)= qb(op1,i,j,k,l,m,n)* prob(op1,i,j,k,l,m,n); } } } } } ff1=0.0; qbefore(op1,0)=sum(qb(op1,0)); qbefore(op1,1)=sum(qb(op1,1)); qafter(op1,0)=sum(qa(op1,0)); qafter(op1,1)=sum(qa(op1,1)); dvariable sumqa=sum(qafter(op1)); qa(op1)/=sumqa; qafter(op1,0)/=sumqa; qafter(op1,1)/=sumqa; ff1-=log(1.e-50+sumqa); for (t=op1+1;t<=nperiods;t++) { // notice that the t loop includes 2 for (i=0;i<=1;i++) { // i,j,k,l,m blocks for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) { for (n=0;n<=1;n++) { //!!5 qb(t,i,j,k,l,m,n).initialize(); // here is where having 6 dimensional arrays makes the // formula for moving the state distributions form period // t-1 to period t easy to program and understand. // Throw away n and accumulate its two values into next // time period after multiplying by transition matrix P for (p=0;p<=1;p++) qb(t,i,j,k,l,m,n)+=P(i,j)* qa(t-1,j,k,l,m,n,p); } } } } } } for (i=0;i<=1;i++) { for (j=0;j<=1;j++) { for (k=0;k<=1;k++) { for (l=0;l<=1;l++) { for (m=0;m<=1;m++) { // !!5 for (n=0;n<=1;n++) qa(t,i,j,k,l,m,n)=qb(t,i,j,k,l,m,n)* prob(t,i,j,k,l,m,n); } } } } } qbefore(t,0)=sum(qb(t,0)); qbefore(t,1)=sum(qb(t,1)); qafter(t,0)=sum(qa(t,0)); qafter(t,1)=sum(qa(t,1)); dvariable sumqa=sum(qafter(t)); qa(t)/=sumqa; qafter(t,0)/=sumqa; qafter(t,1)/=sumqa; ff1-=log(1.e-50+sumqa); } ff+=ff1; ff+=.1*norm2(f); REPORT_SECTION dvar_matrix out(1,2,op1,nperiods); out(1)=trans(qbefore)(1); out(2)=trans(qafter)(1); { ofstream ofs("qbefore4.tex"); for (int t=5;t<=nperiods;t++) { ofs << (t-4)/100. << " " << qbefore(t,0) << endl; } } { ofstream ofs("qafter4.tex"); for (int t=5;t<=nperiods;t++) { ofs << (t-4)/100. << " " << qafter(t,0) << endl; } } report << "#qbefore qafter" << endl; report << setfixed << setprecision(3) << setw(7) << trans(out) << endl; RUNTIME_SECTION maximum_function_evaluations 20000 convergence_criteria 1.e-6 TOP_OF_MAIN_SECTION arrmblsize=500000; gradient_structure::set_GRADSTACK_BUFFER_SIZE(400000); gradient_structure::set_CMPDIF_BUFFER_SIZE(2100000); gradient_structure::set_MAX_NVAR_OFFSET(500); GLOBALS_SECTION #include // !!5 dvariable phi(const dvariable& a1,const dvariable& a2,const dvariable& a3, const dvariable& a4,const dvariable& a5,const dvar_vector& f) { return a1*f(1)+a2*f(2)+a3*f(3)+a4*f(4)+a5*f(5); }

Chapter 4 Advanced Featrues of AD Model Builder

Section 1 Using other class libraries in AD Model Builder programs

A useful feature of C++ is its open nature. This means that the user can combine several class libraries into one program. In general this simply involves including the necessary header files in the program and then declaring the appropriate class instances in the program. Instances of external classes can be declared in AD Model Builder program in several ways. They can always be declared in the procedure or report section of the program as local objects. It is sometimes desired to include instances of external classes in a more formal way into an AD Model Builder program. This section describes how to include them into the DATA_SECTION or PARAMETER_SECTION. After that they can be referred to as though they were part of the AD Model Builder code (except for the technicalities to be discussed below).

AD Model Builder employs a strategy of late initialization of class members. The reason for this is to allow time for the user too carry out any calculations which may be necessary for determining parameter values etc. which are used in the initializatin of the object. Because of the nature of constructors in C++ this means that every object declared in the DATA_SECTION or the PARAMETER_SECTION must have a default constructor which takes no arguments. The actual allocation of the object is carried out by a class member function named allocate which takes any desired arguments. If the users external class satisfies these requirments (which is unlikely in general) the desired object can be include into the DATA_SECTION or PARAMETER_SECTION by using the a line of code like


 !!USER_TYPE fooclass myobj( .... )
For most classes which do not satisfy this requirement the desired class object can be included with the statment

 !!USER_TYPEA fooclass myobj( .... )
In this case a pointer to the object is included in the appropriate AD Model Builder class. This pointer has the prefix pad_ inserted before the name of the object. The pointer to myobj would have the form pad_myobj. From the users's point of view the difference between these two apporaches is that in the former approach the object would be referred to by its name myobj while in the latter case the object would be referred to as *pad_myobj. This extra complexity can be avoided by inserting the statment

myclass& myobj = *pad_myobj;
into any function in the procedure section. After that the object can simply be referred to by its name.

Section 2 Appendix 1 -- The regression function

The robust_regression function calculates the log-likelihood function for the standard statistical model of independent normally distributed errors with mean 0 and equal variance. The code is written in terms of Autodif objects such as dvariable and dvar_vector. They are described in the Autodif User's Manual.


dvariable regression(const dvector& obs,const dvar_vector& pred)
{
  double nobs=double(size_count(obs));  // get the number of
                                        // observations
  dvariable vhat=norm2(obs-pred);  // sum of squared deviations
  vhat/=nobs;                      //mean of squared deviations
  return (.5*nobs*log(vhat));     //return log-likelihood value
}

Section 3 Appendix 2 -- AD Model Builder types

The effect of a declaration depends on whether it occurs in the DATA_SECTION or in the PARAMETER_SECTION. Objects declared in the DATA_SECTION are constant, that is like data. Objects declared in the PARAMETER_SECTION are variable, that is like the parameters of the model which are to be estimated. Any objects which depend on variable objects must themselves be variables objects, that is they are declared in the PARAMETER_SECTION and not in the DATA_SECTION.

In the DATA_SECTION the prefix init_ indicates that the object is to be read in from the data file. In the PARAMETER_SECTION the prefix indicates that the object is an initial parameter whose value will be used to calculate the value of other (non initial) parameters. In the PARAMETER_SECTION initial parameters will either have their values read in from a parameter file or will be initialized with their default initial values. The actual default values used can be modified in the INITIALIZATION_SECTION. From a mathematical point of view objects declared with the init_ prefix are independent variables which are used to calculate the objective function being minimized.

The prefixes bounded_ and dev_ can only be used in the PARAMETER_SECTION. The prefix bounded_ restricts the numerical values which an object can take on to lie in a specified bounded interval. The prefix dev_ can only be applied to the declaration of vector objects. It has the effect of restricting the sum of the individual components of the vector object to sum to 0.

The prefix sdreport_ can only be used in the PARAMETER_SECTION. An object declared with this prefix will appear in the covariance matrix report. This provides a convenient method for obtaining estimates for the variance of any parameter which may be of interest. Note that the prefixes sdreport_ and init_ can not both be applied to the same object. There is no need to do so since initial parameters are automatically included in the standard deviations report. AD Model Builder also has three and four dimensional arrays. They are declared like


3darray dthree(1,10,2,20,3,10)
4darray df(1,10,2,20,3,10)
init_3darray dd(1,10,2,20,3,10)    // data section only
init_4darray dxx(1,10,2,20,3,10)   // data section only

Section 4 References

Bard, Yonathan. Nonlinear Parameter Estimation. Academic Press. N.Y. 1974

Gelman, Andrew., John B. Carlin, Hal S. Stern, and Donald B. Rubin. Bayesian Data Analysis. Chapman and Hall.

Hilborn, Ray and Carl Walters. Quantitative Fisheries Stock Assessment and Management: Choice, Dynamics, and Uncertainty. 1992.

Section 5 How to order AD Model Builder

AD Model Builder bundled with Autodif is available for a wide variety of compilers on 80386 computers including Borland C++, Zortech (Symantec) C++, Visual C++ (32 bit) and the ``GNU'' C++ compiler DJGPP Other compilers are supported at present on SUN, HP, and MIPS UNIX workstations. It is important that you tell us the exact form of your hardware and version of your compiler. Multi-user and site licenses are available. Contact Otter Research Ltd PO Box 265, Station A Nanaimo, B.C. V9R 5K9 Canada Voice or Fax (250)-756-0956 Email otterisland.net