c PROGRAM: FD2.F c c DESCRIPTION: implementation of a 2-D finite-difference method for c solving the dynamic program associated with the c continuous-time stochastic growth model c c ECONOMIC the continuous-time stochastic growth model with c MODEL: agents choosing c(t) to maximize c c E[ int exp^(-rho*t) u(c(t)) dt | k(0),K(0)] c c s.t. dk = ((1-tau)*r(K)+w(K)-c)*dt + sig(k,K)*dz c dK = (F(K)-C(K)-G(K))*dt + sig(K,K)*dz c c FUNCTIONAL u(c) = c^omega/omega c FORMS: F(K) = lambda*K^theta-delta*K c G(K) = tau*F'(K)*K c c PROBLEM: find v(t,k,K) for t=infinity that solves the functional c equation: c c -dv/dt + rho*v = (1-omega)/omega c * (dv/dk)^(omega/(omega-1)) c + dv/dk*[(1-tau)*(theta*lambda*K^(theta-1) c -delta)*k)+lambda*K^theta*(1-theta)] c + dv/dK*[lambda*K^theta-delta*K c -(dv(K,K)/dk)^(1/(omega-1)) c -tau*(theta*lambda*K^theta-delta*K)] c + 1/2*[d^2v/dk^2 * sig(k,K)^2 + 2*d^2v/dkdK * c *sig(k,K)*sig(K,K) + d^2v/dK^2 *sig(K,K)^2] c c ALGORITHM: finite difference -- point implicit method c c TEST CASE: theta=omega=1/2,tau=0 c ==> back to case in fd.f c c INPUT FILE: (a) fd2.inp: the following need to be specified c c Input Description Format c ------------------------------------------------------ c cfl CFL number real > 0 c delta depreciation rate real in [0,1] c lambda production scale real > 0 c omega utility parameter real <= 1 c rho discount rate real in (0,1) c sigma technology parameter real >= 0 c tau tax rate real >= 0 c theta capital share real >= 0 c initialize option to provide 1 (=yes) or 0 c initial guess for c c v initial guess for il x 1 real c value function vector c (if initialize=1) c c OUTPUT FILES: (b) fd2.nxt: new input file with solution for v(k,K) c (c) fd2.dat: output file -- read into Matlab file fd2.m c c REFERENCES: Candler (notes) c c Graham Candler, 8-30-96 c Revised, 2-23-97, ERM c program fd2 implicit real(a-h,o-z) parameter(nx=81,ny=81,itmax=10000) real x(nx),y(ny),v(nx,ny),a1(nx,ny),a2(nx,ny),source(nx,ny), & dv(nx,ny) c c Open input and output files. c open(unit=5, file='fd2.inp') open(unit=7, file='fd2.nxt') open(unit=8, file='fd.dat') c c Read in parameters and, if specified, the initial guess for the c value function. c read(5,*) cfl read(5,*) delta read(5,*) lambda read(5,*) omega read(5,*) rho read(5,*) sigma read(5,*) tau read(5,*) theta read(5,*) initialize crit = 0.00001 il = nx jl = ny ilm1 = il-1 jlm1 = jl-1 c1 = (1.0-omega)/omega c2 = 1.0/(omega-1.0) c3 = omega*c2 tau1 = 1.0-tau rks = ((rho+delta*tau1)/(theta*lambda*tau1))**(1.0/(theta-1.0)) vks = (lambda*rks**theta-delta*rks)**(omega-1.0) c c Initialize v(.). c if (initialize.eq.0) then vv = (c1*vks**c3+vks*lambda*rks**theta-vks*delta*rks)/rho vv = 12.0 do 10 j=1,jl y(j) = 2.0*rks*float(j-1)/float(jl-1) do 10 i=1,il x(i) = 2.0*rks*float(i-1)/float(il-1) dv(i,j) = 0.0 v(i,j) = vv*sqrt(x(i)) 10 continue else do 20 j=1,jl y(j) = 2.0*rks*float(j-1)/float(jl-1) do 20 i=1,il x(i) = 2.0*rks*float(i-1)/float(il-1) dv(i,j) = 0.0 read(5,*) v(i,j) 20 continue endif c c Iterate in time until v(.) converged. c dt = cfl do 30 it=2,itmax resid = 0.0 do 40 j=2,jlm1 do 40 i=2,ilm1 vx = (v(i+1,j)-v(i-1,j))/(x(i+1)-x(i-1)) vxk = (v(j+1,j)-v(j-1,j))/(x(j+1)-x(j-1)) a1(i,j) = -(c1*vx**c2+tau1*(lambda*theta*y(j)**(theta- & 1.0)-delta)*x(i)+lambda*y(j)**theta*(1.0- & theta)+tau*(theta*lambda*y(j)**theta- & delta*y(j))) a2(i,j) = -(lambda*y(j)**theta-delta*y(j)-vxk**c2) source(i,j) = -rho*v(i,j) 40 continue do 50 j=2,jlm1 do 50 i=2,ilm1 vxp = (v(i+1,j)-v(i,j))/(x(i+1)-x(i)) vxm = (v(i,j)-v(i-1,j))/(x(i)-x(i-1)) vyp = (v(i,j+1)-v(i,j))/(y(j+1)-y(j)) vym = (v(i,j)-v(i,j-1))/(y(j)-y(j-1)) dx = 0.5*(x(i+1)-x(i-1)) dy = 0.5*(y(j+1)-y(j-1)) alpx = 0.5*(a1(i,j)-abs(a1(i,j))) almx = 0.5*(a1(i,j)+abs(a1(i,j))) alpy = 0.5*(a2(i,j)-abs(a2(i,j))) almy = 0.5*(a2(i,j)+abs(a2(i,j))) rLHS = 1.0+dt/dx*almx-dt/dx*alpx+dt/dy*almy- & dt/dy*alpy+dt*rho RHS = dt*(source(i,j)-almx*vxm-alpx*vxp- & almy*vym-alpy*vyp) dv(i,j) = RHS/rLHS 50 continue do 60 j=2,jlm1 do 60 i=2,ilm1 v(i,j) = v(i,j)+dv(i,j) resid = resid+dv(i,j)*dv(i,j) 60 continue do 70 j=1,jl v(il,j) = v(ilm1,j)+vks*(x(il)-x(ilm1)) 70 continue do 80 i=1,il v(i,1) = 0.0 v(i,jl) = v(i,jlm1) 80 continue resid = sqrt(resid)/float((il-2)*(jl-2)) if (mod(it,100).eq.0) write(6,*) resid if (resid.lt.crit) go to 999 30 continue 999 if (it.ge.itmax) then write(*,*) 'Maximum number of iterations reached' write(*,*) 'Current residual = ',resid endif c c Write out results. c write(7,'(2X,F12.6,10X,''cfl >= 0 '')') cfl write(7,'(2X,F12.6,10X,''delta in [0,1] '')') delta write(7,'(2X,F12.6,10X,''lambda >= 0 '')') lambda write(7,'(2X,F12.6,10X,''omega <= 1 '')') omega write(7,'(2X,F12.6,10X,''rho >= 0 '')') rho write(7,'(2X,F12.6,10X,''sigma >= 0 '')') sigma write(7,'(2X,F12.6,10X,''tau in [0,1] '')') tau write(7,'(2X,F12.6,10X,''theta in (0,1) '')') theta write(7,'(13X,''1'',10X,''read initial guess for v (1=yes)'')') write(7,'(2X,F12.6,10X,''value function by column'')') v(1,1) do 120 i=2,il write(7,1) v(i,1) 120 continue do 130 j=2,il do 130 i=1,il write(7,1) v(i,j) 130 continue phi = -rho*c2-delta*c3 do 140 i=2,ilm1 vx = (v(i+1,i)-v(i-1,i))/(x(i+1)-x(i-1)) write(8,*) x(i)/rks,vx**c2,phi*x(i) 140 continue do 150 i=2,ilm1 write(8,*) x(i)/rks,v(i,i),phi**(omega-1.0)*x(i)**omega/omega+ & lambda*phi**(omega-1.0)/rho 150 continue 1 format(2X,F12.6) stop end