An Algorithm To Solve Dynamic Models

Board of Governors of the Federal Reserve System International Finance Discussion Papers Number 351 May 1989

AN ALGORITHM TO SOLVE DYNAMIC MODELS Wilbur John Coleman II

Note: International Finance Discussion Papers are preliminary materials circulated to stimulate discussion and critical comment. References in publications to International Finance Discussion Papers (other than an acknowledgment that the writer has had access to unpublished material) should be cleared with the author or authors.

ABSTRACT

This paper presents an algorithm to solve recursive systems, formulated in discrete or continuous time, which have an endogenous state variable. The basis of the algorithm is

a fixed point equation in the function from the state variables to the control variables.

An Algorithm to Solve Dynamic Models Wilbur John Coleman II

Consider a system, evolving through time (or space), whose state at time ¢ can be completely described by a finite dimensional vector s t In general, this state consists of variables which evolve exogenously, and those which evolve endogenously, that is, their evolution is determined such that certain equilibrium conditions are satisfied. In addition to this state vector, the system consists of endogenous variables which are not part of the state of the system, but which nevertheless influence the evolution of the endogenous state variables (e.g., control variables). The following pages present an algorithm which may solve this type of system. Three contexts will be developed to describe the algorithm: the stochastic discrete time system, the deterministic continuous time system, and the stochastic continuous time system. Each section will also contain a specific example, and Table 1, which is towards the end of this paper, summarizes the results of applying this algorithm to these examples.

There certainly exist many other methods which have been employed to solve these systems, some of which have clearly been quite successful for particular problems. Value function iteration for discrete time problems, for example, has proven to be quite robust, although its handling of equilibrium conditions (i.e., properties which the control functions must satisfy, but which are not directly part of the maximization problem) is less than

satisfactory. For example, how does one go about solving for an equilibrium price function

*The author is a staff economist in the International Finance Division of the Federal Reserve Board. This paper represents the views of the author and should not be interpreted as reflecting the views of the Board of Governors of the Federal Reserve System or other members of its staff.

(in a dynamic non—Pareto Optimal economy), as there is no natural iterative method which may be used to solve for these functions? The algorithm proposed here essentially integrates the equilibrium conditions with the optimality conditions on the original control functions, thus obtaining a new set of "controls" which has a likewise recursive specification, and thus suggests a natural iterative method for its solution. A benefit of this approach is that often the new set of endogenous functions is defined on a much lower dimensional state space (e.g., in a representative agent economy, individual and aggregate variables are equated in equilibrium). Even without the problem of equilibrium conditions, though, the algorithm seems to compare favorably with other methods (e.g., see the discussion in Taylor and Uhlig, 1989).

The approach in discrete time has apparently been successful for a wide variety of problems (e.g., the multisector optimal growth and monetary growth models, and various other models of asset prices and trade flows), so it seemed natural to attempt to duplicate this success for continuous time models. For these problems the finite difference methods for solving differential equations (e.g., Runge—Kutta) has proven to be quite robust, and in fact it is not clear, for these continuous time setups, if the algorithm proposed here is any better. It is nevertheless a general method to construct control functions (and not simply a particular path), and since it is a much different approach, which works for the specific

examples described here, it may at least be useful when these other methods fail. The Discrete Time Formulation

Let z€ ZC R” denote the exogenous state variables. If z, is Markovian, then its

evolution can be written as

7441 = Spey, 1),

where ¢€ t+1 € ~ consists of iid random variables drawn from a fixed cdf uw - [0,1], and

f:ZE 3 Z. Let z,€ XC R™ denote the endogenous state variables and ue UC R?

denote the "control" variables, which evolve according to

T= Wx ZU); g:XxZx U3 X,

“= h( 24,24), h:XxZ = U. The task at hand is to find an h€ H which satisfies a set of equilibrium conditions. For many problems these equilibrium conditions take one of two forms. Temporal

conditions can often be written as

B(2p2,u,) = 0, B:XxZxU 3 RP~4, whereas intertemporal conditions can be written as

MfpZyU,) = LMC y 4 2a yey 1) Mules, 1); N,M: Xx Zx U3 RY.

Using these equations, and the laws of motion defined above, define the function A(h),

A:H- A(#), such that (1) | Biz,z,(Ah)(z,z)] = 0, (2) Mz,z,(Ah)(2,2)] = [Mf g[2,2,(Ah)(2,2)], f(z,€),A[ 9 2,2,(Ah)(2,2)),f(2€)]}du(e).

The equilibrium control function is the fixed point Ah = A(h). Note that A(h) is an

argument of h, that A(h) is defined pointwise, and that A’s existence depends upon a

unique Brouwer fixed point. If A is continuous, isotone, strictly concave, and if it maps a nonempty, convex, and compact subset H into itself, then a unique fixed point h € H exists to which the sequence A”h) converges uniformly. The existence of a fixed point follows from Schauder’s fixed point theorem, and the uniqueness and convergence results can be found in Krasnoselskii and Zabreiko (1984).

There exists an intuitive reason for why this algorithm should work. Suppose that the time horizon for the control problem is 7, and that the equilibrium conditions above are required to hold for each period up to T—1. Suppose, also, that the control function is fixed at A at T(ie., up=A(t7z7)). The "optimal" control at T—1 is then A(A), and, given this control at T-—1, the optimal control at T—2 is A*(h), and so on. The further back you go, the less dependent should be the control on A, and thus A”(h) should converge to the optimal control.! This argument is not meant to substitute for a rigorous proof, especially since many algorithms which should work simply don’t. At this point, though, I think proofs are most usefully carried out for particular examples (e.g., see

Coleman (1987) for an application to the optimal growth and monetary growth models).

A Program for a Digital Computer

To implement this algorithm on a digital computer, clearly some modifications need to be made. First of all, H must be approximated by a finite dimensional set. This can be achieved by imposing a grid on (z,z) (heavily concentrated, say, around an estimate of the mean of the steady state distribution) and only requiring that equations (1)—(2) hold at these grid points. Since, though, h must be evaluated at arbitrary points of its domain (see eq. 2), an interpolation routine must be used. I have had success with multi—linear interpolation (in logs), but for small dimensional problems a more sophisticated

interpolater (e.g., splines) might be more efficient. With multilinear interpolation, H is

thus approximated by a finite dimensional set of piecewise linear functions, and if H is compact, this approximation is uniform.”

Another problem is the integration in equation (2). In general, an adaptive quadrature rule will compute the integral to within some prespecified tolerance, but this often proves to be too slow. A much faster route is to rely on a fixed—point quadrature rule based on the weight function dy. If dw takes on a standard form then existing routines may compute the points and weights (such as Hermite—Gauss quadrature for the Normal distribution, and clearly any distribution can be transformed into the Normal). In effect, this imposes a discrete probability model on e€ (but not on 2z), and thus the integration is approximated by a fixed summation.

Some extensions may be possible to speed up the convergence of A” The

generalization of Newton’s method as developed by Kantorovitch (1948),

hay y = hy [A’(h,) ACh);

which uses the Fréchet derivative A’, may speed things up, provided that this derivative (and its inverse) can be easily constructed. There may also exist a satisfactory extension of the multi—dimensional secant method (a discussion of the secant method can be found in

Ortega and Rheinboldt, 1970) which can be applied to this problem.

Ezamople 1 Consider the stochastic optimal growth model,® where a control function h is

sought which solves (some of the notation is switched)

i} t maxE% f#'u(c,), c,= A(z,z,), ner re A A subject to

741 = S(t p24) — ep where the exogenous state variable evolves according to 211 = Pt Epis €yy 1 is iid Normal(0,0°). The state vector is then s = (z,z), and eq. (2) becomes (eq. (1) does not apply)4 u’[(Ah)(z,2)] = Bfw’ {h[f(2,z) — (Ah)(2,2),p2 + ]}F’ [f(2,2) — (Ah)(2,2),p2 + €]du(e). The Continuous Time Formulation: The Deterministic Case

Suppose now that the control and state variables satisfy the ordinary differential

equation system a(t) = gx t),u(t)], g:XxU > R”™,

Ut) = flo(t),u(t)], FX R?,

Embedded in f are the equilibrium conditions which define the problem. Write the

equilibrium control as a function of the state variables,

u(t) = Ala(t)], A:X- U,

where

u(t) = hi [a(t)] a(t) = h’[2( t)]9{2(¢),A[(¢)]}.

Substitute in for u(t), and drop the time argument, to obtain

flz,h(2)] = he (2) g[2,h(2)],

which is a first order quasi—linear partial differential equation in A.

To obtain an approximate solution, note that

h’ (2) g[2,h(2)] » A{z + g[z,h(z)]d(z)} — A(z) 7 d(z)

]

where a(z) is some small increment step, say some fixed A. Define, then, the fixed point

equation as)?

fla(Ah)(a)] = 2 + gle, (4h) (z)]d(z)} — (4h) (2) nC)

Boundary conditions can easily be imposed by requiring that the initial h satisfy these conditions, and never updating A at the boundary. Unfortunately, one shortcoming of the above algorithm is that A(h) +h as d-0

whether or not A isa fixed point. To see this, write

flz,(Ah)(2)] = 2X2 + gles (4h) (2) Jd(z)} — A(z) 4 A(z) — (Ah) (2) ~ d(z) “day

Clearly, to prevent the last term from blowing up, (Ah)(z) must tend to h. (It is necessary that (Ah)(z) determines the derivative directly by entering the last term on the

right hand side, otherwise an incorrect guess of the derivative for h has too much

influence on the iterations.)! To construct a sequence which rapidly converges to the fixed

point, define a sequence of functions d+ 0 and corresponding fixed points A(h,) =h,,

where each set of iterations to find h, begins from the previous solution hoy: Example 2

Consider the deterministic continuous time optimal growth model,® where a control

c is sought which solves

ie] max fe eu c(t)} dt, c 0

subject to

The solution for the control can easily be determined to satisfy

et) = = Ba LF laa) -

Define the control function as c(t) = h[z(t)|, and derive the ordinary differential equation

The functional equation in Ah then becomes”

_w (CAA) (I pe(a) — gy = Me + LF(2) = (AA) (2)] d(2)} = (48) (2) TAY] eC

where the boundary condition is A(0) = 0.

The Continuous Time Formulation: The Stochastic Case

Suppose now that the state evolves according to the stochastic differential equation dx{t) = b[a{t),u(t)]dt + of2(t),u(t)]det), b:XxU4+R™, o:XxU 4 RY, where we wish to find an optimal control function u(t) = Ala(t)], A:X- U,

such that certain optimality conditions are obtained. For many problems, the partial differential equation of dynamic programming can, after some effort, be reduced to the

partial differential equation in the control

P{a,h(2),h’(2),A"(2), - + »,A)(a)} = 0, PXxte R™™M1)/2 RP a )(2), e.g, denotes all derivatives of order r) where boundary conditions are imposed on h. At this point, as expected, there is no clear cut way to proceed. If P corresponds to an analytic non—characteristic Cauchy initial value problem, then the nonlinear th order differential equation can, by differentiation, be reduced to a first order quasi—linear system which can be solved via the method developed above. Often (if not always), though, the problem at hand suggests a natural way to proceed in estimating the directional derivatives, and with these estimates the fixed point equation can be constructed analogously to the one constructed above. This is the route which I will follow in the

example below.

10

Ezample 3 Consider the continuous time stochastic optimal growth model, where a control

function A is sought which solves

max B Je Pale(ae, c(t) = hf2(z),0(2)],

subject to

dat) = {f[2(t),A(t)] — c(t) } dt, d&(t) = —AW(t)dt + odz(t). The differential equation then becomes (the derivation is in Appendix A)

nee eo — 6] = [f(2,0) — Ws, 6)]h, (2, 6)

— oh (2,0) + 220?" TA(2, 9) 11200 9) + 1 62h (2,0), 6 . 6 =o Moe

u"[h(z, 0) |

with the boundary conditions A(0,4) = h{0,6) = 0. This equation can be represented as the system (this does not correspond to the quasilinear system obtained by differentiating

P)

(2,0) 7 a 49),

eye 8) — p] + Oh g( 2,6) Ae Hs ) =

[F(a,6) — h( 2, 6)]h.(,0) + 794

11

The following approximations can then be used:

[F(2,) — h(2,0)]h,(2,0) = He + (F(z, 8) ~ A(2,4)]d(z), 0} — A(z, 0) d(z

h (2,8) — A{z,0 + d(6)} — A(z, 0)

— a9

and similarly for gp(z,0). Construction of the fixed point equation A(h,g) is then

straightforward. Note that, given the solution A above, Ito’s Lemma provides us with

de(t) = areca fat — plat — o acetate + halal t),6(¢)] ode t).

Some Demonstrations of the Algorithm’s Performance

Table 1 demonstrates the algorithm’s convergence for specific parameterizations of the above three examples. Appendix B contains the program to solve Example 1, the discrete time stochastic growth model. The entries in the table correspond to distances between succesive iterates of the control function. As can be seen, the succesive distances monotonically decrease, although the rate of decrease is much faster for the discrete time algorithm.

Concluding Remarks

I have had considerable success in applying the discrete time version of the algorithm described here to a wide diversity of problems. Given the prevalent use of discrete time models in economics and engineering, as well as in other disciplines, this algorithm could be useful in solving a wide variety of problems. The continuous time version, since it involves the additional approximation of a derivative, is probably less robust, but I have included a discussion of it here since I think it offers some new insight

into solving these types of systems.

iter. no. Ex. 12 1 0.7254 2 0.3227 3 0.1938 4 0.1205 5 0.0756 6 0.0471 7 0.0292

0.0178

g 0.0109 10 0.0066 ti 0.0039 12 0.0024 13 0.9014 14 0.0008 15 0.0005 16 0.0003 17 0.0062 18 0.0001 19 0.0000 20 0.0000

Table 1

sup |log ha) — log h,(s)| 8

dz=1.0

dr=0.05

Ex. ob 0.4899 0.1937 0.0964 0.0508 0.0303 0.0191 0.0118 0.0075 0.0048

0.0052 0.0043 0.0040 0.0038 0.0034 6.0030

0.0000 0.0000 0.0000

di= .20

dz=.10

dz=.05

h 0.6135 0.2373 0.1737 0.1303 0.1057 0.0890 0.0777 0.0685 0.0623

0.0013 0.0012 0.0011 0.0011 0.0011 0.0010

0.0007 0.0007 0.0006

Ex. 3° d0=.10

a0=.08

d0=.C5

a. = .95, uc) =—l/c, f(z,z) = e’Vz, p=.7, o =.1, initial A(z,z) = f(z,2).

lla

9 0.2106 0.0857 0.0331 0.0245 0.0210 0.0197 0.0166 0.0144 0.0121

0.0014 0.000-4 0.0003 0.0002 0.0002 0.0002

0.0002 0.0001 0.0001

b. p= .05, uc) = —1/c, f(z) = Jz + 2, initial h(z) = f(z) — 2. The -’s symbolize the

remaining iterations for the indicated values of dz.

c. p= 05, uc) =-I/c, f(a) =e e+ 2, X= .3, o=.1, initial h(2,0) = (0) —2,

(2,0) = 2h(2,0). The distances under column g correspond to levels, not logs.

V2

12

In terms of economic applications the solution methodology presented here may enhance cur ability to study sophisticated intertemporal general equilibrium models. Indeed, I think this technology can truly operationalize general equilibrium rational expectations models and put them on the same footing as the econometric models which

preceded them, and which are still in wide use today.

13

APPENDIX A

This appendix derives the solution to the continuous time stochastic growth model. Let v(t,z,6) denote the maximum value of the objective function starting at time ¢ with

at) and Ot). The solution solves the partial differential equation of dynamic

programming

—pt 1 2 (Al) “uy = max {ePlu(h) + u,(f—h) UO + 50" Ugg}. The first order condition for h is then

(A2) 0= e Phy: (h) — 0,

Differentiate (Al) with respect to z to obtain

—pt,, 1 2 (A3) —U,, =e Ply (A)A, + VA f— A) + 0(f.— h,) - U,9A9 + 50° 0, 9g

Using (A2), we can derive the following expressions for Yin Yap Yapg and Vee

—pt O=e? ul(h)h, — Von!

0= eP'ul(h)hy — Vp

0 = Pum(hyas + &Pu"(A) hag — Yapp

‘O= —pe Pty (h) — Uap

14

Substituting these derivatives of v into (A3) yields the desired differential equation in the

control h.

To obtain the equation for dc(t), use Ito’s Lemma to obtain _ 1 2 dc(t) = [A(f— h) — hp + 50 hoglat + hgodz{t). Use, then, the differential equation in the control to rewrite the term in brackets as

h,(f— h) — hy\8 + 50° hoy = - p) — yo ie

15 APPENDIX B

C

CRRA RIK K KKK EEK ERE KARE KIKI AK KEE RIERA IBREIIKIK III IIIS IIIS IS IIIS IIIS C PROGRAM TO SOLVE THE STOCHASTIC OPTIMAL GROWTH MODEL,

WRITTEN BY WILBUR JOHN COLEMAN II (APRIL, 1989).

* HRKK KKK K KEKE EER EERE RE RE KIERERE RRR EERE a casas accsssssscssesesese

*~

THE PROBLEM IS TO FIND THE FIXED POINT OF A, WHERE HP = A(H) IS DEFINED BY

UP[HP(S)] = BETA*E{UP[H(SP) ]*FP(SP) )

SP(1) = LOG[F(S) - HP(S)]

SP(2) = RHO*S(1) + QX, QX DISTRIBUTED NORMAL(0,SIMGA**2), LET NH = # OF FUNCTIONS THAT ARE BEING SOUGHT; HERE NH = 1. LET (NX,NZ) # OF (ENDOGENOUS , EXOGENOUS) STATE VARIABLES, AND NS = NX + NZ: HERE NX = = 1. LET {NSD(I),I=1,NS} = # OF GRID POINTS FOR CORRESPONDING STATE VARIABLES (S(I),I=1,NS}, SO THERE ARE NSP = NSD(1)*NSD(2)* * *NSD(NS) POSSIBLE SATES INDEXED BY {N=1,NSP}. LET {SMIN(I),SMAX(I)} CORRESPOND TO THE (MINIMUM,MAXIMUM) VALUES OF THE GRID POINTS FOR S(1). PLACE EQJIDISTANT GRID POINTS FOR S(1) INTO {SM(I),I=1,NSD(1)}, THOSE FOR $(2) INTO {SM(I) , I=NSD(1)+1,NSD(1)+NSD(2)}, AND SO ON. FOR EACH INDEX N THERE CORRESPONDS A POINTER K INTO SM SUCH THAT S(1) = SM(K(1)), S(2) = SM(NSD(1) + K(2)), AND SO ON, AND THERE CORRESPONDS A VALUE OF H AT EACH SUCH S WHICH IS STORED IN {HM(N,J),N=1,NSP} FOR EACH FUNCTION J=1,NH, AND SIMILARLY FOR HPM. DEFINE H AT AN ARBITRARY STATE BY MULITLINEAR INTERPOLATION BASED ON HM. THE ALGORITHM THUS COMPUTES HPM FOR A GIVEN H, COMPUTES THE DISTANCE BETWEEN HPM AND HM, PLACES HPM INTO HM, AND ITERATES UNTIL CONVERGENCE IS OBTAINED.

NZ

QAAaAaanannanananaanananaanaaaa

IMPLICIT REAL*8 (A-H,0-Z)

DIMENSION SM(1000) ,HM(10000,1) ,HPM(10000,1) , EPS(2) ,NSD(10) DIMENSION SUPH(10)

COMMON /COMN/SM,HM,HPM, EPS, BETA,NX,NZ,NS,NSD,NSP,NH,MAXIT WRITE(*,’(A35)') ' STOCHASTIC OPTIMAL GROWTH MODEL’

CALL INIT DO 300 NUMIT = 1, MAXIT CALL A DO 100 I = 1, NH 100 SUPH(I) = SUP(HM(1,1),HPM(1,I),NSP) WRITE(6,’(A10,15,A10,10F10.4)') * ‘A ITER. #',NUMIT, ‘NORM = ', (SUPH(1),I=1,NH) DO 200 I = 1, NH 200 IF (SUPH(I) .LE. EPS(1)) GOTO 400 300 CONTINUE

WRITE(*,'(A40)’) ' MAXIMUM NO. OF ITERATIONS EXCEEDED’ 400 CONTINUE

DO 500 I = 1, NH 500 WRITE(11) (HM(N,I),N=1,NSP)

STOP

END

C

CRRKK KKK K EKA RE EEE KEI KER EERE RE II OI II I I II I Esa SUBROUTINE A

C INPUT: FUNCTION H, ETC. OUTPUT: HPM

CHAKRA KHKK KKK KEELER ERK RE KKK EIEIO IO II I asa eesedese IMPLICIT REAL*8 (A-H,0-Z) DIMENSION SM(1000) ,HM(10000,1) ,HPM(10000,1) , EPS(2) ,NSD(10) DIMENSION K(10) ,S(10) ,HGUESS(10) ,HROOT(10) ,NSMPRD(10) , NSMSUM(10) COMMON /COMN/SM,HM,HPM, EPS, BETA,NX,NZ,NS,NSD,NSP,NH,MAXIT COMMON /CNSM/NSMPRD , NSMSUM

16 COMMON /CZ/S

EXTERNAL SC DO 100 I = 1, NS K(I) =1 100 S(1) = SM(1 + NSMSUM(1))

DO 200 N = 1, NSP HGUESS(1) = DLOG(HM(N,1)/(F(S) - HM(N,1))) CALL DNEQNF(SC,EPS(2) ,NH,MAXIT,HGUESS , HROOT, ZNORM)

C DNEQNF IS AN IMSL ROUTINE TO SOLVE NH NONLINEAR EQS. IN NH UNKNOWNS. C INPUT: SUBROUTINE SC(HGUESS,ZVEC,NH) WHICH EVALUATES ZVEC AT HGUESS, C CONV. CRITERION EPS(2), MAX. # OF ITER. MAXIT, GUESS OF ROOTS HCUESS. C OUTPUT: SOLUTION IN HROOT SUCH THAT ZVEC = 0, SUM OF SQUARES IN ZNORM. HPM(N,1) = F(S)/(1.0D0 + DEXP(-HROOT(1))) CALL UPDS(K,S) 200 CONTINUE RETURN END C

CR RRR RIK RIK IKK KEK EK EERE IR III III III III IO I I IIIT ae seak SUBROUTINE SC(HGUESS , ZVEC,NH)

C INPUT: GUESS OF TRANSFORMED HP(S) AS HGUESS, DIMENSION NH, STATE S. Cc QUPUT: ZVEC = RHS - LHS OF EACH FOC-EQUIL. COND. CR RRR RR EKK KK IK KKK LER ER ERIK IKE KI EKER I IK I I a a Ss csesek

IMPLICIT REAL*8 (A-H,O-Z) DIMENSION HGUESS(10) ,ZVEC(10) ,$(10) ,HP(10) COMMON /CZ/S 4P(1) = F(S)/(1.0D0 + DEXP(-HCUESS(1))) CALL Z(S,HP, ZVEC) 2ETURN END

c

CARH RR IK EK K RK KEK ERE RIK REIKI I ab ic iib babi ikakseseskskskt SUBROUTINE Z(S,HP,ZVEC)

C INPUT: STATE S, GUESS HP. OUTPUT: ZVEC.

CR RRER HK EIK KKK RK KEKE ERR KIRKE ROI kaise eK KE IMPLICIT REAL*8 (A-H,0-Z) JIMENSION SM(1000) ,HM(10000,1) ,HPM(10000,1) , EPS(2) ,NSD(10) DIMENSION $(10),SP(10),HP(10) ,ZVEC(10) ,QX(100) ,QWw(100) COMMON /COMN/SM,HM,HPM, EPS, BETA,NX,NZ,NS,NSD,NSP,NH,MAXIT COMMON /CPDF/SIGMA, RHO,QX,QW,NGAUSS

5SP(1) = DLOG(F(S) - HP(1)) ! SP IS NEXT PERIOD’S STATE = 0.0D0 DO 100 I = 1, NGAUSS ! INTEGRATE TO COMPUTE EXPECTED VALUE

SP(2) = RHO*S(2) + QX(1) W = W + UP(H(SP, 1))*FP(SP)*QW(I)

100 CONTINUE ZVEC(1) = BETA*W - UP(HP(1)) ! RHS - LHS RETURN END

C

CRAKE KKK KAR K KKK EK EKER KKK IKKE EERE EERE EER REED RE R RIM I A DOUBLE PRECISION FUNCTION UP(C)

C INPUT: C. OUTPUT: MARGINAL UTILITY.

CHK KARR AKER EKER KEKE RANKER KEKE REE EK ER II III IOI III IEE I EI E E ar ab ai sk sesk IMPLICIT REAL*8 (A-H,0-Z) COMMON /UTIL/TAU UP = 1.0D0/C**TAU RETURN END

C CHER RKRAERE KKK KK KKK EERE KERR EERE RE RRER ERR EREREREREEEREREREREREERERREERRER DOUBLE PRECISION FUNCTION F(S) C INPUT STATE S. OUTPUT: PRODUCTION. CER KRRKEKE KEKE ERR EERE ER REE ERE RERERERR ERR ER RR ERRERREREERREREREERERERERERERE IMPLICIT REAL*8 (A-H,0-Z) DIMENSION $(10) COMMON /PROD/ALPHA, DELTA F = DEXP(S(2))*DEXP(S(1))**ALPHA + DELTA*DEXP(S(1)) RETURN END C CRERKKEK EKER KKK KEE KEKE EERE BERK REE RE EKER RRR IKE RE REE ERR RRERERRRREREEERE DOUBLE PRECISION FUNCTION FP(S) INPUT: STATE S. OUTPUT: MARGINAL PRODUCTIVITY. EFDA EEE IER IKK KEKE KEKE KER ERK ERE KEKE ERE ERERERERERERERRRBRERR ER IMPLICIT REAL*8 (A-H,0-Z) DIMENSION S(10) COMMON /PROD/ALPHA , DELTA FP = ALPHA*DEXP(S(2))*DEXP(S(1))**(ALPHA - 1.0D0) + DELTA RETURN END

sles t. RWW

C Cx

EAE HE IEA HK HIE IE IE IEF IE EIEN TERK IEE IED IE ERIE TE EE IEE IE ETE EER HIRI DOUBLE PRECISION FUNCTION H(S,1)

C INPUT: STATE S, FUNCTION # I. OUTPUT: INTERPOLATED FUNCTION VALUE.

CHK AHHH AANA AK KAKA EKA K ANKE EEE K KARA EERIE EINER RIKI EE ETEK ET TE

IMPLICIT REAL*8 (A-H,0O-Z)

DIMENSION SM(1000) ,HM(10000,1) ,HPM(10000,1) , EPS(2) ,NSD(10)

DIMENSION S(1)

COMMON /COMN/SM,HM,HPM,EPS , BETA,NX,NZ,NS,NSD,NSP,NH,MAXIT

H = GRID(S,SM,HM(1,1),NS,NSD)

RETURN

END

vt

CREEK KKK EKER KKK KRE KR REEREKEERERRERERERERRERR ERR ERE ER RERRERERERRERREREERE DOUBLE PRECISION FUNCTION SUP(FO,F1,NSP)

C INPUT: VECTORS FO, Fl. OUTPUT: SUP DISTANCE, AND Fl => FO.

RR KK REE KER ERE KER ERE ERE RRR ER RE RR ER ER ER ERR ER ERR ERE ERE EEREREE IMPLICIT REAL*8 (A-H,0-Z) DIMENSION FO(1) ,F1(1) SUP = 0.0D0 DO 100 N = 1, NSP

SUP = DMAX1(SUP, DABS(DLOG(F1(N)) - DLOG(FO(N))))

100 FO(N) = F1(N) RETURN END

Cc

CREKEKEKKE KR EKRRERERER EERE SUBROUTINE UPDS(K,S) C INPUT: POINTER K. OUTPUT: NEXT K, CORRESPONDING STATE S. CRRRRAKEK ERK REI ERE EERE EERE EE EE ERE EERE REE ERE EERE RE RERE BER IMPLICIT REAL*8 (A-H,0-Z) DIMENSION SM(1000) ,HM(10000,1) ,HPM(10000,1) ,EPS(2) ,NSD(10) DIMENSION K(10) ,$(10) ,NSMPRD(10) ,NSMSUM(10) COMMON /COMN/SM,HM,HPM, EPS, BETA,NX,NZ,NS,NSD,NSP,NH,MAXIT COMMON /CNSM/NSMPRD ,NSMSUM K(1) = K(1) +1 po 100 I=1, NS - 1

KAKA HAHAHA AKAIKE IKI III III IIIA IIE

18 IF (K(I) .LE. NSD(I)) GOTO 200

K(I) = 1 K(I+1) = K(I+1) +1 100 CONTINUE 200 CONTINUE DO 300 I = 1, NS

300 S(I) = SM(K(I) + NSMSUM(I)) RETURN END

Cc

CORBI EIR EIR EIS IDO OSES I bbbbbbiabsbisisksksbsese SUBROUTINE INIT

C INITIALIZE PARAMETERS, ETC.

CARRIE IRI IRSA RAI I II I bbbbbbbskbbssesessk IMPLICIT REAL*8 (A-H,0-Z) DIMENSION SM(1000) ,HM(10000,1) ,HPM(10000,1), EPS(2) ,NSD(10) DIMENSION QX(100) ,QW(100) ,NSMPRD(10) , NSMSUM(10) DIMENSION K(10),S(10) , SMIN(10) , SMAX(10) COMMON /COMN/SM,HM,HPM, EPS, BETA,NX,NZ,NS,NSD,NSP,NH,MAXIT COMMON /UTIL/TAU COMMON /PROD/ALPHA, DELTA COMMON /CPDF/SIGMA ,RHO,QX,QW,NGAUSS COMMON /CNSM/NSMPRD , NSMSUM

C READ(10,'(10X,110)') NH,NX,NZ NS = NX + NZ READ(10,'(10X,F10.4)') BETA, TAU,RHO,SIGMA,ALPHA, DELTA, EPS, * (SMIN(I) , SMAX(I) ,I=1,NS) READ(10,'(10X,110)') (NSD(I),I=1,NS) ,NGAUSS ,MAXIT C NSMPRD(1) = 1 ! PARTIAL PRODUCTS OF NSD NSMSUM(1) = 0 ! PARTIAL SUMS OF NSD

DO 100 I = 2, NS NSMPRD(I) = NSMPRD(I - 1)*NSD(I - 1) NSMSUM(I) = NSMSUM(I - 1) + NSD(I - 1) 100 CONTINUE NSP = NSMPRD(NS)*NSD(NS)

C SCALE = 1.0D0/DSQRT(DCONST('PI’)) CALL DGQRUL(NGAUSS ,4,A1,B1,0,QXFIX,QX, QW) C DGQRUL IS AN IMSL ROUTINE TO COMPUTE HERMITE-GAUSS QUADRATURE POINTS C AND WEIGHTS FOR A NORMAL(0,1) DISTRIBUTION. C INPUT: # OF QUAD. POINTS NGAUSS. OUTPUT: POINTS QX, UNSCALED WEIGHTS QW. DO 200 I = 1, NGAUSS QX(1I) = SIGMA*QX(1) 200 QW(I) = SCALE*QW(1) Cc N= 0 ! INITIALIZE SM DO 400 J = 1, NS N=N+1 SM(N) = SMIN(J) DO 300 I = 2, NSD(J) N=N+41 SM(N) = SMIN(J) * + (DFLOAT(I - 1)/DFLOAT(NSD(J) - 1))*(SMAX(J) - SMIN(J)) 300 CONTINUE 400 CONTINUE C ! INITIALIZE HM

DAPC = 1.0D0 - ALPHA*BETA/(1.0D0 - (1.0D0 - ALPHA)*DELTA*BETA)

19

DO 500 I = 1, NS K(I) = 1 500 S(I) SM(1 + NSMSUM(I))

DO 600 N = 1, NSP HM(N,1) = DAPC*F(S) CALL UPDS(K,S)

600 CONTINUE RETURN END

C

CEKKKEKEKAK KEKE KEKE KKK KEKE KK KEE KKK HEE KKK KBE EKER HH IEEE EE EEE EEE EEE DOUBLE PRECISION FUNCTION GRID(S,SM,GM,NS,NSD) C INPUT: STATE S$, STATE GRID SM, FUNCTION VALUES GM, DIMENSIONS NS, NSD. C OUTPUT: MULTILINEAR INTERPOLATION OF FUNCTION AT S BASED ON GM AT SM. CERKKEKEKAKEK KKK EERE AKER ER EKER ERE EERE EEE ERE RE EEE KBE EE IE IMPLICIT REAL*8 (A-H,0-Z) DIMENSION S(1),SM(1),GM(1) ,NSD(1) ,KS(10) ,KS2(10) ,T(10) , IP(10) DIMENSION NSMPRD(10) ,NSMSUM(10) COMMON /CNSM/NSMPRD , NSMSUM

Cc DO 100 I = 1, NS ! FIND SURROUNDING GRID POINTS J = 1 + NSMSUM(I) CALL LOCATE(SM(J) ,NSD(I) ,S(I) ,KS(I)) IF (NSD(I) .EQ. 1) THEN T(I) = 1.0D0 ELSE ! COMPUTE POSITION OF STATE RELATIVE TO GRID POINTS T(I) = (SM(J+KS(I)) - S(1))/(SM(J4+KS(I)) - SM(J+KS(I)-1)) ENDIF IP(I) = 0 100 CONTINUE Cc GRID = 0.0D0 DO 400 N = 1, 2**NS ! LOOP OVER ALL CORNERS OF NS DIMENSIONAL CUBE TP = 1.0D0 DO 200 J = 1, NS ! FIND WEIGHT FOR A CORNER IF (IP(J) .EQ. 0) THEN TP = TP*T(J) ELSE TP = TP*(1.0D0 - T(J)) ENDIF IF (NSD(J) .EQ. 1) THEN ! PLACE CORNER INTO KS2 KS2(J) = KS(J) ELSE KS2(J) = KS(J) + IP(J) ENDIF 200 CONTINUE ! ADD WEIGHT TIMES VALUE OF CORNER GRID = GRID + TP*DLOG(GM(KMAP(KS2 ,NS) )) C IP(1) = IP(1) +1 ! UPDATE CORNER DO 300 J = 1, NS - 1 IF (IP(J) .LE. 1) GOTO 400 IP(J) = 0 300 IP(J +1) = IP(J +1) +1 400 CONTINUE GRID = DEXP(GRID) RETURN END C

CHAKA KKAKAA KEE KEKE ERIKA KEK ERR ERE ERIE ER ER ER ERR RRR ER EERIE ERE RR ERE IEE

C

20 FUNCTION KMAP(K,NS)

INPUT: POINTER K, DIMENSION NS. OUTPUT: INDEX N AS KMAP.

CR REI IIE I I I I ISS SS sbi sbksbieaboeaksiesk

100

Cc

IMPLICIT REAL*8 (A-H,0-Z) DIMENSION K(1),NSMPRD(10) ,NSMSUM(10) COMMON /CNSM/NSMPRD , NSMSUM KMAP = K(1) DO 100 I = 2, NS KMAP = KMAP + NSMPRD(I)*(K(I) - 1) RETURN END

CARRE RI RRR ERI RR ROE I I OS I SS abbibbibsssok esses

C C C

100

HHKEK

SUBROUTINE LOCATE(XV,N,X,J)

INPUT: INCREASING VECTOR XV, DIMENSION N, VALUE X.

OUTPUT: LARGEST J SUCH THAT XV(J) < X, WHERE 0 < J <N.

REAR RE KKK ERK I REE RRR IO CIO IIS bbb sbeskakssleak IMPLICIT REAL*8 (A-H,0-Z)

DIMENSION XV(1)

JL =1 ! SET BOUNDS JU=N CONTINUE IF (JU - JL .GT. 1) THEN JM = (JU + JL)/2 ! COMPUTE MIDPOINT IF (X .GT. XV(JM)) THEN ! LOCATE X RELATIVE TO MIDPOINT JL = JM ! MIDPOINT BECOMES LOWER BOUND ELSE JU = JM ! MIDPOINT BECOMES UPPER BOUND ENDIF GOTO 100 ENDIF J = JL ! RETURN LOWER INDEX RETURN END

SAMPLE DATA FILE

NH

NX

NZ BETA TAU RHO SIGMA ALPHA DELTA EPS1 EPS2 XMIN XMAX ZMIN ZMAX NSD1 NSD2 NGAUSS MAXIT

SF OOMNrFRFRFOODOOONSO co} (o>)

ENDNOTES

1. Robert Lucas suggested something like this to me a few years ago, and subsequent to this discussion I realized that a version of this story works for the finite horizon optimal growth setup. Marianne Baxter’s paper then clarified the discussion for me (also in the context of the optimal growth setup), which led to the generalization of the story presented here.

2. That is, within the set of piecewise linear functions a sequence can be constructed which converges uniformly to any given he H.

3. This problem (with iid shocks) is described in Brock and Mirman (1972).

4. In Coleman (1987) I construct a set H which is convex and compact, and a function A:H+H which is continuous, thus guaranteeing the existence of a fixed point. In the deterministic setting, and with a specific assumption on the utility function, I essentially prove A’s monotonicity and strict concavity, thus proving the uniform convergence of A"(h) to the unique fixed point.

5. The intuitive argument for convergence given above also has a loose analogue here. In this case, h fixes the control for tomorrow, and A(h) is the optimal control which takes into account both the effect on the state today and the direction of the state into tomorrow.

6. Approximating A’(z)g[z,A(z)] via partial derivatives results in an unnecessary interaction between A(h) in the derivative and A(h) in g (eg., if g is linear then a quadratic in A(h) appears).

7. In the following example, if fl) = 0 and if h’(1) does not equal the true value

exactly, then (Ah)(1) =0 if A(h) does not enter the derivative as above. The correct

21

22

solution is A(1) > 0, and the above algorithm results in (Ah)(1) > 0. 8. This problem was solved by Cass (1965) and Koopmans (1965). 9. Suppose —u’(c)/u"(c) is increasing in c, then, if A is increasing, and if A is small,

A(h) is uniquely defined since

A y,t) = —" (y) Lf’ (2) — p] —- h{x + (f(z) - yJA} -y t y rs rs

u(y) is increasing in y (these statements need only be true for some large upper bound on 4g, say Z such that f(z) = 0). Provided, again, that A is small, then A(h) will be

increasing since 2(y,z) is decreasing in z. These results, though, do not guarantee

convergence or existence, but only that A(h) is well defined at each step of the iteration.

REFERENCES

Baxter, Marianne, 1988, Approximating Suboptimal Dynamic Equilibria: An Euler Equation Approach, Rochester Center for Economic Research Working Paper 139.

Brock, William A. and Leonard J. Mirman, 1972, Optimal Economic Growth and Uncertainty: The Discounted Case, Journal of Economic Theory 4, 497-513.

Cass, D., 1965, Optimum Growth in an Aggregative Model of Capital Accumulation, Reveiw of Economic Studies 32, .

Coleman, Wilbur John II, 1987, Money, Interest, and Capital, unpublished Ph.D. Dissertation, The University of Chicago.

Courant, R. and D. Hilbert, 1962, Methods of Mathematical Physics, Vol. II (John Wiley and Sons, New York).

Fleming, Wendell H. and Raymond W. Rishel, 1975, Deterministic and Stochastic Optimal Control (Springer—Verlag, Berlin).

John, Fritz, 1971, Partial Differential Equations (Springer-Verlag, New York).

Kantorovitch, L. V., 1948, Functional Analysis and Applied Mathematics, Uspekhi Matematicheskekh Nauk 3, no. 6, 89-185, translated from Russian by: Curtis D. Benster, 1952, National Bureau of Standards Report 1509.

Koopmans, T. C., 1965, On the Concept of Optimal Economic Growth, in: The Econometric Approach to Development Planning, Pontificiae Academiae Scientiarvm Scriptvm Varia (North—Holland, Amsterdam).

Krasnosel’skii, M. A. and P. P. Zabreiko, 1984, Geometrical Methods of Nonlinear Analysis (Springer—Verlag, Berlin).

Ortega, J. M. and W. C. Rheinboldt, 1970, Iterative Solution of Nonlinear Equations in Several Variables (Academic Press, New York).

Taylor, John B. and Harald Uhlig, 1989, Solving Nonlinear Stochastic Growth Models: A

Comparison of Alternative Solution Methods, unpublished working paper, The University of Minnesota and The Institute for Empirical Macroeconomics.

23

IFDP NUMBER

351

350

349

348

347

346

345

344

343

342

341

340

339

338

337

Please: address requests for co Papers, Division of International F Federal Reserve System, Washington, D.C.

International Finance Discussion Papers

TITLES 1989

An Algorithm to Solve Dynamic Models

Implications of the U.S. Current Account Deficit

Financial Integration in the European Community

Exact and Approximate Multi-Period Mean-Square Forecast Errors for Dynamic Econometric Models

Macroeconomic Policies, Competitiveness, and U.S. External Adjustment

Exchange Rates and U.S. External Adjustment in the Short Run and the Long Run

U.S. External Adjustment:

Progress and Prospects

Domestic and Cross-Border Consequences of U.S. Macroeconomic Policies

The Profitability of U.S. Intervention Approaches to Managing External Squilibria: Where We Are, Where We Might Be Headed, and How We Might

Get There

A Note on "Transfers" A New Interpretation of the Coordination Problem and its Empirical Significance

A Long-Run View of the European Monetary System

1988

The Forward Exchange Rate Bias: Explanation

A New

Adequacy of International Transactions and Position Data for Policy Coordination

20551.

24

AUTHOR(s

Wilbur John Coleman II

David H. Howard

Sydney J. Key

Neil R. Ericsson Jaime R. Marquez

Peter Hooper

Peter Hooper

William L. Helkie Peter Hooper

Ralph C. Bryant John Helliwell Peter Hooper

Michael P. Leahy

Edwin M. Truman

David B. Gordon Ross Levine

Matthew B. Canzoneri Hali J. Edison

Hali J. Edison Eric Fisher

Ross Levine

Lois Stekler

pies to International Finance Discussion inance, Stop 24, Board of Governors of the

IFDP NUMBER 336

335

334

333 332

331

330 329 328

327.

International Finance Discussion Papers TITLES

Nominal Interest Rate Pegging Under Alternative Expectations Hypotheses

The Dynamics of Uncertainty or The Uncertainty of Dynamics: Stochastic J-Curves

Devaluation, Exchange Controls, and Black Markets for Foreign Exchange in Developing Countries

International Banking Facilities Panic, Liquidity and the Lender of Last Resort: A Strategic Analysis Real Interest Rates During the

Disinflation Process in Developing Countries

‘International Comparisons of Labor

Costs in Manufacturing

Interactions Between Domestic and Foreign Investment

The Timing of Consumer Arrivals in Edgeworth’s Duopoly Model

Competition by Choice

25

AUTHOR(s)

Joseph E. Gagnon Dale W. Henderson

Jaime Marquez

Steven B. Kamin

Sydney J. Key Henry S. Terrell

R. Glen Donaldson

Steven B. Kamin David F. Spigelman

Peter Hooper Kathryn A. Larin

Guy V.G. Stevens Robert E. Lipsey

Marc Dudey

Marc Dudey