PREV
ContentsINDEX
Chapter 4 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_function_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 \ttfoot 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 \chapno.1 and \chapno.2.
The TPL file ham5.tpl for the fifth order autoregressive model
is reproduced here. By employing higher dimensional arrays the conversion of
the TPL file from a fourth order autoregressive
process to a fifth order one
is largely formal. 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);
}
NEXT