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.

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

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

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); } }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// 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

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-6maximum_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 #includeThe 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.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); }

# 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.281342The 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.

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); }

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.

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 }

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

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.