The full text on this page is automatically extracted from the file linked above and may contain errors and inconsistencies.
orKing raper beries
3
Algorithms for Solving Dynamic Models
with Occasionally Binding Constraints
Lawrence J. Christiano and Jonas D.M. Fisher
Working Papers Series
Macroeconomic Issues
Research Department
Federal Reserve Bank of Chicago
April 1994 (WP-94-6)
FEDERAL RESERVE BANK
OF CHICAGO
A lg o rith m s fo r S o lvin g D yn a m ic M odels w ith
O cca sio n a lly B in d in g Constraints*
Lawrence J. Christianot
Jonas D.M. Fisher1
April, 1994
A bstract
We describe several methods for approximating the solution to a model in which in
equality constraints occasionally bind, and we compare their performance. We apply
the methods to a particular model economy which satisfies two criteria: it is similar to
the type of model used in actual research applications, and it is sufficiently simple that
we can compute what we presume is virtually the exact solution. We have two results.
First, all the algorithms are reasonably accurate. Second, on the basis of speed, accu
racy and convenience of implementation, one algorithm dominates the rest. We show
how to implement this algorithm in a general multidimensional setting, and discuss the
likelihood that the results based on our example economy generalize.
JEL Classification: C6, C63, C68
Keywords: Occasionally binding constraints. Parameterized Expectations, Collocation
Send correspondence to: Prof. Jonas D.M. Fisher, Department of Economics, Social
Science Centre, University of Western Ontario, London, ON, Canada, N6A 5C2.
We than k G raham C andler, K enneth Ju d d , A lbert M arcet, David M arshall, and Ellen M cG rattan for
helpful com m ents. C h n stia n o is grateful to the N ational Science Foundation for su p p o rt. T he views expressed
herein are those of the authors and not necessarily those of the Federal Reserve B ank of Chicago or the Federal
Keserve System .
'N o rth w estern University, N B ER, Federal Reserve B anks of Chicago and M inneapolis.
+university of W estern O ntario.
1. I n t r o d u c t i o n
In recent years there has been substantial interest in studying the quantitative properties of
dynamic general equilibrium models. For the most part, exact solutions to these models are
unobtainable and so in practice researchers must work with approximations. An increasing
number of the models being studied have inequality constraints that occasionally bind. The
main examples of this are heterogeneous agent models in which there are various kinds of
constraints on the financial assets available to agents. 1 Other examples include multisector
models with limitations on the intersectoral mobility of factors of production, and models
of inventory investment. 2 An important consideration in selecting algorithms for solving
models like these is the quantity of computer and programmer time required to achieve an
acceptable degree of accuracy. The purpose of this paper is to help shed light on these issues.
We describe six algorithms, and evaluate their accuracy in solving the one-sector infinite
horizon optimal growth model with random productivity disturbances. In this model the
nonnegativity constraint on gross investment is occasionally binding. We chose this model
for two reasons. First, its simplicity makes it feasible for us to solve the model by doing
dynamic programming on a very fine capital grid. Because we take the dynamic programming
solution to be virtually exact, this constitutes an important benchmark for evaluating the
four algorithms considered. Second, the one sector growth model is of independent interest,
since it is a basic building block of the type of general equilibrium models analyzed in the
lite r a tu r e .3
All the methods we consider work with the Euler equation associated with the necursive
e'XaT l en’o ^ 1>TaS an
AlyaSa n and G ertler ( 1991)> den Haan (1993), H eaton and Lucas
J 9 ' 1 uSSe tt (1993) k iy o ta k i and Moore (1993), M arcet and K etterer (1989), M arcet and M arimon
(1992), M arcet and Singleton (1990), Telm er (forthcom ing), and M cC urdy and R icketts (1992).
-F or an exam ple of the form er, see Atkeson and Kehoe (1993), and B oldrin, C hristiano and Fisher (1994).
exam ples of the la tte r include C hristiano and F itzgerald (1991) and K ahn (1992).
For exam ple solving th e heterogeneous agent m odels of A iyagari (1993), A iyagari and G ertler (1991)
ant H uggett (1993) requires repeatedly solving a p a rtia l equilibrium asset accum ulation problem for an
me lvidual agent , for different values of a particu lar m arket price. A solution to the general equilibrium
problem is obtained once a value for the m arket price is found which implies a solution to the partial
equilibrium problem th a t satisfies a certain m arket clearing condition. T he p artial equilibrium m odel solved
m these exam ples is sim ilar to the grow th m odel we work w ith in this paper.
n J ? ' £ *
1
representation of the model. Thus, a solution is viewed as a policy function relating decisions
to a small number of state variables. All but one of the algorithms considered work with a
version of the model in which the nonnegativity constraint is incorporated by the method of
Lagrange multipliers. These include suitably modified versions of the algorithms emphasized
by Bizer and Judd (1989), Coleman (1988), Danthine and Donaldson (1981), Judd (1992a),
and Marcet (1988). The sixth algorithm, an example of the Finite Element Method, works
with a version of the model in which the nonnegativity constraint is incorporated by a penalty
function method. This algorithm has been advocated by McGrattan (1993).4
Our main finding is that, for the model economy studied, one algorithm dominates the
others in terms of speed, accuracy and programmer time. This algorithm approximates
the solution indirectly by parameterizing the conditional expectation in the Euler equation
using an exponentiated polynomial, as in Marcet (1988).5 We show that conventional imple
mentations of parameterized expectations have some shortcomings, and document that our
preferred algorithm dominates on these dimensions.
In our example, there are two principal advantages in parameterizing a conditional ex
pectation.
First, the conditional expectation function is smoother than other functions
characterizing the solution, such as the policy function. In general, it is easier to obtain an
accurate approximation, the smoother the function being approximated. A second advan
tage is that working with parameterized expectations is efficient from the point of view of
programmer time. In the context of methods based on Lagrange multipliers, the require
ment that the Euler equations and Kuhn-Tucker conditions be satisfied implies a convenient
mapping from a parameterized expectation function into policy and multiplier functions.
This obviates the need to separately parameterize the latter. Methods which focus on the
policy function must jointly parameterize the policy and multiplier functions. Doing this
4See the chapter in Ju d d (1992b) on rational expectations m odels for references to earlier analyses of
m odels w ith nonnegativity co n stra in ts. As the m aterial in th a t chapter indicates, several of th e m ethods
used in this p aper actually correspond to approaches taken by G ustafson and o th er agricultural econom ists
decades ago.
5For o th er applications of th e P E A when there are occasionally binding co n strain ts, see den H aan (1993),
M arcet and K etterer (1989), M arcet and M arim on (1992), M arcet and Singleton (1990), and M cC urdy and
R icketts (1992).
9
in a way that the Kuhn-Tucker conditions are satisfied is tricky and adds to programmer
time. For methods that focus on policy functions, an alternative to working with Lagrange
multipliers is to work with a penalty function. However, these methods require searching
for a parameter in the penalty function, which can add substantially to programmer and
computer time. Although we carefully document these statements for our model economy,
we expect them to be true in a broader class of models as well.
The paper is organized as follows. In the following section the model to be solved is de
scribed. This is followed by a review of how the six algorithms can be used to approximate
the unconstrained version of the model in which the nonnegativity constraint on invest
ment is ignored. In the next section we describe the way these algorithms are modified to
accommodate the nonnegativity constraint on investment. Results from implementing the
algorithms for a particular parameterization of the model are discussed in section 5. In the
final section we offer some concluding remarks.
2. T h e M o d e l
We examine a simple version of the stochastic growth model with inelastic labor supply.
At date 0 the representative agent values alternative consumption streams according to
/:°
^
=
o
where
ct
denotes time t consumption and
f3 <E (0,1)
is the agent’s discount
iactor. The aggregate resource constraint is given by
Ct + ^*h -i _ (1 ~ $)kt < f{kt,0t) = exp(0t)k°,
(2.1)
where kt denotes the beginning-of-period-t stock of capital, and S, a <E (0,1). Here, 6 is the
late of depreciation of capital, and
exogenous with respect to
kt
a
is capital’s share in production. We assume
0t
e 0 is
and has a first order Markov structure with the density of 0 t+1
conditional on 6t given by pe>{0t+l | 6t). In our computational experiments, we assume 0t is
i.i.d.
with 0 = {<7, - a } , and that the probability associated with each of <r and -<r is 1/2.
The initial stock of capital, k0, is given.
3
In the irreversible investment version of the model, we require that gross investment be
non-negative,
i . e .:
k t+1
-
(1
-
S)kt >
0.
(2.2)
In the reversible investment version, (2.2) is ignored.
We formulate the planner’s problem in recursive form. In doing so we drop
subscripts
t
and use ' to denote next period’s value of a variable. The planner’s dynamic program is then
given by
W ( k , 9) =
’
v
max
f(k,e)+(i-$)k>k'>o
U ( c ( k , k ', 9 ) )
+ j3 f
Je'ee
IE(A-', 9 ’) p e, ( 9 ' | 9 ) d 9 '.
(2.3)
Equation (2.2) must also be satisfied in the irreversible investment version of the model.
Assumptions we will place on
U(-)
guarantee that (2.1) always binds for this economy. In
(2.3) we have used this fact to replace consumption in U(-) with c(-), the function implicit
from (2.1). Finally, IE (•) is the planner’s value function.
To solve the planner’s problem we introduce a Lagrange multiplier, A, on constraint (2.2).
The solution to the planner’s problem is a set of time invariant functions
that determine
k'
and A, respectively, given values of
k
and
9.
g(k,9)
and
h(k,9)
These functions must satisfy
an Euler equation,
U c { k. g ( k , 9). 9)
—h ( k ,
9 ) — f3
f
m ( k , 9 , 9 ' ; g , h ) p 6.(P ' \ 9 ) d 9 '
= 0,
(2.4)
J 01
and a set of Kuhn-Tucker conditions
g{k, 9 ) - {
1-
8)k >
0, h ( k ,
9) >
0, and
h ( k , 9 ) [ g { k , 9) -
(1 -
6)k] =
0.
(2.5)
Also,
m(k.0,O':g,h) = U c(g{k,0),g{g(k,0),$'),9')\fi;(g(k,8),O') + 1 - 6 j
-h(g(k,e)J')(l-S).
In (2.4)-(2.6),
fk
denotes the derivative of / , while
4
Uc
combines the derivative of
U
with the
function c( k ,
k' , 0 ) .
If the standard deviation of the technology shock, a, is small enough, (2.2) will never
bind and A = 0 for all
6
and
k
in the ergodic set for capital. Methods for approximating the
solution to the planner’s problem are well known for this case. Here we consider the case
where cr is large enough so that A > 0 with nonzero probability.
3. S o l v i n g t h e U n c o n s t r a i n e d M o d e l
It is convenient for us to begin by reviewing how the six algorithms studied in this paper
are implemented in the reversible investment version of the model. To have a consistent
terminology for discussing and comparing the algorithms, we use the framework in Reddy
(1993) s numerical analysis text, which corresponds closely to the framework presented in
Judd (1992a, 1992b). With one class of exceptions, the algorithms considered in this paper
are what Reddy calls
of Marcet (19SS)’s
The exceptions, standard implementations
A l g o r i t h m (PEA), fail to be
weighted residual
weighted-residual methods.
P aram eterized Expectations
methods only because of a technicality.
In the reversible investment case, there is only (2.4) with
h =
0 to solve, or, more
compactly,
P{h\ ,
foi
k >
, _
U c(hg{k.0),O)
r
0 )d 0 ' =
0,
(3-1)
0 and all 0 € 0 . In (3.1), the 0 argument in m reflects that in the reversible investment
case, the multiplier on gross investment,
E u le r r e s id u a l
h,
is identically zero. We refer to
function. Solving the model amounts to finding a function
functional equation,
R(g)
g
R(k,0;g)
as the
that solves the
= 0, he., sets the Euler residual function to zero everywhere.6
This problem is complicated by the fact that
k
can take on a continuum of values. This
implies that solving (3.1) is a problem of solving a continuum of equations (one for each
5
k,0)
in a c o n tin u u m o f u n k n o w n s (o n e
k1
v a lu e fo r e a c h
k, 0).
A p a r t fro m a few c a s e s , in w h ic h
R
has a convenient structure, solving this problem is computationally infeasible.
Instead, we select a function, g a, parameterized by a finite set of coefficients,
choose values for
a , a *,
to make
R ( g a)
a,
and
‘small’. Weighted-residual methods compute a* as
the solution to what Reddy (1993, p. 580) refers to as the
w e i g h te d - r e s i d u a l f o r m
of (3.1):7
(3.2)
The choice of weighting functions in (3.2),
practice, the range of integration over
k
w(k,6),
operationalizes the notion of ‘small’. In
in (3.2) is finite, with
k <
k <
k.
The boundary
points of this interval are chosen to ensure that the values of the capital stock generated by
the model always lie in the interior of the interval,
initial guess of
k
and
k
{k,k).
Computationally, we obtain an
by finding the interval ( k , k ) that contains
k0
and the ergodic set of
the log-linear approximation to the policy function in its interior.8
To apply the weighted-residual method, one has to select a family of approximating
functions. ga,a set of weighting functions, w{k, 0), and strategies for evaluating the integrals
in (3.2) and the integral implicit in the expectation operator in
R.
The procedures we
consider make different choices on these three dimensions. Two general types of g a functions
include spectral and finite element functions. In the former, each component of
ga(h\cr),
a
or </a(A\ —ct), over the whole range of
has influence over only a limited range of
a
influences
while in the latter, each component of
k
k 's .9
Regarding the weighting functions, a
necessary condition for (3.2) to pin down the parameter vector a, is that there be a number
of weighting functions equal to the dimension of
a.
We consider three types of weighting
functions. In one, the ic(k, 6)'s are related to the basis functions generating ga(k, 9), in which
case the algorithm is an example of the
G a le r k in
method. In another, the basis functions are
"In our case, (3.2) reduces to f k>0 R ( k , cr;ga) w(k, a ) d k + f k>0 R ( k , —<r;ga)w(k, —a ) d k = 0.
8See Christiano (1991, Appendix] for details about solving the model studied here using a log-linearization
method.
9See Judd (1992a,1992b), McGrattan (1993), and Reddy (1993) for more detailed discussions of spectral
and finite-element functions, respectively.
6
particular kinds of dirac delta functions, in which case the algorithm is an example of the
c o llo c a tio n
method. The versions of the PEA that we consider also choose an a to solve an
expression of the form (3.2) . 10 However, technically, standard implementations of the PEA
do not fall in the class of weighted-residual methods as defined by Reddy. This is because,
as we will see below, they work with weighting functions in which the parameter vector, a, is
an argument. Finally, two numerical procedures are used to evaluate the integrals in (3.2):
quadrature methods and Monte Carlo integration. We now turn to a detailed discussion of
the six algorithms considered.
3.1. Two S pectral, W eighted-R esidual M ethods
3.1.1. Parameterized Expectations
What distinguishes the class of Parameterized Expectations Algorithms is that they approx
imate the function, </, indirectly by approximating the conditional expectation in (3 . 1 ) as
fo llo w s:
J0I
where
ea(k.O)
g,0)pe>(0'
| 0 ) d 0 ' % exp[ea(£, 0 )],
(3 .3 )
is a function, parameterized by a finite set of parameters,
a.
The purpose of
the exponential in (3.3) is to guarantee nonnegativity. The PEA’s approximation of
obtained by solving
U c ( k . k ' J ) = , 3 e x p ( e a ( k , 0 ))
for
k'
given each
k,6,
U~'[-}
denotes the inverse function of
U c.
is
yielding:
£a( M ) = exp(0)A*° + (1 - 6)k - U - 1 [0exp (e„(M))],
where
g
(3.4)
The PEA approximation of
when the nonnegativity constraint on investment is ignored. We simply set
all A\0. It remains to describe how the various PEA’s go about computing
h,
/)„, is trivial
h a( k ,0 ) =
0 for
a*.
One way to view the three versions of the PEA that we consider, is to think of them as
sohing a particular fixed point problem. A given value of a induces, via (2.6), (3.4), and
10For a related discussion, see Marcet and Marshall (1994).
7
p e>,
a distribution on
for any fixed ( k , 9 ). A new set of parameters values, a',
m ( k , 9 , 9 ' - , g a, h a)
is found which makes
R ( k , 9 - , g ai) =
exp(eai(k,9))— f e, m ( k , 9 , 9 ' ] g a , h a ) p g i ( 9 ' |
zero in the sense of (3.2).11 Denote this mapping from a to a' by a' = ^(a;
N p ).
9 )d 9 '
close to
The value of
a selected by the PEA is the fixed point, a*, such that a* = S(a*; N p ). We computed a* by
applying a standard nonlinear equation solving routine.12 The versions of the PEA that we
consider differ in the form of the weighting functions used in (3.2) and in the computational
strategy for evaluating the integrals.
C o n v e n tio n a l P E A
In our implementation of what we call
c o n v e n tio n a l P E A
, we parameterize the expecta
tion function as follows:
ea(k ,9 ) = * j r
(3.5)
a t( 9 ) P t( v ( k ) ) ,
t=0
for 9
€ 0.
l.13 Here,
We define
The basis functions for
a
ea
are the
N p
Legendre polynomials, -P.(-), ? = 0, . . . , Arp-
is the 2 N p x 1 dimensional vector of parameters, {a,-(0), i
<p : ( k , k ) —*
=
0, . . . ,
N p - 1,9
€ 0}.
( —1,1) to ensure the polynomials in (3.5) are of similar orders of
magnitude. That is,
v-(*-) = 2 | 5 | - l .
(3.6)
11 Note, this definition of R coincides exactly w ith the one in (3.1), since the PEA policy function, ga,
implies U c(k\ ga{k, 0)i 0) = l3exp (ea(k\ 0 )).
12T he sta n d a rd im plem entation of the PEA finds a m by a m ethod of successive approxim ation, as the
limit of a, £((/; Ar^), S~{cr, Ar^ ) , . . . . As noted by Ju d d (1992b, chapter 13, pp. 11-14) and M arcet (1988),
this algorithm can yield explosive, oscillatory sequences, a, a ' , . . . , particularly for higher values of N p . One
altern ativ e is to instead ite ra te on th e o p erator 5 , where £*(0 ) = (1 —n)a -h /i£ (a; Arp), for a sm all fixed value
o f //. A problem w ith this approach is th a t it m ay require tim e-consum ing experim entation w ith altern ativ e
values of //. In our experience, th e equation-solving a ltern ativ e described in the te x t finds a* m ore quickly and
reliably. T he nonlinear equation solver we used is NLSYS in GAUSS. See M arcet (1988), M arshall (1992),
and M arcet and M arshall (1994) for a discussion of th e existence of a* and of the properties of exp[ea. (fe, 0)],
<7a*(fei0) as N p
00 .
13T he polynom ials, Pi , have dom ain and range ( —1,1), and are defined as follows. T he ith polynom ial is
P^x ) = 1 + q \x + . . . + q \x \ w ith th e a ’s defined by th e requirem ent Po(^) = 1 and
Pi(x)Pj(x)dx = 0
for j = 0 ......../ — 1 and / > 1. T he orthogonality p ro p erty of these polynom ials is designed to m itigate
m ulticollinearity problem s associated w ith step 2 of th e conventional PE A and P E A w ith exogenous oversam pling, which is discussed later in th e tex t. T his co nstruction of our basis polynom ials m ay m itigate
m ulticollinearity, but does not elim inate it, since th a t requires th a t the integrand in the above o rthogonality
condition be w eighted by a p robability density for x.
S
For any given value of a, a' is found by running a nonlinear regression of m ( h :, 0,
on the space of functions generated by exp (e~(k, 0 ) ) { Ov a €
W NP.
$'■ g a ,
%a)
To specify the regression,
we need to indicate how many observations of every possible type, (jfc,0,6>'), were used. This
is accomplished by specifying a density function, p(k, 0,0'-a). This density has the following
structure. Let P i(M ;a ) denote the marginal density of (k,0), which may depend on a.
Then, p{k,0,0'-a) = Pi{k,0;a)pe,(0> |0). The nonlinear regression is:
a’ =
argmin/M *, {exp(e~(£, 0 ) ) 3f2,VP
m ( k , 0 , 0 ' - g a, Aa)}2p(A-, 0 ,0'; c i) d k d 0 d 0 '
(3.7)
S(a;NP).
The hist oidei conditions associated with this regression are
JkAe, {e x P ( e ° ' ( ^ 0 ) )
- m ( k , 0 , 0 ' ; g a J i a )}
exp(ga, ( M ) ) ^ .( M )-P ( k , 0 , 0 ' ; a ) d k d 0 d 0 ' = 0,
da/
for / = 1, . . . , 2A> Taking into account the structure o f t h e fixed point,
(3.8)
is easily seen
to solve the version of (3.2) with weighting matrices
u-'(A-, 0: a-) =
p , { k , 0, a )
exp(e„.(fc,
^ ,
(3 .9 )
for / = 1, . . . ?2Arp.
Lnder conventional PEA, all three integrals in the function
S
lowing Monte Carlo method. First simulate a series of length T,
random number generator, and compute an initial value of
1. Simulate
{ k i , k 2,
given initial value,
..,A-r+ i} recursively using
{ 0 O, 0 U . . . , 0 T } ,
using a
Then:
k t+1 = g a ( k t , 0 t ), t =
0,1, . . . T and the
k0.
14,
PFA°T " m *
l° °btainia startinS value
l
a . 14
are evaluated by the fol
^
« is to generate the data according to step 1 of the conventional
baS6d °n l0S'Hnear aPPr0Ximatl0"’- d taki“*
9
value of
2. Find a', in the nonlinear least-squares regression problem : 15
a' =
1 T_1
argmin exp(e~(fct, 0 <)) - m (k t, 0 t , 0 t+i; £ o,fta)j
a€S 2NP 1
(3.10)
t=0
for 1 = 1 , . . . , 2 JVP.
For T large, the function being minimized in (3.10) coincides with the one being mini
mized in (3.7).
The PEA specification of the density of (fc, 6 ) concentrates observations on points of high
probability. Our com putational experim ents suggest that greater dispersion in ( k , 6 ) may be
desirable. Marcet and Marimon (1992) have made this observation in the context of a study
of the far-from-steady-state properties of a model. However, we find that this may be true
even when the objects of interest are properties of the steady state distribution of functions
of k. In our example, by increasing dispersion relative to conventional PEA , one gets a more
accurate estim ate of properties of the steady state distribution of k using a lower value of
T.
Presumably, this reflects the well-known fact that high variance in explanatory variables
implies greater precision in regression estim ates . 16 W ith these considerations in mind, we
studied two perturbations of the conventional PEA which imply greater dispersion in (&, 0).
PEA
with E x o g e n o u s O v e r s a m p l i n g
Under Marcet and Marimon’s (1992) P E A w i t h e x o g e n o u s o v e r s a m p l i n g , p i ( k , Q ; a ) is
modified so that extra mass is exogenously placed in particular regions of the state space,
15T he nonlinear least squares problem in step 2 was handled using a version of th e procedure applied in
M arshall (1992). Let 5 (a ; N p) denote S (a; N p ) w ith (3.10) in step 2 replaced by
l T~l
a'
2
= argm in — ^
e~(Art ,0*) - log f ?n(fct , 0tl 0*+i; ga,ha))\
Zest2NP 1 *=o L
v
=
0,
for / = 1 ,.. .,2 A rp. T h is is ju s t a linear regression and is easy to solve. We first com puted a** such th a t
a mM = S ( a “*; N p) using a non-linear equation solver. F inding a** takes little com puter tim e because of the
sim plicity of the m odified step 2. We then used a** as a sta rtin g value for solving a* = S(a*; N P). In our
experience, a* and a** are very close.
16In our context, th ere is an offsetting effect. Namely, w ith too m uch dispersion, accuracy of the p aram e
terized expectation in th e neighborhood of steady sta te is sacrificed.
10
( M ) - This m ethod is implemented by adding J terms to the criterion function in step 2
of the conventional PEA. These terms require simulating J sequences of length f each, of
the technology shock: { 6 { , t = 0 ,.. ., T - 1 } and of the capital stock, { k j , t = 1
j
= 1 ,..., J , where
for
is a value of the capital stock close to the region of interest. The
additional terms are:
T
E
j=i
-1
E {exp (ea'(ki,0ij) -m(ki,e{,ei+1-,gaj la)y
t=o
J
(3.11)
PEA-Collocation
Once the PEA is expressed as a weighted residual m ethod, it is clear that there are many
other ways to find a ’ . One could evaluate all integrals using one of a variety of available
quadrature formulas.1' Also, there are a variety of different weighting schemes that one can
use, some of which are discussed below. Finally, there are a great many alternative classes
of finite paiam eter functions that one can use to parameterize expectations.
Here, we pursue one particularly promising weighted residual m ethod. It works with a
more dispersed set of ( M ) ’s than does conventional PEA. It converts the nonlinear regres
sion in conventional PEA and PEA with exogenous oversampling into a linear regression on
an ort hogonal set of explanatory variables. There is reason to expect (and this is confirmed
in section 5) that the number of observations required in the regression is very small. Fi
nally. there is some a p r i o r i reason for believing that the method may have good accuracy
properties.
The method we pursue is a collocation m ethod, in which the weighting functions are dirac
delta functions. It is consistent with the use of either Monte Carlo or quadrature methods
to evaluate the integral in the definition of R . The dirac delta functions are constructed to
assign positive weight to the values of k corresponding to the N p zeros of the N p ' t h order
Q u a d ra tu re m ethods approxim ate integrals by th e w eighted sum of the integrand, evaluated a t a rel
a t e e h sm all num ber of points. T his approxim ation to integrals is known to be very accurate in the one
dim ensional case, and recently Ju d d and B ernardo (1994), applying the ideas of Stroud (1971), have argued
th at m ultidim ensional q u ad ratu re integration can also be m ade very efficient.
11
Chebyshev polynomial, TV p-18 That is, given a we compute a' to solve the following problem:
R (k,0',< h > ) =
for i
= 1 ,..., N p ,
6 e
exp ( e a . { k i , 0 ) ) — f m ( k i , O , 0 ' - , g a , h a ) p { e ' \ 9 ) d .O ' - 0
0. Here, A, = 9 _ 1 (ri), where T^p(r.) = 0,
i =
(3.12)
1 , . . . , N p . In addition,
we replaced Pi in ( 3 .5 ) with Chebyshev polynom ials, T,-, i = 0 , 1 , . . . , N p - 1. There is a
slight abuse of notation in ( 3 . 1 2 ), since R is also a function of a, but our notation does not
reflect this. The system of equations, (3.12), is (3.2) with the weighting functions, i v ( k , 0 ) ,
constructed using delta functions, 8 ( k — A:,),
N p
8(k uj' ( M )
l a)
as follows:
for 9 = a
(3.13)
= <
0
for 9 = —a
8 { k - k t)
for 0 — —a
0
for 0 = a
and
(3.14)
u?''(M) = <
i
= 1, . . . , Arp.
Our choice of Chebyshev polynomials as basis functions for the parameterized expec
tations and for determining the grid on k was influenced by the following two consider
ations. First, the discrete orthogonality property of Chebyshev polynomials greatly facili
tates the computations in (3.12) when N p is large .1819 This property implies that the mapping,
a1 -
S ( a ; N p ),
defined by the solution to (3.12), has a particularly simple analytic form:
ai(0)'
=
A
p
5 3 T/ { ^ ( k i ) ) log
i=1
f
m(ki, 9 , 0
g a, h a)p(6' \ 9)d9'
ue>
J
,
(3.15)
18T he C hebyshev polynom ials are defined as follows: T q (x ) = 1, T\(x) = x, and Ti(x) = 2xTi-\(x) —
Ti-i(x), for i > 2. T he dom ain and range of these polynom ials is ( —1,1).
19T he discrete orthog o n ality p ro p erty is th a t, for i j < N p :
Np
T
Ti(rk)Tj(rk) =
k=i
(
0 ,for i ^ j
<{ N p , for i = j = 0
[ N p / 2, for i = j
0,
where rk, k = 1 , . . . , N p are th e roo ts of T,vi>(-)-
12
for / = 0 ,1 , . . . , N p — 1, 9 € 0 . Here, fi = 2 for / > 0 and fi = 1 for / = 0. In obtaining (3.15),
we made use of the fact that (3.12) holds if, and only if, it holds for the log of the terms
on each side of the minus sign. The parameters in (3.15), a', is the set of coefficients in a
linear regression in which the explanatory variables are all orthogonal. As a result, there is
no multicollinearity problem, even if N p is quite large. For example, we have applied the
algorithm with N p as high as 100. In contrast, we had difficulty executing the regression
step in conventional PEA (see (3.10)) for N p larger than 5, because of multicollinearity
problems.20 Second, the Chebyshev interpolation theorem suggests that it is a good idea to
select grid points using the roots of a Chebyshev polynomial (see Judd (1992a, 1992b) for a
formal statem ent of the theorem). Suppose we have a given value of a, based on some fixed
value of N p . According to the Chebyshev polynomial approximation theorem, if N p —►oo in
the com putation of a ', then
sup
|| R ( k , 0; g a,)\\ -4 0 as l \ Tp -» oo,
ke(k.k).6£Q
(3.16)
when a' is given by (3.1o). Thus, the theorem suggests that with large N p, the function
exp ( e al( k , , 6 ) ) will not display pathological behavior between grid points.21 This is an at
tractive property that is not satisfied by polynomial interpolation schemes generally.
Other weighted residual methods can also be used to find ««. One such m ethod, Galerkin,
is discussed below.
We chose to go with collocation because it allows us to convert the
nonlinear regression step in the conventional PEA into a linear regression step. For example,
this conversion is not possible with Galerkin, which sets weighted averages of R to zero. This
impossibility reflects the fact that the log of an average is not equal to the average of the
log.
To summarize, the PEA can be viewed as a weighted residual method. In this context, the
2°den H aan and M arcet (1990) report sim ilar difficulties.
(3 ,2
^ T ? ! 311* t0 em Pllaf ze " h a t (3 1 6 ) does not say. Let us m ake th e dependence on a of R in equation
,x p „ c , bv u n lin g
T „«n, (3.16) does n o t s,.v
as A» -
oo. We are currently working on a proof of this la tte r proposition.
13
three versions of the PE A are differentiated according to the weighting functions used and the
manner of evaluating the integrals. The conventional PE A puts relatively heavy weight on
(jfc, 0 ) pairs with high probability and evaluates all integrals in the analysis using M onte Carlo
integration.22 The PEA with exogenous oversampling shifts more weight into exogenously
specified regions, but otherwise pursues the same com putational strategy as the conventional
PEA . We also described PEA-collocation. This appears to have several advantages relative
to conventional PEA: the nonlinear regression step with multicollinear explanatory variables
in conventional PEA is transformed into a linear regression with orthogonal explanatory
variables; the number of observations in the regression step is very small, and equals the
number of parameters in the parameterized expectation function, i.e., T = 2Arp, which is no
greater than 16 in our experiments (in conventional PE A , T can be in the tens of thousands);
and the distribution of ( k , 0 ) is more dispersed, thus ameliorating the conventional P E A ’s
problem that it tends to concentrate observations too much.
3 .1 .2 . G a lerk in
Judd (1992a) has discussed approximating policy functions by Chebyshev polynomials and
applying the Galerkin method. In our version of this approach, we proceed as follows.23 The
decision rules are:
;Y(0)-1
</(M)*&(M) =
£
a t( 0 )T t( v ( k ) ) .
i=0
(3.17)
The basis functions for g a are the N J Chebyshev polynomials. The [ N { c r ) + N ( - a ) ] x 1
vector
a = {ai(6)
| z = 0, 1 , . . . , N ( 6 ) —1 ,0 £ 0 }
22In his com m ent on conventional P E A , Ju d d (1993) expresses concern ab o u t the absence of a solid
rationale for sam pling at high probability p o ints, or for using M onte C arlo ra th e r th a n q u ad ratu re in teg ratio n .
O ur co m putational results in section 5 below have nothing to say ab o u t the la tte r p o in t, b u t do suggest th a t
sam pling at high p ro b ab ility p oints is inefficient.
23For a case stu d y com paring th e m ethod discussed in this subsection w ith a log-linearization procedure,
see C hari, C hristiano, and Kehoe (forthcom ing).
14
contains the as-yet-undetermined scalar coefficients and <^>(-) is defined in (3.6). For now, we
suppose that N ( a ) = N ( - a ) = N J . The 2A^weighting functions, w ( k , 6 ) , are constructed
from the basis functions. They are:
1
d g a (lc, 9 )
(1 - v?(fc)2) 1/2
da,
w ‘( k , 6 )
(3.18)
’
for / = 1
We evaluate (3.2) using M -p o in t Gauss-Chebyshev quadrature.
com pute the M
To do this, we first
> N J roots, r„ i = 1 , . . . , 3 / , of the M th order Chebyshev polynomial
and use these to construct a grid of capital stocks that is stored in the M x 1 vector k ,
k
= [^_1(?’i), v?_1(?-2) , . . . , ^ _1(rA/)]'. Second, we form the N J x M matrix A of rank N J :
To(r,)
T 0 ( r 2)
•
T o { vm )
Ti(ri)
T x( r 2)
• ••
T i ( i'm )
(3.19)
\ fJ- \
(7'l )
T N j - i { r 2)
■■■■
(r M )
lism g this notation, the Gauss-Chebyshev quadrature approximation of (3.2) is written com
pactly. in matrix form, as follows:
A R ( k , 0 ;a)
u lie ie R ( f c , 0 ; a ) =
[ R( k ' i , 0; a ) , R ( k 2 , 0;
= 0 , 0 <E 0 ,
a ) , . . . , i?(A'^/, 9; a) ] ' .
(3.20)
Equation (3.20) represents a
nonlinear system of 2 N J equations in the 2 N J unknowns, which can solved using standard
computational routines. Below, we refer to this method as Spectral-Galerkin.
3 .2 .
T w o F in ite E le m e n t, W e ig h te d R e s id u a l M e t h o d s
We consider the simplest class of finite element functions, those which approximate the
policy function with a g a ( k , 6 ) taken from the space of functions that is piecewise linear and
15
continuous in k for each fixed 8 24 The parameters of this function are the values of k' = g
at each point on a grid of N f capital stocks, for each 8. Denote this capital grid by a vector
K,
with elements ordered from sm allest to largest, K = (fci, k 2, . . . , k N j ) \ Here, k x > k and
k N j < Tc.
Also, denote the value of k ' at each ( k x, 8 ) by <z,-(0), for 8 € 0 , i = 1 ,2 , . . . , N * . The
x 1 vector a denotes the set of these parameters. The formal representation of g a is:
2N f
N>
& ( M ) = ;c«i(W(*),
i=i
are the
(3.2i)
functions,
follows:
Mi(k)
k-kt_j
kt-kt-i ’ k i . i < k < ki
= { kt+i—k
ki+i-ki' A,.2 ^ A. ^ A.2 i
elsewhere,
0,
for i = 2 , 3 , . . . Ar/ - 1,
Mi(k) =
k7-k
k2~ki ’ Ai < k < A'2
elsewhere,
0,
Bind
k~ k K , ~ l
Ms'j(k) = <
,
kNJ~kKl-1
kNj _ x < k <
k Nf
elsewhere.
0,
After specifying a set of 2A’^ weighting functions, i o ( k , 8 ) , equation (3.2) is used to pin
down values for a . An advantage of finite elem ent m ethods is com putational speed. The
fact that the parameters of the finite elem ent m ethod have only local impact implies that
the number of operations needed to solve this system of equations is smaller by orders
of magnitude than is the case in, for exam ple, the spectral m ethods.25 In the context of
the finite element m ethods with collocation, existing efforts to realize this com putational
24Reddy (1993) describes system atic procedures for expanding the space of finite elem ent functions to
include more th an one dim ension, and piecewise polynom ials of order higher th a n one.
25Let 7? denote the dim ension of a. By order of m agnitude of the operation count, we m ean an integer
j such th a t c( 7i ) / / 7J —‘a non-zero constant as n —►oo, w here c(n) is th e num ber of operations needed to
com pute a. For exam ple, j = 3 in th e S p ectral-G alerkin m eth o d because of th e m atrix inversion involved in
applying N ew ton-R aphson.
16
efficiency have centered on a particular t i m e s t e p p i n g algorithm (see, for example, Bizer and
Judd (1989), Coleman (1988), and Danthine and Donaldson (19S1)). In the context of finite
element methods with Galerkin, Judd (1991, p. 12) and McGrattan (1993) point out that
sparse matrix inversion m ethods can cut the order of m agnitude of the operation count.26
3 .2 .1 .
C o llo c a tio n
The finite elem ent, weighted residual method with collocation chooses a so that
R{ki,0;ga) = 0,
(3.22)
for i = 1 , 2 , . . . , A'' and 0 € 0 . This is (3.2), with the weighting functions, u>(A-,0), con
structed using dirac-delta functions analogous to those used for Quadrature PEA. Equation
(3.22)
is a nonlinear system of 2Ar' equations in 2Ar4 unknowns. Coleman and others apply
the following t i m e s t e p p i n g method for solving (3.22):
1. Fix a .
(We use starting values based on a log-linear approximation to the m odel’s
solution.)
2. For each element of the capital grid k, find the k ’( 0 ) that solves
U c( k h k'i ( 0 ) , 0 ) =
(3.23)
^ { ( } ) U c ( ^ { O ) . g a ( k ' ( 0 ) , a ) , a ) [ f k( k ' ( O y , a ) + 1 - 6 }
+ C 2 ) W , ( 0 ) , g a ( k ’( e ) , - a ) , - a ) [ f k( k ' ( e y , - a ) + l -
(5]},
for 6 = —<j, <r.
3. Set a' = {*,'(0), * = 1 , 2 , . . . N f , 6 e Q } .
4. If the maximum deviation of a' and a exceeds some chosen tolerance, set a = a' and
go to step 2. Otherwise, a 1 is the value of a sought.
r Z 1" ™
appHcatl0n ° f the finite eIem ent m ethod " i t h G alerkin, we did not apply a sparse m atrix inversion
17
Below, we refer to this algorithm as FEM-collocation.
3.2.2. Galerkin
The finite elem ent, weighted residual m ethod with Galerkin has been advocated by McGrattan (1993). In our exam ple, the m ethod works to select the value of a that solves the version
of ( 3 .2 ) in which the 2 Ar/ weighting functions, w ( k , 0 ) , are constructed from the basis func
tions, M i ( k ) , i = 1,2, . . . N ^ in a manner exactly analogous to (3.18) or (3.13)-(3.14). The
integral in (3.2) is approximated by a finite sum using M -p o in t, Gauss-Legendre quadrature.
The algorithm then solves the analog equations to (3.20) by Newton-Raphson m ethods .27
Below, we refer to this algorithm as FEM-Galerkin.
4. S o l v i n g t h e C o n s t r a i n e d M o d e l
This section describes modifications to the algorithms discussed in the previous section,
which are designed to accom m odate the irreversible investment version of our model. We
pursue two types of m odifications. One is based on the Euler equation associated with the
Lagrangian representation of the constrained problem.
The other is based on the Euler
equation associated with the penalty function representation of the problem. We apply the
Lagrangian method in the context of both spectral methods and FEM-Galerkin. For reasons
elaborated on below, we apply the penalty function method in the context of the FEM only.
The penalty function method is conceptually straightforward.
cussion.)
Since the Lagrangian m ethod is less straightforward, we begin this section by*i
27Taking in to account the region over which M i is zero, equation (3.2) is:
i=
(See Reddy for a dis
- 1 and / / ,+1 R(k ,0,ga)Mi(k)dk for i = 1, and
R(k,0,ga)Mi(k)dk
R(k,0 ,ga)M>(k)dk
for
for i = N * . These
expressions are (3.2) w ith w eighting functions defined analogous to (3.13)-(3.14), b u t th e dirac d e lta functions
replaced by Mi(k). AZ-point G auss-Legendre q u ad ratu re in teg ratio n of each integral involves selecting a grid
of M capital stocks, say £, , over th e associated range of in tegration using th e following procedure. F irst,
use the G auss-Legendre q u a d ra tu re form ulas (see Press, Teukolsky, V etterling, and F lannery (1992, pp. 140153)) to select M points in th e interval ( - 1 ,1 ) . T hen, th e elem ents of
are obtained by applying y - 1 (')
to these num bers. Second, we com pute the M N f —dim ensional vectors k = [l’i , . . . , kj^j) and R(k, 0; a). An
j\f x M N * m a trix A is com puted so th a t AR(k, 0; a) represents th e finite sum approxim ation of th e integral
in (3.2).
18
presenting our basic strategy for applying it. Under the Lagrangian m ethod, we seek two
functions: one relating the capital decision to the state and the other relating the Lagrange
multiplier to the state. It is com putationally efficient to restrict the space in which these
functions belong. We im pose two types of restrictions. The first type is valid generally,
and simply enforces the Kuhn-Tucker conditions. The second type of restriction involves
assumptions about the properties of the exact functions that solve the problem. We assume
that ( 1 ) the irreversibility constraint is never binding for the high value of the shock, (ii)
the capital policy function is continuous, (iii) the multiplier policy function is continuous,
and (iv) for fixed 6>, if the constraint is binding from some level of k then it is also binding
for all higher levels of k. In practice, the validity of these assumptions can be verified e x
post
by studying the Euler residual function associated with a proposed numerical solution.
We do this in section 5, where we report our numerical results. There we also evaluate the
validity of our assumptions by studying the solution to our problem based on a dynamic
P rogra mini ng a lgori t hm.
Figure 1 depicts hypothetical policy functions in k x k' space. The policy functions are
drawn for a case in which gross investment has n o t been constrained to be nonnegative, and
m fact does take on negative values for some values of k. The points at which the g ( k , a )
and g ( k , - a ) functions cross the 45° line mark the ergodic set of capital for this case. Notice
that g ( k . a ) never crosses the (1 - 6 ) k line within the ergodic set. This implies that when
°
= ^ inVCStment is »evei' n a t i v e . If. when we im pose the nonnegativity constraint, g ( k , a )
retains its basic shape, then we can infer that A will always be zero when 0 = a . On the
other hand, we see from the figure that g ( k , - a ) crosses the (1 - S ) k line at a point within
the ergodic set. Based on this result, we conjecture that when we im pose the nonnegativity
constraint, the exact policy function has the following property: when k > k , for some k,
and 6 = - a the nonnegativity constraint binds and A > 0, and when k < k and 6 = - a the
nonnegativity constraint does not bind and A = 0.
We use these observations and the requirement that the Kuhn-Tucker conditions be
satisfied to constrain the space of approximate policy functions. Our task is to find functions
19
to approximate two policy functions: g ( - ) as before, and the function determining A, h ( - ) .
We restrict the space of approximating functions for g ( - ) as follows:
g(k,a)
(4.1)
« g a(k,(r),
max{<7a(A:), (1 — £)A’},
k < k
(1 - S)k,
k < k
g ( k , - a ) & g a( k , - c r ) = <
(4.2)
where % is a function of a and is defined by the property
g a( k)
(4.3)
= (1 - S ) k .
We restrict the space of approximating functions for /*(•) as follows:
h(k,cr)
(4.4)
« 7ia ( k , a ) = 0 ,
and
h ( k , —a )
« h a{k,
a)
0,
k< k
m a x { /ia(A'),0},
k < k
= <
(4.5)
with the property
h{k) = 0.
(4.6)
The first expression just says that we use some convenient function ya(-,<r), i . e . a poly
nomial or a piecewise linear function parameterized by the vector a , to approximate the true
rule g { - , c r ) . Expression (4.2) embodies the restriction that the approximating function for
—a ) , g a ( •,
—a ) , must force gross investm ent to be zero for k > k . A polynom ial or a piece-
wise linear function g a( •) is used to approximate the true policy function for k < k . The max
operator in (4.2) ensures the constraint on non-negative gross investm ent is never violated.
Also, by solving for k using (4.3) we make sure g a ( - , - c r ) is continuous. In (4.4) we use our
conjecture that the gross investment constraint never binds when 9 = cr to force A to be zero
for all k in this case. For the case 6 = —a we force A to be zero for k < k and we propose a
20
polynomial or piecewise linear function h a (-) be used to approximate h ( •, - a ) for k > k . The
max operator is used to ensure the Lagrange multiplier is always non-negative and condition
(4.6) forces the function that approximates h ( •, - a ) , h a (-, - a ) , to be continuous.
It should be clear from expressions (4.1)-(4.6) that our space of approximating functions
forces the Kuhn-Tucker conditions to be satisfied exactly. Our only task is to find approxi
mating functions within this class that set the Euler residuals to zero for admissible k and
This task is not unlike the one we encountered when we ignored constraint (2.2), which
6.
suggests we can apply similar methods to solve the constrained model.
Assumptions (i)-(iv), stated at the beginning of this section, play an important role
in our Lagrange multiplier procedure.
In section 5 below, we report evidence that these
assumptions aie valid in our model, so that imposing them is innocuous in our application.
However, they may be difficult to impose in higher dimensions, or there may be models in
which the assumptions are not valid. Because of this, it is useful to note that there exist
versions of our Lagrange multiplier method which do not involve assum ptions (i)-(iv). For
example, we can modify our procedure so that (iii)-(iv) are dropped by replacing (4.2)-(4.3)
by g { k , - a ) ss g a ( k , - a ) = m a x { g a { k ) , { l - S ) k ) for all k > 0, and (4.5)-(4.6) by
H t .-*) * U < ' . -a) = |
>(l-6)k
_°’
l max{/?a(A’), 0},
g a(k, - a )
< (1 - S ) k
Tins perturbation on our method does not involve computing the variable, k. We have
conducted several experiments with this procedure and found it to be practical.
We now turn to the description of our modifications to the algorithms discussed in the
previous section.
4 .1 .
A L a g ra n g ia n M o d ific a tio n to t h e P E A
The modification to the PEA to accommodate the case where (2.2) binds occasionally is
remarkably straightforward. We must now be careful to allow the function m in (2.6) to
accommodate a potentially nonzero multiplier, h. To do this, we work with the following
21
modified version of (3.3):
[m(k,6,6'-,g,h)p0'(0'|0)dd'taexp[ea(M)]>
(4.7)
J6'
where e a ( k , 0 ) is defined in (3.5) and m is defined in (2.6). The P E A ’s approximation to the
decision rule is:
g a ( k, 6)
= m ax { ( 1 - 6 ) k , e x p ( 6 ) k a + (1 - 6 ) k - U ~ l \fi exp (ea(fc, 0 ))]} ,
(4.8)
and its approximation to the multiplier function, h a ( k , 0 ) , is:
M M ) = U c ( k , g a { k , 0 ) , 0 ) - /?exp[ea( M ) ] -
(4 *9)
W ith these modifications to g a and M the three versions of the PEA can be implemented
as described in the previous section.
4 .2 .
A L a g ra n g ia n M o d ific a tio n to t h e S p e c tr a l-G a le r k in M e t h o d
We choose functional forms for </a(-,cr), g a {- ) i and /».„(•) as follows:
A'(CT)-1
«;(<7 )r ;M k ) ) i
(4.10)
ai(~<r)Ti(<p(k)),
(4.11)
g a{ k , < r ) =
1= 0
Ar(—cr)—1
9 »(*)
=
H
t=0
Nx
h , ( k ) = ' £ b iT M k ) ) -
(4.12)
i=0
An Euler residual function can be constructed in the manner used before to form R ( k , 0 ; a ) ,
where a is the [Ar(cr) + N ( - c r ) + N x] x 1 vector of unknown polynomial coefficients.
To apply the Spectral-Galerkin m ethod we m ust find a grid vector k and weighting matrix
.4 that can be used to form the system A R ( k , 6 ' , a ) . We can use a Chebyshev grid as before
22
to construct k. When constructing the grid we made sure it was fine enough so that there
were a substantial number of grid points to the right of k. Presumably this is needed to
ensure a good approximation t o h( - ) .
:
We set N ( c r ) = N ( —a ) + l \ TX = N J, N x = N ( —a ) — 1, and selected M , the number of
elements in fc, such that M > N J . We construct the matrix A as shown in (3.19). The
approximation problem is then identical to the one described before: find a that solves
A R ( k , 6 ;a)
= 0, for 6 G O .28
Our Lagrangian modification of Spectral-Galerkin accom m odates nondifferentiable deci
sion rules. This seems appropriate in problems with occasionally binding constraints. We
found it less convenient to accom m odate nondifferentiable decision rules in the context of a
penalty function version of Spectral-Galerkin, so we did not pursue this.
4 .3 . A L a g ra n g ia n M o d ific a tio n to th e F E M -C o llo c a tio n M e th o d
This section describes how we applied the Lagrangian method to FEM -collocation.29 We
choose piecewise linear functions to form g a {- , <r) , <)„(•), and h a (-) and select the capital stock
grid k — ( k i , k 2 , ■.., k j v f ) ' .
The objective is to solve for the coefficients associated with
this grid: a , ( Q) , i = 1 ,2 , . . . , A'G 0 G 0 , as before, and 6,-, i = 1,2 , . . . , i \ r/, where each &,•
corresponds to the value taken by the Lagrange multiplier at the Pth element of k when
Stack the undetermined coefficients in the vector
0 = —a .
a
= (a!(cr),« 2 (or),. . . , a N / ( c r ) , a i ( - a ) , a 2( - a ) , . . . , a N f ( - c r ) , b u b 2, . . . , bN } )'.
28If N ( a ) ^ N ( — ct) + j V \ we can use the M x 1 Chebyshev grid as in the previous case, b u t now we
m ust choose separate w eighting m atrices for the Euler residual functions R(k,cr\a) and R(k, -a\a). These
weighting m atrices, call them A h and A 1, can be constructed in a m anner analogous to the construction of A
in the previous case. T he appro x im ation problem is to find a th a t solves A h R { k , <7 ; a) = A 1R(k, - a ; a) = 0.
29See Colem an, Gilles, and Labadie (1992) for another application of the Lagrangian m ethod in th e context
of th e FEM -collocation algorithm .
23
We m odify the FEM -collocation algorithm as follows. In step 2 of that algorithm equation
(3.23) for 9 = a is replaced by
U c(ki, k \ ( a ) , < j ) =
P{(h)Uc(m<r),ga(ki(<r),<r),<r)[fk(ki(<r),<r) + 1 - $]
(4.13)
+ ( l W c ( k ' i ( a ) , g a (k'i ( a ) , - a ) , -<r)[/jt(^(cr), -cr) + 1 - 8]
- h a( k > ( a ) , - * ) ( 1 - 8 ) ) } .
Equation (3.23) for 9 = —a is replaced by:
U c{ k i ,
A’t' (-c r ), — cr) - A, =
/?{ (|)^ (A -1/( - a ) , ^ Q( ^ ( - < 7 ) ,<T),cr)[/jt(A-'(-cr),cr)+ 1 - 8]
(4.14)
H l W c ( k ' ( - a ) , g a( K ( - ^
-<r), - ° ) l f k W { - * ) , -<r) + l -
«]
- h a (k't ( - a ) , - a ) ( l - 8 ) ) } .
For each i equation (4.13) is solved by choice of k [ [ a ) as before. Equation (4.14) is first
solved by choice of A*'(—a ) with A; = 0 . If A’-(—cr) > ( 1 — S) k { then we proceed to the next
value of i in the sequence i = 1 ,2 , . . . , N * . Otherwise, k'{ ( —cr) is set equal to ( 1 — S) k { and
(4.14) is solved by choice of A;. The only other modification to the algorithm is to add to
the updating rules of step 3 the conditions 6; = A,-, z = 1 , 2 , . . . , ArT 30
4 .4 . A P e n a lt y F u n c tio n M o d ific a tio n to t h e F E M -G a le r k in M e t h o d
We now turn to the penalty function im plem entation of FEM-Galerkin. This is a modified
version of the algorithm applied by M cGrattan (1993). In this approach, a penalty is applied
to violations of the constraint on capital accum ulation. Specifically, we solve a modification
30T here are two technical m a tte rs to be resolved regarding the im plem entation of side conditions (4.3)
and (4.6). T h e function g is defined as a set of linear segm ents which are joined a t ga(ki) = ai(— (T) for
/ = 1 , 2 , . . . , N f and at ga(k) = (1 —S)k, where k is o b tained as follows. Exam ine ai(— a) for i = 1 , 2 , . . . , N f
u ntil th e first /, say ?' occurs w ith a t (—cr) < (1 — S)ki. T hen use the line segm ent defined by ga(kii_2)
and ga{kii_ :l) to linearly ex trap o late a value for k. Form ally, g is defined exactly as g is in (3.21), w ith the
exception th a t k is added in to the list ki, Ao,. . . , kN f . We im pose (4.6) by defining h to be com posed of
linear segm ents joined at h a(k) and at h a(ki) = bi for i = l , 2 , . . . , A r^, w ith h a(k) = 0.
24
to the original planner’s problem as follows:
w (k, 0) =
m ax
(c(*, k \
0))
+ 0 1
W ( k ' , 0 ' ) Pe ,(0> \ 9 ) d 9 ' -
| [max{0, (1 - 6 ) k - * ' } ] 2 .
(4.15)
Here tt is a nonnegative penalty parameter. For tt = 0 , (4.15) describes the problem for the
model when the gross investment constraint is ignored. For positive tt, violations of (2.2)
reduce the planner’s objective function. Intuitively, we might expect that for large values of
tt
the solution to problem (4.15) would be “close” to the exact solution of the constrained
problem.31
We apply the penalty function method by solving the sequence of problems corresponding
to an increasing sequence {^ n}. In a typical experim ent, the sequence contained 31 elements
beginning with 1, 2, 10, 20, 50, ... and ending with 1,200. For each value of trn, it is necessary
to solve, using the FEM-Galerkin m ethod, the Euler equation associated with (4.15):
U c { k , g ( k , 9 ) , 0) -
7rn m ax{0, ( 1 - 8 ) k - g ( k , 9 ) }
- P !8'{L!c{<Ak.o),g{g(k,6),o'),o,)[fk(g{k,e),ei) + (i - s)]
~ ( l - 8 ) 7 T n max[0. (1 - 6 ) g ( k , 0 ) - g ( g ( k , 0 ) , 0 ' ) ] } p e, { 0 ' | 9 ) c W = 0.
The algorithm stops when the maximum violation of the gross investment constraint on the
capital stock grid is smaller than some prespecified tolerance.32 Denote by tt* the value of
the penalty parameter when the algorithm is completed. Then, following Luenberger (1969,
Luenberger (1969, T heorem 1, p . 306) provides a theorem for the case where the solution to a constrained
m axim ization problem is a finite dim ensional vector. In this case, solutions to the penalty function version
con e r ^ f tim C° rrefSp°!ldmS
an lncreasinS sequence of penalty param eters tending tow ard infinity will
Umably k " S traightf° rWard t0 CXtend the ^ e o r e n f t o our e n v h o n m e *
=
tim e intensive
j
11 te n s n e
~ W * ~ MW)]}.
r hen " V3lU; ° f ffn iS encountered w ith 7 ( ^ ) less th an the chosen tolerance.
uq
d ! C°[reSp0ndlng t0 an lncreasing sequence {trn } is com puter and program m er\ \ hen we instead a ttem p ted to in itiate th e algorithm w ith a large value of tt however the
r o T o t l l ' ^ r ” ' our l,,i,ial gu“ s at u,e soiuiio'’ (,,an,"! ’the
25
Theorem 2, p. 307), an approximation to h(k,6) is given by:
h a(k, 6) =
7r* m ax{0, (1 - 8 ) k - g a ( k , 0)}.
It is worth noting that the FEM is particularly suited to working with penalty functions.
This is because it uses a functional form that easily accom m odates nondifferentiabilities in
the exact policy function. As we shall see below, there is a likelihood that the solution to
the constrained problem will involve a policy function for capital that is not sm ooth. If we
were to apply the Spectral-Galerkin method with a penalty function we would encounter
difficulties because that m ethod attem pts to approximate the exact policy function with a
globally sm ooth approximator.
5. E v a l u a t i n g t h e A l g o r i t h m s
The algorithms we have described were used to approximate the solution to a particular pa
rameterization of the model. In addition to applying these algorithms we also approximated
the model solution using dynam ic programming ( DP) applied to a discrete version of the
model. We take the DP solution to be virtually identical to the exact solution and use it
as a benchmark for evaluating the algorithms discussed in this paper. Details about these
computations are reported in appendix 1. One of our findings is that the results of all the
algorithms are reasonably accurate. The differences between solution m ethods are small and
not economically very meaningful in the context of our model economy. Nonetheless, there
are some noticeable differences in accuracy and in com putation tim e, and we think these
are potentially useful as input into decisions about which algorithm to use in more complex
modeling environments.
We study three aspects of the approximate solutions: the Euler residuals, the policy
functions, and the implications of the policy functions for various first and second moment
properties of the model. W ith one exception, the parameter values we chose are standard
in the real business cycle literature and are as follows: a = 0.3, 8 = 0 .02, f3 = 1.03- 25, and
26
a
— 0.22. The exceptional case, cr, was chosen large enough to ensure that the investment
constraint binds occasionally. Finally, we specified U ( c t ) = Inc*.
We study three cases for the PEA and two cases for each of the other algorithms. For
the PEA, A p
3 in all cases except PEA-collocation applied to the irreversible investment
model, in which case N p = 8. The solution labelled N p = 3 corresponds to what we have
called conventional PEA, while the one labelled N p = 3* is PEA with exogenous oversampling
starting near an estim ate of k . The solutions labelled N p = 3** and N p = 8** correspond to
PEA-collocation applied in the reversible and irreversible investment cases, respectively.33
For the Spectral-Galerkin algorithm the cases are: N J = 3 and 8 in the reversible investment
case and (Ar(<r), Ar(-< j), N x ) = (5, 3, 2) and ( 9, 5, 4) in the irreversible investment case. For
the FEM -collocation algorithm the two cases are: A^ = 36 and 72. In the case of the FEMGalerkin algorithm we stud}- solutions based on N * = 18 and 36.34 For each algorithm, the
highest order parameterization reported is the one for which we obtained convergence. Our
convergence criterion was based on the second moment properties reported in tables 2 and
3. For each method we incremented the number of parameters until the change in all second
moment properties was less than one standard deviation. Typically, the last moments to
converge were the ones based on financial variables.
All com putations were carried out on a Gateway 2000 486 D X 2/66 and the programs
33For (lie conventional and m odified PEA algorithm s, we set T = 10,000. T his com pares to a value of
2,500 used by den Haan and M arcel (1990). T hey work with a m odel sim ilar to ours and assum e
the technology shock stan d ard deviation is 0.32, which co ntrasts with a stan d ard deviation of 0.22 in our
m odel (the one-step-ahead conditional stan d ard deviation in the technology shock in their m odel is 0.10).
For conventional PE A w ith oversam pling, we set T = 7,500, J = 100, T = 25, and k{ corresponding to a
num ber (34.0) in the neighborhood of k, for j = 1,. . J.
T -
34A dditional details of how the algorithm s were im plem ented are as follows. For the Spectral-G alerkin
cases we used M = 100. For th e p aram eterization of th e m odel we exam ined, this guaran tees an am ple
supply of grid points on either side of k. T he grids for th e FEM algorithm s were chosen to be equally spaced
bet ween boundaries ju s t outside th e initial guess for th e ergodic set, (k, k). T h e tolerance on violations of the
investm ent co n strain t for the p en alty function version of the FEM -G alerkin algorithm was set a t 5 x 1 0 -5
T his tolerance was reached for tt* = 1,250 w ith N * = 18 and
= 1,200 w ith N * = 36. In th e text,
the approxim ate policy functions were expressed in term s of the level of th e capital stock. VVe did th is to
simplify the presentation. In th e calculations we work in term s of the log of the capital stock. G rids for
the Spectral-G alerkin and PEA -collocation m ethods were constructed based on the log of th e cap ital stock
while grids for the FEM were co nstructed based on the level of the capital stock.
27
are available on request. The com putation tim es are displayed in table l . 35 These tim es
should be interpreted with caution. First, we did not make extensive efforts to optim ize the
computer code. Second, it could be misleading to extrapolate the relative tim e requirements
reported in our experim ents to larger problems. For example, the technology shock in our
model can take on two values only. This biases com putational tim es in favor of m ethods such
as Spectral-Galerkin and PEA-collocation which exploit this fact, and against conventional
PEA which does not. Also, the operation counts of the various algorithms are of different
orders of magnitude. For exam ple, the Spectral-Galerkin algorithm involves a number of
operations that grows at the rate of the cube of the number of decision rule parameters
sought, while the operation count for the FEM grows linearly or with the square.36
Still, there are several observations worth making about the tim ing of the various al
gorithms. First, note that PEA-collocation is faster by orders of m agnitude than all the
other algorithms. Second, in contrast with the other algorithms, the com putational tim e
for the P E A ’s does not increase substantially when the irreversible investm ent constraint is
imposed. This reflects the fact that, in contrast to the other algorithms, taking account of
this constraint adds virtually no com putational burden to a PEA. Third, the slowest algo
rithm applied to the irreversible investm ent model is FEM-Galerkin. This reflects the fact
that this algorithm involves repeatedly solving the model for higher values of the penalty
function parameter.
35W ith three exceptions, all the program m ing was clone in GAUSS. Tw o exceptions were the conventional
PE A and PE A w ith exogenous oversam pling, which com bined FO R T R A N w ith GAUSS. T h e sim ulation p a rt
o f the PEA was program m ed in FO R T R A N and im p o rted in to a GAUSS shell. FO R T R A N program m ing
could have reduced th e com putation tim es for the FE M -collocation algorithm significantly and th e com pu
ta tio n tim es for th e FE M -G alerkin alg orithm m arginally. T h e second exception is th e D P calculations which
were done in F O R T R A N . T he D P algorithm s required hours to achieve convergence, b u t th e exact tim es are
not reported on th e table.
36C hristiano and Fitzgerald (1991) elab o rate on this observation and illu strate it in a com parison of
Spectral-G alerkin and FEM -collocation. It should be noted th a t to achieve a given degree of accuracy in
problem s w ith sm o o th decision rules, fewer param eters m ay be required in th e S pectral-G alerkin proce
dure th a n in the FE M . Also, the assertion ab o u t th e ra te of grow th of th e operatio n count in SpectralG alerkin assum es th a t a sta n d a rd N ew ton-R aphson equation-solving m ethod is used. T h e operation count
for Spectral-G alerkin could be reduced if a m ore sophisticated version of th is algorithm were used.
28
5 . 1.
Euler Equation Residuals
Here, we focus on the Euler residual function (ER F), defined by R ( k , 9; a*), where a* denotes
the solved value of a . We consider the E R F ’s for both the reversible and irreversible invest
ment versions of the model. We study the graphs of the residual functions (see figures 2a and
2 b) and the maximum absolute value (MAV) of
R ( k , 6 ; a*)
over k <E ( k , k ) , for 0 = a and - a ,
respectively (see table 1). We found it useful, in the context of the PEA , to also compute
MAV’s over a narrower interval, containing 90% of the simulated capital stocks. The upper
and lower boundaries of this set are indicated by vertical lines in the PEA component of
figure 2a. For the PEA, MAV numbers not in parentheses in table 1 are based on this 90%
confidence interval, while numbers in parentheses are based on the entire interval, ( k , k ) . In
figures 2-5, results are displayed for a range of capital stocks corresponding roughly to the
set, ( k , k ) = (22.0,40.0).
We now consider the two left hand columns of figures 2a and 2b and the top panel of table
1, which peitain to the ERF s of the reversible investment model. We begin bv summarizing
our findings for the PEA. We found that increasing N p beyond N p = 3 in the context
of conventional PEA has relatively little impact on the results. In the interests of saving
space, we do not document this finding here. Essentially, results based on conventional PEA
converged at A'p = 2 and roughly correspond to our A rp = 3 findings.37 Note from figure
2a that the ERF's for conventional PEA and 0 = a are consistently negative over the 90%
region of capital stocks. This is a sign of inaccuracy in the conventional P E A ’s solution.38
One way to accom m odate this is to increase the length of the Monte Carlo simulations, T .
However, we found that unpractically large values for T are required to achieve a significant
degree of improvement in accuracy.39
Increasing the order of approxim ation to N p — 5 does not change the results significantly. We found it
difficult to increase N p above 5.
38A t th e sam e tim e, the deviations from zero in th e euler errors are not large by one economic m easure.
T he percent increase in consum ption which would move current m arginal utility of consum ption down enough
to close a given gap, R, in the euler error is 1OO0CR/(1 - f3cR), where c is th e level of consum ption. T he
worst MAV in the 90% confidence region for the P E A is th e value of 0.00013 reported for conventional PEA ,
0 = ~o- la consum ption units this is about 0.03 percent.
39We did the N p = 3 calculations for T = 20,000, 40,000, 60,000, and 80,000. For 9 = <t , the MAV’s
29
An alternative to increasing T is to alter the distribution of ( k , 6 ) points at which the
computations are done. This can be seen by noting the significant improvements that are
obtained in the N p — 3* and N p = 3** versions of the PEA . In particular, note that the
E R F’s for 6 = a are closer to zero when N p = 3* or 3**, and by orders of magnitude in
the latter case. We infer from these results that conventional P E A ’s Monte Carlo procedure
for selecting points in the state space for the com putations is not optim al. The alternative
procedures, PEA with exogenous oversampling and PEA -collocation, seem to work better.
In each case the distribution of ( k , 0 ) values used in the calculations is more diffuse relative
to that in conventional PEA . We conjecture that this is the basic reason that they do better.
The idea is that they do better for the same reason that regression coefficients are more
precisely estim ated, the greater is the dispersion in the explanatory variables.
Notice from figure 2 that the performance of the conventional PEA and PEA with exoge
nous oversampling deteriorates significantly in the outer 5% tail areas of the interval ( k , k ) .
The dramatic improvement evident with PEA-collocation over the entire interval ( k , k ) is
quite striking in comparison (visually, it is hard to distinguish from zero in the figure), espe
cially since this improvement is achieved by requiring only that the Euler residuals be zero
at three points in this interval.
For the other three solution algorithms, increasing the number of parameters in the
decision rule is very effective at driving the E R F’s toward zero. In each case, convergence to
zero is roughly uniform over the range ( k , k ) . Note how sm ooth the E R F’s corresponding to
Spectral-Galerkin are, in contrast to those based on the two FEM m ethods. This reflects the
fact that, in our example, the sm oothness in the Spectral decision rule mimics more closely
the properties of the exact decision rule.
We now consider the two right hand columns of figure 2 and the bottom panel of table
1, which pertain to the E R F ’s of the irreversible investm ent m odel. We again begin by sum
for these cases are 9.5 x 1 0 "5 (9.5 x 1 0 "5), 7.4 x 1 0 "5 (3.0 x 10~4), 5.3 x 10~5 (4.5 x 1 0 "4), 3.2 x 10“ 5
(3.4 x 1 0 -4). respectively. For 0 = - a th e M AV’s are 1.0 x 1 0 "4 (1.0 x 1 0 "3), 2.6 x 10“ 5, (1.2 x 1 0 "4),
2.3 x 10_5(1.5 x 10- 4 ), 3.5 x 10-5 (9.7 x 10- 5 ). As in table 1, n u m bers not in parentheses correspond to
M AV’s based on an interval containing 90% of th e realizations of th e cap ital stock. N um bers in parentheses
correspond to MAV's based on an interval th a t contains 100% of th e realizations.
30
marizing our findings for the PEAs. The performance of the N p = 3 and 3* versions of this
algorithm is roughly comparable to what it is in the reversible investm ent case. In particular,
we found that increasing N p does not contribute much to accuracy in the conventional PEA ,
but PEA with exogenous oversampling does help.40 The PEA-collocation ERFs deteriorate
somewhat (particularly around the point at which the investm ent constraint begins to bind:
k — k =
33.40, 6 = —o , according to D P), but still dom inate the other im plementations of
the PEA.
We now summarize our results for the other algorithms. Relative to the reversible in
vestment case, the two Galerkin methods have a harder tim e driving the E R F’s zero. The
MAV for the best Spectral-Galerkin solution (i.e., ( 9, 5, 4) ) is approximately 9 x 10- 5 , as
opposed to 4 x 10“7 in the reversible investment case.41 Similarly, the MAV for the best
FEM-Galerkin m ethod is 4 x 10- 4 , as opposed to 7 x 10-5 in the reversible investment case.
By contrast, the rate of convergence for FEM -collocation is comparable across the reversible
and irreversible investm ent models. The reason the Galerkin m ethods have problems is sim
ilar to that emphasized in the case of PEA-collocation. It has to do with the difficulty they
have in driving the Euler residuals to zero in the neighborhood of k = k.
5.2 . P o lic y F u n c tio n C o m p a r iso n s
We now com paie the polic\ function approximations obtained for the two versions of the
model and for the four solution algorithms. In figure 3, policy and multiplier functions based
on the highest order Spectral-Galerkin m ethod are compared with those based on the DP
40We did the N p - 3 calculations in the irreversible investm ent m odel for T = 20,000, 40,000, 60,000,
and 60, 000. For # = <t, the M AV’s for these cases are 1.1 x 10-4 (2.3 x 10- 4 ), 8.5 x 10-5 (5.2 x 10~4),
6.1 x 10-5 (6.6 x 10- 4 ), and 3.9 x 10~5 (5.4 x 10- 4 ), respectively. For 9 = -c r the MAV’s are 1 2 x 10-4
(2.0 x 1 0 -4), 8.2 x 1 0 - 5 (7.8 x 10~4), 8.2 x 1 0 "5 (7.2 x 10~4), and 9.4 x 10~5 (7.9 x 10~4). As in table 1,
num bers not in parentheses correspond to MAV’s based on an interval containing 90% of th e realizations of
th e capital stock. N um bers in parentheses correspond to M AV’s based on an interval th a t contains 100% of
the realizations.
We considered higher degree approxim ations for the Spectral-G alerkin m ethod b u t were unable to achieve
anything resem bling the convergence to zero noted in the reversible investm ent case. For exam ple, w ith
A' = 8 in the reversible investm ent case, the MAVs are ab o u t 2 x 1 0 " 11 and 7 x 1 0 "11 for 6 = a and
0 = -cr, respectively, while in the irreversible investm ent case w ith (1 5 ,8 ,7 ) the MAVs are ab out 4 x 1 0 "5
and 5 x 10 5 for 6 = cr and 6 = — cr, respectively
31
m ethod, for both reversible and irreversible investm ent versions of the m odel.
There are three main results in figure 3. First, there is very little difference between the
solution based on Spectral-Galerkin and DP. Over m ost of the range of k the functions are
identical to the eye. We infer from this that the Spectral-Galerkin m ethod provides a highly
accurate approximation to the solution. Second, the shape of the policy and multiplier func
tions validate the four assumptions we made when constructing the space of approximating
functions in section 4. Third, for k < k the constrained and unconstrained policy functions
are virtually identical, while investm ent is (slightly) lower in the constrained economy with
0
== cr, when k > k. Presumably, this reflects in part a rate of return effect: the payoff from
capital investment is lower in the irreversible investment economy, since there is some chance
that O' = —<7, in which case the lim itation against consuming capital is binding. The net
effect of the irreversibility constraint on the average capital stock is quite small, since the
impact on investment of the irreversibility constraint is positive when k > k and 6 = —<7.42
The FEM decision rules are indistinguishable from Spectral-Galerkin, and so we do not
graph them. It is worth comparing the PEA and Spectral-Galerkin rules, however. Figure
4 compares the PEA and highest order approximation Spectral-Galerkin investment policy
functions for the reversible investm ent economy. Qualitatively, the findings here are con
sistent with our analysis of the Euler residual functions. First, the PEA with exogenous
oversampling appears to do better in the 90% confidence region for capital, than conven
tional PEA. Also, the greatest inaccuracy in the conventional PEA and PEA with exogenous
oversampling appears to be in the lower tail of the capital stock distribution. Finally, the
PEA-collocation policy functions appears to be indistinguishable from the Spectral-Galerkin
rule.
Figure 5 compares the PEA and the highest order approximation Spectral-Galerkin in
vestm ent policy functions for the irreversible investm ent economy. Again, the results here
are consistent with our analysis of the Euler residuals in figure 4. Thus, the PEA-collocation
42\Ye found th a t th e average capital stock in the reversible and irreversible versions of the m odel is 31.3
in each case. U ncertainty per se does seem to have an im pact since the steady s ta te cap ital stock is 30.5.
32
decision rule appears to be essentially indistinguishable from what we take to be more or
less the exact decision rule. Also, the decision rule produced by the PEA with exogenous
oversampling represents a definite improvement over conventional PEA. One way to see this
is that conventional PEA seems to more seriously miss identify k than PEA with exogenous
oversampling. This can be seen most clearly in the graphs of the approximate Lagrange
multipliers, computed using (4.9).43
5 .3 .
Approximate M odel Implied M oments
We now examine several model moments computed using our four approximate policy func
tions and the DP algorithm. In table 2 moments related to the unconstrained model are
displayed. The m om ents we computed for this case are as follows: R £ (the mean value of
f ' { k t , 0 t ) + (1 — <$)), R f (the mean return on a one-period-ahead state-uncontingent consump
tion loan), R e - R l (the mean equity premium), a j , j = y , c , i (the standard deviation of
gross output, consumption and gross investment, respectively), p ( y , j ) , j = c , i (correlation
of gross output with consumption and gross investment , respectively), and freq(i < 0) (the
frequency of times that gross investment is negative). The rate of return variables, R s and
R e,
are expressed in annual percent terms. In table 3 moments related to the constrained
model are displayed. In addition to the moments displayed in table 2 we com pute moments
related to Tobin’s q, the price of new capital in terms of consumption goods.
We define
this price as follows: q = 1 - X / U ' { c ) . 44 Also, we replace freq(? < 0) with freq(A > 0) (the
frequency of times that the gross investment constraint binds) in table 3. All statistics are
based on averages from samples of length 114 replicated 500 times. Numbers in parentheses
are Monte Carlo standard errors.45
43W ith N p = 3, calculations using conventional PEA for T =10,000, 20,000, 40,000, 60,000, and 80,000,
resulted in the following estim ates of k : 33.03, 33.15, 33.27, 33.31, 33.34. W ith N p = 3* (i.e., PEA with
exogenous oversam pling) and T = 10,000, we obtained k = 33.29. W ith N p = 8** (i.e., PEA -collocation)
we obtained k = 33.37. We take the_fc implied by DP, which is k = 33.40, to be the exact solution. T hus, a
g i\e n le\el of accuracy (in term s of k) can be achieved w ith a lower value of T by applying alternatives to
conventional PEA .
44See Sargent (1980) for an analysis of T obin’s q in a setting sim ilar to ours.
4oT hese are <rx / \/500, w here <rx is the stan d ard deviation, across our 500 replications, of some statistic, x.
33
Consider table 2. W ith conventional PEA there is slight inaccuracy in its predictions
for financial variables ( e . g . the equity premium is 0.089 % (0.011) with conventional PEA
versus 0.049 % (0.010) with DP, with standard errors in parentheses) and the algorithm over
estim ates the frequency that gross investm ent is negative (9.89 % (0.38) with conventional
PEA and 8.89 % (0.36) with DP.) PEA with exogenous oversampling and PEA -collocation
both show improvement along these dimensions. Spectral-Galerkin shows convergence at
N J =
5, FEM collocation shows convergence at N * = 72, and FEM Galerkin at N * = 36.
Notice FEM collocation with N * = 36 is unable to achieve convergence while FEM Galerkin,
using exactly the same approximating function, does achieve convergence. This reflects the
relative com putational efficiency of sm oothly weighting the Euler residuals, which Galerkin
does, in the context of the FEM.
Now consider table 3. The primary differences between m ethods, which in any case are
small, lie in their implications for statistics involving financial variables. For example, the
main difference between conventional PEA and the other versions of PEA is that the former
over predicts R e, R e — R f , and a q.
6. C o n c l u d i n g R e m a r k s
Our purpose in this paper is to provide researchers working with more complex model
economies than the one studied here, with some guidance to help select from among the
many available solution algorithms. We expect that in these more complex problems, compu
tational speed and programming convenience will be im portant, desirable characteristics, in
addition to accuracy. W ith this in mind, we compared and evaluated six com putational algo
rithms for solving models with occasionally binding inequality constraints. These algorithms
include: three versions of Marcet’s parameterized expectations algorithm (PEA ); a version
of Judd's Spectral-Galerkin algorithm, extended here to include a Lagrange m ultiplier func
tion as one of the objects sought; two finite elem ent m ethods, Coleman’s FEM -collocation
algorithm, modified to accommodate a Lagrange multiplier, and M cG rattan’s FEM-Galerkin
algorithm, which accom m odates inequality constraints by including a penalty function in the
34
objective. In addition, to provide a benchmark solution, we also did dynamic programming
on a discretized version of our model with a very fine grid. A unique feature of our analysis
is that we illustrate the use of the Euler residual function in evaluating the accuracy of a
solution algorithm.46
To our initial surprise, all the algorithms worked quite well. We were particularly sur
prised at the accuracy with which several of the algorithms predict the Lagrange multiplier.
Even algorithms such as M cG rattan’s and a version of M arcet’s, which com pute the mul
tiplier indirectly, provide reasonable estim ates of this function.
Also, for the most part
all of the algorithms are reasonably accurate for computing particular statistics involving
endogenous variables from the example model economy.
Still, we have developed information we believe is useful for discriminating among these
algorithms. By far the easiest algorithm to implement is M arcet’s PEA. As Marcet (1988)
points out, the algorithm requires virtually no modification to accom m odate inequality con
straints. In the case of the other algorithms, accom m odating inequality constraints involves
substantial complications.
For example, implementation of M cG rattan’s method requires
considerable ‘baby sittin g’ of the computer program, as one tries out various values of a
penalty function parameter. The Spectral-Galerkin and FEM -collocation m ethods also en
tail additional complications to accom m odate inequality constraints. This reflects the fact
that they require directly parameterizing a Lagrange multiplier function, in addition to the
policy functions.
While Marcet’s PEA seems to be the easiest to im plem ent, we had difficulties with
conventional versions of it. A key component of those versions is a cumbersome nonlinear
regression step, potentially involving tens of thousands of observations. One reason for the
large number of observations is that the explanatory variables are inefficiently concentrated
in a narrow range. We devised an alternative (PEA -collocation), in which the regression step
is linear, the explanatory variables are orthogonal, and the required number of observations
4bFor an altern ativ e procedure for evaluating the accuracy of a solution algorithm , see den Haan and
M arcet (1994).
35
in the regression is very small: no more than sixteen in our experim ents. This m ethod
produced results as accurate as the best other m ethod, and is orders of m agnitude faster.
Although it is clear that PEA-collocation is the best solution m ethod for our exam ple,
that does not guarantee that it will dom inate in higher dimensional cases. Here, there are
at least two considerations. First, do the linearity and orthogonality properties of PEAcollocation survive into multidimensional settings? In appendix 2, we define m ultidim en
sional PEA-collocation and show that these properties do indeed survive in general. The
second consideration involves the mapping from a parameterized expectation function to
policy and multiplier functions, which is at the heart of any PEA. In our exam ple, this map
ping is trivial, but in higher dimensions it involves solving nonlinear equations. In principle,
there could be examples in which this is very costly in programmer and/or computer time,
in which case perhaps an alternative method might dominate. Here, it should be born in
mind that P E A ’s have been applied routinely in high dimensional models (see footnote 5).
36
References
[1] Aiyagari, S. Rao. 1993. Frictions in Asset Pricing and Macroeconomics. Federal Reserve
Bank of Minneapolis Working Paper No. 518, September.
[2] Aiyagari, S. Rao and Mark Gertler. 1991. Asset Returns W ith Transactions Costs and
Uninsured Individual Risk. J o u r n a l o f M o n e t a r y E c o n o m i c s , 27: 311-331.
[3] Atkeson, Andy and Patrick Kehoe. 1993. Industry Evolution and Transition: The Role
of Information Capital. Manuscript.
[4] Bizer, David S., and Kenneth L. Judd, 1989. Uncertainty and Taxation. A m e r i c a n
E c o n o m i c R e v i e w , P a p e r s a n d P r o c e e d i n g s , May.
[5] Boldrin, Michele, Lawrence J. Christiano, and Jonas D.M. Fisher. 1994. Habit Persis
tence and the Equity Premium in the Neoclassical Growth Model. Manuscript.
[6] Chari, V .V ., Lawrence J. Christiano, and Patrick Kehoe. Forthcoming. “Policy Analysis
in Business Cycle Models'’, in Cooley, editor, F r o n t i e r s o f B u s i n e s s C y c l e R e s e a r c h ,
Princeton University Press, Princeton, NJ.
[7] Christiano, Lawrence J. 1991. Modeling the Liquidity Effect of a M onetary Shock. Fed
eral Reserve Bank of Minneapolis Q u a r t e r l y R e v i e w , Winter.
[S] Christiano, Lawrence J. and Terry Fitzgerald. 1991. The M agnitude of the Speculative
Motive for Holding Inventories in a Business Cycle Model. Federal Reserve Bank of
Minneapolis Working Paper 415.
[9] Coleman, Wilbur John. 1988. Money, Interest and Capital in a Cash-in-Advance Econ
omy. International Finance Discussion Paper No. 323. Washington: Board of Governors,
Federal Reserve System.
[10] Coleman, Wilbur John, Christian Gilles, and Pamela Labadie. 1992. The Liquidity
Premium in Average Interest Rates. J o u r n a l o f M o n e t a r y E c o n o m i c s , 30: 449-465.
[11] Danthine, Jean-Pierre and John B. Donaldson. 1981. Stochastic Properties of Fast vs.
Slow Growing Economies. E c o n o m e t r i c a 49: 1007-1033.
[12] den Haan, Wouter. 1993. Heterogeneity, Aggregate Uncertainty and the Short Term
Interest Rate: A Case Study of Two Solution Techniques. University of California at
San Diego Manuscript.
[13] den Haan, Wouter and Albert Marcet. 1990. Solving the Stochastic Growth Model by
Parameterizing Expectations. J o u r n a l o f B u s i n e s s a n d E c o n o m i c S t a t i s t i c s 8: 31-34.
[14] den Haan, Wouter and Albert Marcet. 1994. Accuracy in Simulation. R e v i e w o f E c o
n o m i c S t u d i e s . 15*
[15] Heaton. John and Deborah Lucas. 1992. The Effects of Incomplete Insurance Markets
and Trading Costs in a Consumption-Based Asset Pricing Model. J o u r n a l o f E c o n o m i c
D y n a m i c s a n d C o n t r o l . 16: 601-620.
37
[16] Huggett, Mark. 1993. The Risk-Free Rate in Heterogeneous-Agent Incomplete-Insurance
Economies. J o u r n a l o f E c o n o m i c D y n a m i c s a n d C o n t r o l , 17: 953-969.
[17] Judd, Ken. 1991. Minimum W eighted Residuals M ethods for Solving Aggregate Growth
Models. Institute for Empirical Macroeconomics Discussion Paper No. 49.
[18] Judd, Ken. 1992a. Projection M ethods for Solving Aggregate Growth Models. J o u r n a l
o f E c o n o m i c T h e o r y , 58: 410-452.
[19] Judd, Ken. 1992b. N u m e r i c a l M e t h o d s i n E c o n o m i c s . Manuscript.
[20] Judd, Ken. 1993. Comments on Marcet, Rust and Pakes. Forthcoming comments made
at the World Congress of the Econometric Society, Barcelona.
[21] Judd, Iven and Antonio Bernardo. 1994. Projection M ethods for Computing Conditional
Expectations. Stanford University Working Paper.
[22] Kahn, James. 1992. Why is Production More Volatile than Sales? Theory and Evi
dence on the Stock-out Avoidance M otive for Inventory Holding. Q u a r t e r l y J o u r n a l o f
E c o n o m i c s , 107: 481-510.
[23] Kivotaki, Nobuhiro and John Moore. 1993. Credit Cycles. Manuscript.
[24] Luenberger, David. 1969. O p t i m i z a t i o n b y V e c t o r S p a c e M e t h o d s . New York: John W ilev
& Sons.
[25] M aicet, Albert. 1988. Solving Nonlinear Stochastic Growth Models bv Parametrizing
Expectations. Carnegie-Mellon University, mimeo.
[26] Marcet, Albert and Juan Ketterer. 1989. Introducing Derivative Securities: A General
Equilibrium Approach. Carnegie-Mellon University Working Paper.
[27] Marcet, Albert and Ramone Marimon. 1992. Com m unication, Com m ittm ent and
Growth. J o u r n a l o f E c o n o m i c T h e o r y , 58: 219-249.
[28] Marcet, Albert and David Marshall. 1994. Solving Non-linear Rational Expectations
Models by Parameterizing Expectations: Convergence to Stationary Solutions. North
western University Working Paper.
[29] Marcet, Albert and Ken Singleton. 1990. Equilibrium Asset Prices and Savings of H et
erogeneous Agents In the Presence of Portfolio Constraints. Unpublished manuscript.
[30] Marshall, David. 1992. Inflation and Asset Returns in a Monetary Economy. J o u r n a l o f
F i n a n c e , 47: 1315-42.
[31] McCurdy, Tom and Nicholas Ricketts. 1992. An international econom y with countryspecific money and productivity growth processes. Queen’s University Manuscript.
[32] McGrattan, Ellen. 1993. Solving the Stochastic Growth Model with a Finite Element
Method. Federal Reserve Bank of Minneapolis mimeograph.
[33] Press, W illiam H., Saul Teukolskv, W illiam Vetterling, and Brian Flannery. 1992. N u
m e r i c a l R e c i p e s in F o r t r a n . New York: University of Cambridge Press. 3
4
[34] Reddy, J.N. 1993. A n I n t r o d u c t i o n t o t h e F i n i t e E l e m e n t M e t h o d . New York: McGrawHill.
38
[35] Sargent, Thomas. 1980. Tobin’s q and the Rate of Investment in General Equilibrium.
Carnegie-Rochester Conference Series on Public Policy 12.
[36] Stokey, Nancy and Robert E. Lucas, Jr., with Edward Prescott. 1989. R e c u r s i v e M e t h o d s
i n E c o n o m i c D y n a m i c s . Cambridge, MA: Harvard University Press.
[37] Stroud, A.H. 1971. A p p r o x i m a t e C a l c u l a t i o n o f M u l t i p l e I n t e g r a l s . Englewood Cliffs, NJ:
Prentice-Hall.
[38] Telmer, Chris. Forthcoming. Asset Pricing Puzzles and Incomplete Markets. J o u r n a l o f
Finance.
39
Appendix 1: The Dynamic Programming Algorithm
Our DP algorithm is standard. It involves first iterating to convergence on a value
function and then deriving a decision rule from the converged value function. The mapping
that we iterated on is:
uj+ i(*> °)
{«(*> k'i°) +
=
<r) +
-< r )]j,
for 0 £ 0 and A' £ A = {Ax, A2, . . . , A*a/} . Also,
u(A', k , 0 ) = ln[exp(0)A° + (1 — 6)A — A']
and
A(A, 0) = A 0 {A' : ( 1 - 6 ) k < A' < exp { 9 ) k a + ( 1 - 6 ) k }
for the constrained problem, and
A ( k , 0) = k r \ { k ‘ :
0 < k ‘ < exp ( 6 ) k a + (1 - 6 ) k }
for the unconstrained problem. Here, i’j(-,<r) and v j ( - , — a ) are points in 3fJA/, j = 1 , 2 , __
Also v 0 ( k \ 0 ) = 0, for 0 E 0 and k' G k. The points in k are equally spaced with A*,- < A-1+1, i =
1 , 2 , . . . . M — 1, A'i = 16.9, k t\ f = 55.1, and A I = 20,000. We iterated on the above mapping
until reaching a fixed point which was assumed to be achieved when | (e_, — c 7_ 1)./T J_ 1 | <
1 x 10_ ‘, here | x \ is the largest element of x in absolute value and x . / y represents element
by element division of the vectors x and y . Denote the fixed point by v . We then computed
the two decision rule vectors G'(-,cr), G ( - , —a ) £
as follows.
G(k,0) =
argmax | t i ( k , k \ 6 ) +
k<dA(k,9) t
l
[t>(A*', <r) + v ( k \ - c r )]l ,
J
where 0 £ 0 and A' £ A.
The DP investment decision rules graphed in section 5 and the DP second moment
properties are based on G ( k , 9 ) . The DP version of the multiplier reported in section 5 is
computed as follows.
A(Ai,0)
u 1( k i , G ( k i , 0 ) , 0 ) - v l ( k i , 0 )
1-6
,i = 1,2,..., M.
Here, Hi is the derivative of u with respect to its first argument. Also V\ is our estim ate
of the derivative of v with respect to its first argument. We obtained this estim ate by first
fitting, by least squares, a seventh order polynomial to c(A,, 9 ) , i = 1 ,2 , . . . , M for 9 £ 0:
c(A„ 9) = /3O(0) + /91(0)S5(Ai) + •••+ /?7(0)b(A,)]7,i = 1,2,..., M.
40
Here <p :[&i, 1z m \ — * [0,1]. Then,
v 1( k i , d ) = 0 , ( 0 )
+ 2p 2{ e ) v { k i ) + • • • + 7 0 T( 0 ) [ (p ( k i )}6 , t = 1 , 2 , . . . , M .
Appendix 2: Multidimensional Applications of PEA-Collocation
In this appendix we describe how PEA-collocation is applied in models with an arbitrary
finite number of endogenous and exogenous state variables. We show that the principle
qualitative features of PEA-collocation (e.g., linearity of the regression and orthogonality of
the regressors) survive in a multidimensional setting. We also show that PEA-collocation
encounters a ‘curse of dim ensionality’ problem in very high dimensional system s. We propose
an alternative, PEA-Galerkin, for dealing with circumstances like this.
The P ro b le m
Let k £ A C 3?' denote a vector of endogenous state variables. We suppose that the
exogenous variables, 0 6 0 C 3?m, are a first order, stationary Markov process with transition
density p e, ( 0 ' | 0 ) . Since we restrict m only to be finite, this is equivalent to assuming the
exogenous variables have an arbitrary finite ordered Markov representation. In contrast to
the analysis in the main text, here we suppose 0 is a vector of continuous random variables.
Let u : A x A x 0 —►3? denote the one-period return function. This may be an indirect utility
function that iesults after static decisions, such as labor supply in standard business cycle
models oi the sectoial allocation of capital in a multisector model, have been maximized
out. We consider model economies which lead to the analysis of the following functional
equation:11'
W(M)
max
u ( k , k r, 6 ) + 0 f
W ( k ' , 9 ' ) p e, ( 0 ' | 0 ) d 9 \ for all k 6 K , 0 £ 0
k'er{k,0),G{k,k',6)>ouJo's®
where G is a w x 1 vector-valued function, G : I \ x A' x 0 -► 3^’. Also, IT : A' x 0 —►3? is a
value function, and T : A x 0 —> A and G ( k , k ' , 9 ) > ()„, characterize the constraints on the
choice of k in the maximization. Here, 0U, denotes a w x 1 vector of zeros. The correspondence
T summaiizes constraints that either bind always or never, while the function G summarizes
restrictions that bind occasionally.*48 In what follows it is convenient to use the notation
•s = vec ( k , 6 ) , where s is a q x 1 vector, q = / + m.
Let g : A x 0 -» I \ be the (single-valued) policy function which attains the maximum in
the functional equation, and let h : I \ x 0 - > $ t w denote the i v x 1 vector-valued multiplier
4< For a detailed discussion of m odel economies like this, see Stokey and Lucas w ith P resco tt (1989).
48By a constraint never binding, we m ean th a t its m ultiplier is zero for alm ost all (k,0) £ I\ x 0 .
An exam ple of this is the one sector grow th m odel w here T sum m arizes th e nonnegativity co nstraints on
consum ption and capital, and the ap p ro p riate In ad a conditions on u are satisfied. By a constraint always
binding, we m ean th a t its m u ltiplier is alm ost always nonzero. An exam ple of this is the resource constraint
in the m odel ju st described. In the m odel analyzed in the body of the p aper, T sum m arizes the resource
constraint and nonnegativity co n strain ts on capital and consum ption, and th e co n straint, G(k,k',8) > Ou,,
is the nonnegativity constraint on gross investm ent.
41
f u n c t i o n c o r r e s p o n d i n g t o t h e c o n s t r a i n t s , G ( f c , k',6) > 0 ^ . W e s u p p o s e t h a t t h e p o l i c y a n d
m u lt ip lie r fu n c t io n s m u s t s a t is f y t h e E u le r e q u a t io n s ,
u2 { k, g{k, 0), e) + G 2 ( k , g ( k , 0 ) , e ) T h{k,e) + l3 j m{k, 9, ff-, g, h)pe,{ff \G)dff = %
J 81
(6.1)
and the Kuhn-Tucker conditions,
G( k , g ( k , 0 ) , 6 ) > 0w, h(k, 0) > On,, and h(k,6). * G { k , g ( k , 0 ) , 0 ) = 0W,
(6.2)
where .* denotes element-by-element multiplication. Also,
m(k, M 'l g, h) = «, ( g ( k , 0), g ( s ( M ) , «') , $') + G, (g(k, 6 ) , g ( g ( k , « ),« '), 0')T h (g(k, 6), O' ) ,
(6.3)
where m : I\ x 0 x 0 —*• 3^, is an l x 1 vector-valued function. In (6.1)-(6.3), tq denotes the
/ x 1 vector of derivatives of u with respect to the ith argument, G{ denotes the iv x / matrix
of derivatives of G with respect to its ith argument, and T denotes the matrix transposition
operator. The problem is to approximate g and h, solutions to the functional equations
(6.1)-(6.3).
The PEA
The PEA approximates g and h indirectly by parameterizing the j-th conditional expec
tation in (6.1) by a function exp(eaJ(s)), j — 1 , 2 , . . . , / :
exp(eaJ(s)) ~ [
m 3(s,6';g,h)pe>{0' \ 0)d$’, j = 1 , 2 , . . . , / ,
Je'eo
(6.4)
for all
A'x 0 . Here, a3 G
is a finite set of parameters, and m 3 is the j th element of the
function, m. Let a = vec(«1, a2, . . . , a1). We define a mapping from a to policy and multiplier
functions, ga and ha, as follows. For any given «, replace the conditional expectation in (6.1)
by the parameterized expectation in (6.4) and let ga and ha denote the policy and multiplier
functions which satisfy the Euler equation and Kuhn-Tucker conditions.49 These policy and
multiplier functions then imply a conditional expectation function, f g,e& m(s, O'; ga, h a )pe>{0' \
0)d0' for all s £ A’ x 0 . A new vector of parameters, a' = S( a]l ■N p), is then chosen to
make the function exp(ea/(s)) as close as possible to this conditional expectation. Here, INP
denotes the number of elements in a. A PEA seeks an a* such that a* = 5(a*; / • N p).
PEA-Collocation
As in the text, we construct ea(s) using Chebyshev polynomials. We begin by defining
what these are in a multidimensional setting, and by displaying their discrete orthogonality
property. As in the text, it is this property which accounts for the orthogonal regressor
property of PEA-collocation.
49In general, finding ga and h a for any (A, 9) (E A’ x 0 requires solving a system of equations using num erical
m ethods. In co n trast, in the m odel econom y of the m ain tex t, th e solution to these equations has a simple
analytic form (see (4.8) and (4.9)).
42
Chebyshev polynomials with q-dimensional domain are constructed from the one dimen
sional Chebyshev polynomials studied in the text. Let
= {T0{x), Ti(x) , . . . , Tn- i ( x ) }
denote the basis functions for one dimensional Chebyshev polynomials of degree n — 1, for
n > 1, where x E (—1,1). The tensor product basis for degree n — 1 Chebyshev polyno
mial functions of q variables is constructed by taking all possible <?-term products of the n
elements in
Accordingly, the resulting basis is
^
=
{ T i 1( x i ) T i 2 { x 2) - - - T iq( x q)
| ij
= 0 , 1, . . . , n
-
1,j
= 1 , 2 , . . . , 9}.
Here, (xj, ar2, ..., x 9) is an element of the q-(o\d Cartesian product of (—1,1), which we denote
by (—1, l ) 9. Notice that
contains nq elements.
A convenient feature of this tensor product basis is that it inherits the discrete orthog
onality properties of
(see Judd (1992b, chapter 5) and the references given there). Let
<I>1 , 4>2 , ■• •, <£n« be a list of the elements of
where <f>v : ( —1, l ) 9 —> ( —1,1) for v = 1,..., n9.
Then the discrete orthogonality property in the multidimensional setting is, for i , j < nq,
n<>
0, for i ± j
Ci(n,q), for i = j
<Pi { r v ) <T>j ( r v)
v =l
’
where C{(n,q) are constants that depend on the basis and rv E ( —1,1)7 is composed of a
selection of q elements from the set of n zeros of Tn, v = 1 , 2 , . . . , nq. In particular, the set of
rVs is defined by the nq ways of choosing q of the n zeros of Tn. The zeros of Tn are given by
rk =
/ ( 2 A'—1 )tt\ ,
COS I ---- —---- 1 , k = 1 ,...,
n.
Also, for / = 1 , 2 , . . . , nq
ct{n ,q ) = t
$ > { r v ) 2-
V—\
We construct the parameterized expectation function using the elements of
as follows:
nq
e«u(s) =
^2ai^i((p(s)),j
i=
= 1, 2 ,...,/.
1
Notice that l \ p — nq here. The function is the multidimensional version of the analogous
function used in the main text (see (3.6).) That is <p : n L iU ^ s ;) C J?9 —> (—1, l ) 9, where
(s,-, Si), i = 1 , 2 , . . . , q bound the exogenous and endogenous variables.
We now derive the multidimensional version of the orthogonal regressor result, (3.15).
With PEA-collocation, a' = S( a ; / • nq) is defined by:
R j ( s v -,ga, X > )
- exp(eay(v>(s,,))) - /
m j (sv, $'■, ga, ha)pe>{d' \ 9)dff = 0, j =
J 6*G©
1
9
/
(6.5)
43
for v = 1 , . . . , n9, where s v = <^-1 (rv), v = 1 , . . . , n9 .50 Now (6.5) holds if and only if it holds
for the log of the terms on each side of the minus sign. That is, for each v = 1 , 2 , . . . , n9 and
each j = 1 , 2 , . . . , / , (6.5) holds if and only if
Multiply both sides of each equation in (6.6) by <}>i(sv) and for fixed j sum over v =
1 , 2 , . . . , nq. By the discrete orthogonality property, all terms on the left side of the equality,
except those involving 4>i(<p(sv))2, v = 1, . . . , n 9, are zero. Repeating this procedure for
4>2,4>3,...,(f>n<,i one finds that, analogous to (3.15), the mapping a' = S(a-,l ■nq) has the
simple analytical form,
PEA-Galerkin
A disadvantage of a tensor product basis is that the number of elements in the basis
grows exponentially as the dimension increases. One could instead work with a strict subset
of the number of elements in the tensor product basis. For example, Judd (1992b, chapter
5) suggests working with the following subset:
Notice that C ^ C $J,9\ since
simply deletes high-order cross product terms in
For
very large problems the computational burden of finding the a* that solves a* = S( a*; / • n9)
may be unacceptable. In these circumstances, a useful alternative may be to use a smaller
basis. However, if one continues to work with the Chebyshev zeros, rv, v = 1 , 2 , . . . , n9, then
PEA-collocation is no longer implementable in general. This is because PEA-collocation now
attempts to solve the n9 equations, (6.6), using less than nq unknowns. There are several
options. One is to apply PEA-collocation to a reduced number of equations. Another is to
maintain the number of equations and apply a different weighted residual method.
From the text, it is clear that there are many such methods. One such method is a
modified version of the Galerkin procedure discussed in the text, PEA-Galerkin, which is
applied as follows. First, select a value for n, and choose a subset of Arp basis functions from
$Jj9\ where N p < nq. Then, compute rv,u = 1,2, . . . , n 9 as above and form the N p x nq
50T he integral in (6.5) could be ap p roxim ated by q u ad ratu re or M onte C arlo m ethods.
44
m a t r i x A a n d t h e nq x 1 v e c t o r FLj(s,ga,ha) a s f o l l o w s :
A =
<Mrt)
M pi)
^i(r2)
<f>2(r2)
•
• ■■
<f>2(rn«)
i Rj(S') (Jai
. 4>Np{rx) ^Arp(fi*2) •
—
1?Qai ha)
2i Qai ha)
Qa') ha)
where s = [$i,. . . , s n9], and s v, v = 1 , 2 , . . . , n9, is as defined above. By the discrete orthog
onality property, the rows of A are orthogonal. Finally, find the value of a that solves the
system of l • N p nonlinear equations:
ARj { s , g a, h a) = 0 ,j = 1 , 2 , . . . , / .
45
T a b le 1
C o m p u ta tio n T im e s a n d M a x im u m A b s o lu te V a lu e s
o f t h e E u le r R e s id u a ls fo r t h e V a r io u s A lg o r it h m s
FEM-Collocation:
FEM-Galerkin:
Constrained
PEA:
N p = 3*
151.9
N p = 3**
0 .6
NJ = 3
NJ = 5
N f = 36
N* = 72
N* = 18
N* = 36
3.5
6.3
253.1
557.6
Np
=
3
17.9
174.4
181.7
b*
/o?
N p = 3*
8 .8
MAV of Euler Residuals
0 = —cr
6= a
9.2 x 10~ 5
1.3 x 10~ 4
( 1 . 0 x 1 0 "4) (1.6 x lO"3)
3.3 x lO" 5
1 .6 x 1 0 ~ 4
(5.1 x 10~4) (1.5 x lO"3)
1.7 x 10~ 6
2 . 1 x lO" 6
(3 . 2 x 1 0 ~6) (4.1 x 10~6)
4.5 x lO” 4
1 .2 x 1 0 ~ 4
3.7 x lO" 7
6 .2 x 1 0 - 7
1 .0 x 1 0 ~ 4
7.8 x 10~ 5
3.3 x 10" 5
2.7 x 10" 5
3.7 x 10" 4
1.3 x lO" 4
7.2 x lO" 5
4.8 x 10“ 5
1 .1 x 1 0 ~ 4
1.7 x 10" 4
(1.6 x lO"4) (8 . 8 x lO"4)
3.3 x lO" 5
1.9 x 10" 4
(7.2 x 10"3)
2 .6 x 1 0 " 5
9.9 x 10" 5
(2 . 6 x 1 0 ~5) (9.9 x 10-5)
6 .2 x 1 0 - 4
7.3 x 10" 4
8 .8 x 1 0 - 5
8 .0 x 1 0 - 5
1 .0 x 1 0 “ 4
7.8 x 10“ 5
3.4 x 10" 5
2 .8 x 1 0 " 5
2.4 x lO" 2
7.3 x 10- 3
3.7 x 10- 4
1.5 x lO" 4
Np =
Spectral-Galerkin:
FEM-Collocation:
FEM-Galerkin:
8 **
(5,3,2)
(9,5,4)
N f = 36
N* = 72
N f = 18
N f = 36
2.5
17.6
27.9
522.9
996.2
212.5
1783.9
1
Spectral-Galerkin:
Time
(Seconds)
155.1
•fe.
Approximation Method
Np = 3
PEA:
X
o
Model
Unconstrained
Notes: (i) For the S pectral-G alerk in and FE M -G alerkin approxim ations to the constrained m odel we used
the u nconstrained approxim ations for sta rtin g values. S tartin g values for the unconstrained calculations were
based on the log-linear ap p ro x im ate decision rules. T he tim es displayed include co m p u tatio n tim es for these
sta rtin g values, (ii) In practice, the FE M -G alerkin algorithm as applied to the constrained m odel involves
solving p enalty function versions of the m odel for several values of the penalty p aram eter 7r. T he com putation
tim es for this ap p licatio n reflect this fact, (iii) T he ordered trip lets for the S pectral-G alerkin approxim ations
to the constrained m odel correspond to (N(cr), N ( — a), N x ). (iv) See the te x t for an explanation of the
asterisks associated w ith the PE A entries, (v) T he MAV num bers corresponding to the PEA have the
following in te rp re ta tio n . T hey are based on confidence intervals from sim ulations of length 10,000 based
on im plied decision rules for the approxim ation in question. T he num bers not in parenthesis correspond
to MAVs based on an interval containing 90% of the realizations of the capital stock. T he num bers in
parenthesis correspond to MAVs based on an interval th a t contains 100% of the realizations.
T a b le 2
S t a t is t ic s fr o m V a r io u s A p p r o x im a t io n s o f th e U n c o n s t r a in e d M o d e l
Statistic
Re
Rf
Re - Rf
ay
p(y,c)
p(y,
*)
Freq(z < 0)
PEA: N p =
3*
3
DP
3.092
3.099
3.120
(0.007) (0.006) (0.006)
3.052
3.040
3.031
(0.015) (0.015) (0.015)
0.047
0.089
0.051
(
0
.
0
1
1
)
(
0 .0 1 1 )
(0 .0 1 0 )
62.2
62.2
62.2
(0.04)
(0.04)
(0.04)
7.21
7.22
7.25
(0.09)
(0.09)
(0.09)
59.9
60.0
60.1
(0.04)
(0.04)
(0.04)
0.354
0.362
0.361
(0 .0 0 2 ) (0 .0 0 1 ) (0 .0 0 1 )
0.993
0.993
0.993
(0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 )
9.89
9.17
8.89
(0.37)
(0.36)
(0.38)
FEM-G: N * =
FEM-C: AT/ =
72
3**
18
36
5
36
3.102
3.090
3.088
3.092
3.090
3.073
(0.007)
(0.007) (0.007) (0.007) (0.007) (0.007)
3.042
3.043
3.043
3.039
3.043
3.049
(0.015)
(0.015) (0.015) (0.015) (0.015)
(0.015)
0.047
0 .0 1 2
0.047
0.024
0.045
0.063
0.049
(0 .0 1 0 ) (0 .0 1 1 ) (0 .0 1 0 ) (0 .0 1 0 ) (0 .0 1 0 ) (0 .0 1 0 ) (0 .0 1 0 )
62.2
62.1
62.2
62.2
62.2
62.2
62.3
(0.04)
(0.04)
(0.04)
(0.04)
(0.04)
(0.04)
(0.04)
7.13
7.31
7.13
7.13
7.12
7.15
7.13
(0.09)
(0.09)
(0.09)
(0.09)
(0.09)
(0.09)
(0.09)
60.1
60.1
60.1
59.8
60.3
60.0
60.1
(0.04)
(0.04)
(0.04)
(0.04)
(0.04)
(0.04)
(0.04)
0.358
0.370
0.358
0.357
0.358
0.358
0.358
(0 .0 0 1 ) (0 .0 0 1 ) (0 .0 0 1 ) (0 .0 0 1 ) (0 .0 0 1 ) (0 .0 0 1 ) (0 .0 0 1 )
0.993
0.993
0.993
0.993
0.993
0.993
0.993
(0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 )
8.64
8.72
10.50
8.23
8.59
8.97
8.67
(0.37)
(0.39)
(0.37)
(0.35)
(0.36)
(0.37)
(0.36)
S-G:
3
3.146
(0.006)
3.022
(0.016)
N J =
Notes: (i) Statistics shown are averages from samples of length 114 replicated 500 times, (ii) Freq(t < 0) indicates the per cent rate at which gross investment is
negative across samples, (iii) Numbers in parenthesis are Monte Carlo standard errors, (iv) See the text for a description of the asterisk notation used for the
PEA entries in the table, (v) Finally, S-G stands for Spectral-Galerkin, FEM-C stands for FEM-Collocation and FEM-G stands for FEM-Galerkin.
T a b le 3
S t a t is t ic s fr o m V a r io u s A p p r o x im a t io n s o f t h e C o n s t r a in e d M o d e l
Statistic
Re
Rl
Re - Rf
°y
<7c
0i
0q
p(y,c)
p(y,i)
p{ y, <i )
Freq(A > 0)
PEA: N p =
r
3
DP
3.126
3.160
3.126
(0.005) (0.005) (0.005)
3.052
3.043
3.055
(0.017) (0.018) (0.017)
0.074
0.071
0.117
(0.015) (0.016) (0.015)
62.2
62.2
62.2
(0.04)
(0.04)
(0.04)
7.20
7.08
7.15
(0.08)
(0.08)
(0.09)
59.7
59.7
59.6
(0.03)
(0.03)
(0.03)
0.319
0.341
0.318
(0.014) (0.014) (0.014)
0.404
0.411
0.400
(0 .0 0 2 ) (0 .0 0 2 ) (0 .0 0 2 )
0.994
0.994
0.994
(0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 )
0.236
0.248
0.229
(0.008) (0.008) (0.008)
10.92
9.67
1 0 .1 1
(0.41)
(0.40)
(0.39)
Spectral- Galerkin
FEM-G: N* =
FEM-C!: N f =
(5,3,2)
(9,5,4)
72
36
18
36
3.122
3.121
3.125
3.588
3.125
3.128
3.109
(0.005) (0.006) (0.006) (0.005) (0.006) (0.017) (0.006)
3.053
3.054
3.053
3.051
3.058
3.069
3.053
(0.017) (0.017) (0.017) (0.017) (0.017) (0.027) (0.017)
0.072
0.072
0.076
0.051
0.068
0.519
0.070
(0.015) (0.015) (0.015) (0.013) (0.015) (0.043) (0.015)
62.2
62.2
62.2
62.2
62.2
62.3
61.9
(0.04)
(0.04)
(0.04)
(0.04)
(0.04)
(0.04)
(0.04)
7.09
7.09
7.08
7.10
7.07
8.96
7.07
(
0
.
1
1
)
(0.08)
(0.08)
(0.08)
(0.09)
(0.08)
(0.08)
59.7
59.7
59.8
59.8
59.8
58.3
59.8
(0.03)
(0.03)
(0.03)
(0.03)
(0.03)
(0.05)
(0.03)
0.312
0.302
0.283
0.310
0.295
1.070
0.303
(0.014) (0.014) (0.014) (0.014) (0.014) (0.040) (0.014)
0.400
0.409
0.399
0.396
0.398
0.445
0.398
(0 .0 0 2 ) (0 .0 0 2 ) (0 .0 0 2 ) (0 .0 0 2 ) (0 .0 0 2 ) (0.004) (0 .0 0 2 )
0.994
0.994
0.994
0.994
0.994
0.994
0.990
(0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 ) (0 .0 0 0 2 )
0.226
0.217
0.157
0 .2 0 1
0 .2 2 0
0 .220
0.218
(0.008) (0.007) (0.008) (0.007) (0.007) (0.006) (0.007)
9.04
4.53
9.17
9.73
8.08
9.31
9.31
(0 .2 2 )
(0.39)
(0.39)
(0.40)
(0.36)
(0.39)
(0.38)
8 **
Notes: (i) Statistics shown are averages from samples of length 114 replicated 500 times, (ii) Freq(A > 0) indicates the per cent rate at which the constraint
binds across samples, (iii) Numbers in parenthesis are Monte Carlo standard errors, (iv) See the text for a description of the asterisk notation used for the PEA
entries in the table, (v) FEM-C stands for FEM-Collocation and FEM-G stands for FEM-Galerkin. (vi) The ordered triplets is the Spectral-Galerkin columns
correspond to (N(cr), N ( —a), N x).
0
>rNC^m0
C oO Si
n^irN^"
Figure 2a.
Euler
residuals
Reversible Invest ment
Reversible Inves t ment
PEA: 75 = a
PEA: $ = - a
f o r PEA a n d S p e c t r a 1 - G a l e r k i n
Irreversible Inves t ment
PEA:
Irreversible Invest ment
= a
PEA:
= —a
o«-
- 0.0001
0.0001
0.0003
0.0005
to
o
Capital Stock
1
S p e c t r a l - Ga l e r ki n:
1
Capital Stock
Spectral —Galerkin:
= -a
0.0006
Spect ral —Galerkin: 1$ = cr
Capital Stock
= a
Capital Stock
S p e c t r a l —G a le rk in :
oo
-0.0006
- 0.0002
0.0002
o
1
22
26 30 34 38
Capital Stock
42
22
26 30 34 38
Capital Stock
42
Capital S to c k
= —cr
Figure 2b.
Reversible Investment
Euler r e sid u a ls
for
F E M - c o l 1o c a t i o n and F E M - G a l e r k i n
Irreversible Invest ment
Reversible Invest ment
FEM-Crelocation:
FEM—Collocation: 'O = a
= —a
FEM —Collocation: ^
Irreversible Invest ment
FEM—Collocation:
= —a
l
O
C
M
ln
-0.000025 0.000025 0.000075
c\j
0.0004
-0.0002
-0.0000
0.0002
FEM—Galerkin:
=a
FEM-Galerkin: a? = —a
cn
FEM-Galerkin:
$ = a
FEM-Galerkin:
= —a
Figure 3.
=
im p lie d by Dynam ic Program m ing and S p e c t r a l - G a l e r k i n
Unconstrained Investment Rules:
a
=a
Constrained Investment Rules: $ =
o
1.14
Investment Rule:
Policy functions
Capital Stock
Unconstrained Investment Rules:
lO
Capital
= -a
St oc k
Constrained Investment Rules: # = -a
lO
- 0 . 1 5
-0 .0 5
0.05
0 .15
0.25
Capital Stock
Investment Rule:
= -cr
0.03
A Rule:
1
1
1
1
7?
1
= -a
p
A Rule:
ro
0.02
--- Unconstrained
-- Constrained
-0 .0 1
0.00
0.01
X'
21
23
25
27
29
31
3 3
Capital Stock
35
37
39
41
= —cr
Figure
PEA v e r s u s S p e c t r a l - G a l e r k i n
Investm ent Rules
in t h e R e v e r s i b l e
Un co n st ra in ed Investment Rule
In vestm en t econom y
U n c o n s t ra in e d I nv es tm ent Rule
= —a
1.14
-0.15
1.22
1.30
i9 = <J
21
23
25
27
29
31
33
Capital Stock
35
37
39
41
39
41
39
41
U n c o n s t r a in e d I nv es tm en t Rule
= —a
1.14
-0.15
1.22
1.30
0.05 0.20
-i? = a
21
23
25
27
29
31
33
Capital Stock
35
37
U n c o n s t r a in e d I nv es tm en t Rule
= a
0.05 0.20
1.28
-0.15
1.22
i? = —a
Capital
Stock
21
23
25
27
29
31
33
Capital Stock
35
37
F ig u r e 5-
PEA v e r s u s S p e c t r a l - G a l e r k i n
C o n s t ra in e d I n v e s tm e n t Rule
rd
=
o
I n v e s t m e n t and Lam bda R u l e s
in th e
Co n st ra in ed In ve s tm e n t Rule
45 = - a
Irreversible
I n v e s t m e n t E conom y
C o n st ra in e d X Rule
45 = —a
.14
CN
Capital Stock
Capital Stock
Capital Stock
Co n st ra in ed I n v e s tm e n t Rule
1? = a
Co n st ra in ed I nv es tm en t Rule
45 = —a
C o n s t ra in e d X Rule
CN
45 =
—a
CN
Capital Stock
Capital Stock
C o n s t ra in e d In v e s tm e n t Rule
C on st ra in ed In v e s tm e n t Rule
45 = - a
45 = <7
C o n st ra in e d X Rule
45 = —<7
CM
23 25 27 29 31 33 35 37 39 41
Capital Stock
C apital
Stock
Capital S to ck
W orking P ap er Series
A series ofresearch studies on regional economic issuesrelating tothe Seventh Federal
Reserve District,and on financial and economic topics.
REGIONAL ECONOMIC ISSUES
Estimating Monthly Regional Value Added by Combining Regional Input
With National Production Data
WP-92-8
Philip R. Israilevich and Kenneth N. Kuttner
Local Impact of Foreign Trade Zone
WP-92-9
D a vid D . Weiss
Trends and Prospects for Rural Manufacturing
WP-92-12
William A. Testa
State and Local Government Spending-The Balance
Between Investment and Consumption
WP-92-14
R ichard H. M attoon
Forecasting with Regional Input-Output Tables
WP-92-20
P.R. Israilevich , R. M ahidhara, and G J D . H ewings
A Primer on Global Auto Markets
Paul D . B allew and R obert H. Schnorbus
WP-93-1
Industry Approaches to Environmental Policy
in the Great Lakes Region
W P-93-8
D a vid R. Allardice, R ichard H. M attoon and William A. Testa
The Midwest Stock Price Index-Leading Indicator
of Regional Economic Activity
W P-9 3-9
William A. Strauss
Lean Manufacturing and the Decision to Vertically Integrate
Some Empirical Evidence From the U.S. Automobile Industry
WP-94-1
Thomas H. K lier
Domestic Consumption Patterns and the Midwest Economy
WP-94-4
R obert Schnorbus and Paul B allew
1
Working paperseriescontinued
ISSUES IN FINANCIAL REGULATION
Incentive Conflict in Deposit-Institution Regulation: Evidence from Australia
WP-92-5
E dw ard J. Kane and G eorge G. Kaufman
Capital Adequacy and the Growth of U.S. Banks
WP-92-11
H erbert B aer and John M cElravey
Bank Contagion: Theory and Evidence
WP-92-13
G eorge G. Kaufman
Trading Activity, Progarm Trading and the Volatility of Stock Returns
WP-92-16
Jam es T. M oser
Preferred Sources of Market Discipline: Depositors vs.
Subordinated Debt Holders
WP-92-21
D ouglas D . Evanoff
An Investigation of Returns Conditional
on Trading Performance
WP-92-24
Jam es T. M oser and Jacky C. So
The Effect of Capital on Portfolio Risk at Life Insurance Companies
WP-92-29
Elijah B rew er III, Thomas H. M ondschean, and Philip E. Strahan
A Framework for Estimating the Value and
Interest Rate Risk of Retail Bank Deposits
WP-92-30
D a vid E. Hutchison, G eorge G. Pennacchi
Capital Shocks and Bank Growth-1973 to 1991
WP-92-31
H erbert L. B aer and John N. M cElravey
The Impact of S&L Failures and Regulatory Changes
on the CD Market 1987-1991
WP-92-33
Elijah B rew er and Thomas H. Mondschean
Junk Bond Holdings, Premium Tax Offsets, and Risk
Exposure at Life Insurance Companies
WP-93-3
Elijah B rew er III and Thomas H. Mondschean
2
Working paperseriescontinued
Stock Margins and the Conditional Probability of Price Reversals
W P -93-5
Paul Kqfinan and Jam es T. M oser
Is There Lif(f)e After DTB?
Competitive Aspects of Cross Listed Futures
Contracts on Synchronous Markets
WP-93-11
Paul Kojman, Tony Bouwman and Jam es T. M oser
Opportunity Cost and Prudentiality: A RepresentativeAgent Model of Futures Clearinghouse Behavior
WP-93-18
H erbert L. Baer, Virginia G. France and Jam es T. M oser
The Ownership Structure of Japanese Financial Institutions
WP-93-19
Hesna Genay
Origins of the Modem Exchange Clearinghouse: A History of Early
Clearing and Settlement Methods at Futures Exchanges
WP-94-3
Jam es T. M oser
The Effect of Bank-Held Derivatives on Credit Accessibility
WP-94-5
Elijah B rew er III, Bernadette A. Minton and James T. M oser
MACROECONOMIC ISSUES
An Examination of Change in Energy Dependence and Efficiency
in the Six Largest Energy Using Countries-1970-1988
WP-92-2
Jack L. H ervey
Does the Federal Reserve Affect Asset Prices?
WP-9 2-3
Vefa Tarhan
Investment and Market Imperfections in the U.S. Manufacturing Sector
WP-9 2-4
Paula R. Worthington
Business Cycle Durations and Postwar Stabilization of the U.S. Economy
WP-92-6
M ark W. Watson
A Procedure for Predicting Recessions with Leading Indicators: Econometric Issues
and Recent Performance
WP-92-7
Jam es H. Stock and M ark W. Watson
3
W orking p aper series continued
Production and Inventory Control at the General Motors Corporation
During the 1920s and 1930s
WP-92-10
Anil K . K ashyap and D a v id W. W ilcox
Liquidity Effects, Monetary Policy and the Business Cycle
WP-92-15
Law rence J. Christiano and M artin Eichenbaum
Monetary Policy and External Finance: Interpreting the
Behavior of Financial Flows and Interest Rate Spreads
WP-92-17
Kenneth N. Kuttner
Testing Long Run Neutrality
WP-92-18
R obert G. K ing and M ark W. Watson
A Policymaker’s Guide to Indicators of Economic Activity
WP-92-19
C harles E vans, Steven Strongin, and Francesca Eugeni
Barriers to Trade and Union Wage Dynamics
WP-92-22
Ellen R. Rissman
Wage Growth and Sectoral Shifts: Phillips Curve Redux
WP-92-23
Ellen R. Rissman
Excess Volatility and The Smoothing of Interest Rates:
An Application Using Money Announcements
WP-92-25
Steven Strongin
Market Structure, Technology and the Cyclicality of Output
WP-92-26
Bruce Petersen and Steven Strongin
The Identification of Monetary Policy Disturbances:
Explaining the Liquidity Puzzle
WP-92-27
Steven Strongin
Earnings Losses and Displaced Workers
WP-92-28
Louis S. Jacobson, R obert J. LaLonde, and D aniel G. Sullivan
Some Empirical Evidence of the Effects on Monetary Policy
Shocks on Exchange Rates
WP-92-32
M artin Eichenbaum and Charles Evans
4
W orking paper series continued
An Unobserved-Components Model of
Constant-Inflation Potential Output
W P -03-2
Kenneth N. Kuttner
Investment, Cash Row, and Sunk Costs
WP-93-4
Paula R. Worthington
Lessons from the Japanese Main Bank System
for Financial System Reform in Poland
WP-93-6
Takeo Hoshi, Anil Kashyap, and G ary Loveman
Credit Conditions and the Cyclical Behavior of Inventories
WP-93-7
Anil K. Kashyap, Owen A. Lam ont and Jerem y C. Stein
Labor Productivity During the Great Depression
WP-93-10
M ichael D . B or do and Charles L. Evans
Monetary Policy Shocks and Productivity Measures
in the G-7 Countries
WP-93-12
C harles L. Evans and Fernando Santos
Consumer Confidence and Economic Ructuations
WP-93-13
John G. M atsusaka and Argia M. Sbordone
Vector Autoregressions and Cointegration
WP-93-14
M ark W. Watson
Testing for Cointegration When Some of the
Cointegrating Vectors Are Known
WP-93-15
M ichael T. K. Horvath and M ark W. Watson
Technical Change, Diffusion, and Productivity
WP-93-16
Jeffrey K. Cam pbell
Economic Activity and the Short-Term Credit Markets:
An Analysis of Prices and Quantities
WP-93-17
Benjamin M. Friedman and Kenneth N. Kuttner
Cyclical Productivity in a Model of Labor Hoarding
WP-93-20
A rgia M. Sbordone
5
W orking paper series continued
The Effects of Monetary Policy Shocks: Evidence from the Flow of Funds
W P -94-2
L a w re n c e J . C h ristia n o , M a rtin Eich en bau m a n d C h a rle s E v a n s
Algorithms for Solving Dynamic Models with Occasionally Binding Constraints
WP-94-6
L a w re n ce J . C h ristia n o and Jo n a s D M . F is h e r
6
@FedFRASER