c PROGRAM: FD.F c c DESCRIPTION: implementation of a 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) c(t)^omega/omega dt | k(0)] c c s.t. dk = (lambda*k^theta-delta*k-c)*dt + sigma*k*dz c c PROBLEM: find v(t,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*(lambda*k^theta-delta*k) c + d^2v/dk^2* sigma^2 * k^2 /2 c c ALGORITHM: finite difference -- explicit, implicit, and c diagonalized implicit method c c TEST CASE: theta=omega=1/2 c ==> v=phi^(omega-1)*k^omega/omega+lambda*phi^(omega-1)/rho c c=phi*k, phi=(rho+delta*omega)/(1-omega) c c INPUT FILE: (a) fd.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 theta capital share real >= 0 c initialize option to provide 1 (=yes) or 0 c initial guess for c c x capital grid il x 1 real vector c v initial guess for il x 1 real vector c value function c (if initialize=1) c c OUTPUT FILES: (b) fd.nxt: new input file with solution for v(k) c (c) fd.dat: output file -- read into Matlab file fd.m c c REFERENCES: Candler (notes) c c Graham Candler, 8-30-96 c Revised, 2-22-97, ERM c program fd implicit real(a-h,o-z) parameter(nx=81,itmax=10000) real x(nx),v(nx),a1(nx),b1(nx),source(nx),dv(nx),lambda common/a/ a(nx),b(nx),c(nx),d(nx) c c Open input and output files. c open(unit=5, file='fd.inp') open(unit=7, file='fd.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,*) theta read(5,*) implicit read(5,*) initialize crit = 0.000001 il = nx ilm1 = il-1 c1 = (1.0-omega)/omega c2 = 1.0/(omega-1.0) c3 = omega*c2 rks = ((rho+delta)/lambda/theta)**(1.0/(theta-1.0)) vks = (lambda*rks**theta-delta*rks)**(omega-1.0) do 5 i=1,il read(5,*) x(i) x(i) = rks*x(i) dv(i) = 0.0 5 continue c c Initialize v(.). c if (initialize.eq.0) then vv = (c1*vks**c3+vks*lambda*rks**theta-vks*delta*rks)/rho do 10 i=1,il v(i) = vv*sqrt(x(i)) 10 continue else do 20 i=1,il read(5,*) v(i) 20 continue endif c c Iterate in time until v(.) converged. c do 30 it=2,itmax resid = 0.0 do 40 i=2,ilm1 vx = (v(i+1)-v(i-1))/(x(i+1)-x(i-1)) a1(i) = -(c1*vx**c2+lambda*x(i)**theta-delta*x(i)) b1(i) = -(0.5*sigma*sigma*x(i)*x(i)) source(i) = -rho*v(i) 40 continue dt = 1.0 do 50 i=2,ilm1 dx = 0.5*(x(i+1)-x(i-1)) dt = amin1(dx/abs(a1(i)),dt) 50 continue dt = cfl*dt c c using the EXPLICIT METHOD ... c if (implicit.eq.0) then do 60 i=2,ilm1 dx = 0.5*(x(i+1)-x(i-1)) vx = (-v(i)+v(i+1))/dx vxx = (v(i+1)-2.0*v(i)+v(i-1))/dx/dx RHS = dt*(source(i)-a1(i)*vx-b1(i)*vxx) dv(i) = RHS 60 continue do 70 i=2,ilm1 v(i) = v(i) + dv(i) resid = resid + dv(i)*dv(i) 70 continue endif c c using the IMPLICIT METHOD ... c if (implicit.eq.1) then do 80 i=2,il-1 dx = 0.5*(x(i+1)-x(i-1)) vxp = (v(i+1)-v(i))/(x(i+1)-x(i)) vxm = (v(i)-v(i-1))/(x(i)-x(i-1)) vxx = (v(i+1)-2.0*v(i)+v(i-1))/dx/dx ap = 0.5*(a1(i)+abs(a1(i))) am = 0.5*(a1(i)-abs(a1(i))) RHS = dt*(source(i)-ap*vxm-am*vxp-b1(i)*vxx) a(i) = dt/dx*am+dt/dx/dx*b1(i) b(i) = -dt/dx*ap+dt/dx/dx*b1(i) d(i) = 1.0 + dt/dx*(ap-am)-2.0*dt/dx/dx*b1(i)+dt*rho c(i) = RHS 80 continue call sy(2,ilm1) c call tridiag(2,ilm1) do 90 i=2,ilm1 v(i) = v(i) + c(i) resid = resid + c(i)*c(i) 90 continue endif c c using the DIAGONALIZED IMPLICIT METHOD ... c if (implicit.eq.2) then do 100 i=2,ilm1 dx = 0.5*(x(i+1)-x(i-1)) vx = (-v(i)+v(i+1))/dx vxx = (v(i+1)-2.0*v(i)+v(i-1))/dx/dx rLHS =-a1(i)*dt/dx+dt*rho-2.0*dt*b1(i)/dx/dx RHS = dt*(source(i)-a1(i)*vx-b1(i)*vxx) dv(i) = RHS/rLHS 100 continue do 110 i=2,ilm1 v(i) = v(i) + dv(i) resid = resid + dv(i)*dv(i) 110 continue endif v(il) = v(ilm1)+vks*(x(il)-x(ilm1))*0.75 resid = sqrt(resid)/float(il-2)/dt if (mod(it,100).eq.0) write(*,*) it,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,''theta in (0,1) '')') theta write(7,'(2X,I12,10X, & ''method? (0=explicit,1=implicit,2=diagonlized implicit)'')') & implicit write(7,'(13X,''1'',10X,''read initial guess for v (1=yes)'')') write(7,*) write(7,'(2X,F12.6,10X,''capital stock grid'')') x(1)/rks do 120 i=2,il write(7,1) x(i)/rks 120 continue write(7,*) write(7,'(2X,F12.6,10X,''value function '')') v(1) do 130 i=2,il write(7,1) v(i) 130 continue phi = -rho*c2-delta*c3 do 140 i=2,ilm1 vx = (v(i+1)-v(i-1))/(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),phi**(omega-1.0)*x(i)** & omega/omega+lambda*phi**(omega-1.0)/rho 150 continue 1 format(2X,F12.6) stop end subroutine sy(il,iu) parameter(nx=81) common/a/ a(nx),b(nx),c(nx),d(nx) lp = il+1 do 10 i=lp,iu r = b(i)/d(i-1) d(i) = d(i)-r*a(i-1) c(i) = c(i)-r*c(i-1) 10 continue c(iu)=c(iu)/d(iu) do 20 i=lp,iu j = iu-i+il c(j) = (c(j)-a(j)*c(j+1))/d(j) 20 continue return end subroutine tridiag(n1,n) parameter(nx=81) common/a/ a(nx),b(nx),c(nx),d(nx) d(n1) = 1.0/d(n1) a(n1) = c(n1)*d(n1) n2 = n1+1 n1n = n1+n do 10 k=n2,n k1 = k-1 b(k1) = b(k1)*d(k1) d(k) = d(k)-a(k)*b(k1) d(k) = 1./d(k) a(k) = (c(k)-a(k)*a(k1))*d(k) 10 continue c(n) = a(n) do 20 k1=n2,n k = n1n-k1 c(k) = a(k)-b(k)*c(k+1) 20 continue return end