c PROGRAM: SGM1DL.F c c DESCRIPTION: implementation of finite-element method for stochastic c growth model, 1-D, linear interpolants c c ECONOMIC the discrete-time stochastic growth model with c MODEL: agents choosing c(t) to maximize c c E[ sum beta^t c[t]^(1-tau)/(1-tau) | k_0] c t c subj. to c k[t+1] = theta[t]*k[t]^alpha+(1-delta)*k[t]-c[t] c ln(theta) drawn from a Markov chain c c[t],k[t] >= 0 c i[t] >= 0 for sufficiently large value of gamma c c PROBLEM: find c=c(x,i) that solves the functional equation: c c c^{-tau} = beta pi(i,j) c~^{-tau}(alpha* c *exp(ln_theta(j))*x~^(alpha-1)+1-delta) c c where c c~ = c(x~,j), c x~ = exp(ln_theta(i))*x^alpha+(1-delta)*x-c(x,i) c pi(i,j)= probability of transiting from state i to j c c(0,i) = 0, all i c x in [0,xmax] c i,j in {1,2,...,I} c c ALGORITHM: the finite element method with piecewise bilinear c shape functions N_a(x), ie., c(x,i) = sum N_a(x)*c^i_a, c where c^i_a's are the coefficients computed here. c c TEST CASE: delta=1, tau=1 ==> c = (1-beta*alpha)*theta*k^alpha. c c INPUT FILE: (a) sgm1dl.inp: the following need to be specified c c Input Description Format c ------------------------------------------------------ c alpha capital share real in (0,1) c beta discount factor real in (0,1) c delta depreciation rate real in [0,1] c tau utility parameter real >= 0 c gamma penalty parameter real >= 0 c ievenx option for evenly 1 (=yes) or 0 c spaced x-partition c initialize option to provide 1 (=yes) or 0 c initial guess for c c nbc number of boundary integer c conditions c x grid for capital if ievenx=1, c 2 x 1 real vector c with x(1),x(nx) c otherwise, c nx x 1 real vector c with entire grid c ln(theta) grid for technology ns x 1 real vector c pi transition matrix ns*ns x 1 real vector c (column by column) c c initial guess for nx*ns x 1 real c consumption function vector; read c(:,1) c (if initialize=1) first, c(:,2) c second, etc. c i,j,c points c(x(i),j) triplets of (integer, c that are to be fixed integer,real) c during computation c c EXAMPLE: suppose the mesh is given by c c c e=1 e=2 e=3 e=4 < element #s c *-------*--------*------------*-----* c 1 2 3 4 5 < node #s c x(1) x(2) x(3) x(4) x(5) < x values c c For this example, the length of the x partition (nx) c is 5 and the number of elements (ne) is 4. c c c Notes: c i. see sgm1dl.inp for an example input file; c ii. if grid size changes, edit this file and c change dimensions nx and ns; c iii. all other inputs are set in this file; c see comment lines below. c c OUTPUT FILES: (b) sgm1dl.nxt: new input file with solution for c(k,i) c (c) sgm1dl.dat: input file for Matlab post-processing c routine -- type sgm1dl in Matlab. c c REQUIRED (BLAS) dgemm,dtrsm,idamax,dger,dscal,dswap c PROGRAMS: (LAPACK) dgesv c (PRESS) qgausl c c REFERENCES: Eischen (notes) c Lapack and BLAS, ftp to netlib2.cs.utk.edu c McGrattan, Federal Reserve Bank of Minneapolis, SR#164 c Press, et. al. (1986), NUMERICAL RECIPES c Taylor and Uhlig (1990), JBES c c SEE ALSO: sgm1dq.f, sgm2dl.f, sgm2dq.f c c Ellen McGrattan, 2-20-91 c Revised, 9-6-94, ERM c * c c Set dimensions: c nx = length of x partition c ns = length of theta c nq = maximum number of quadrature points per element (in x c direction); nq>=nxp(i), i=1:ne. c parameter (nx=10,ns=1, nq=10) implicit double precision (a-h,o-z) double precision f(nx*ns), df(nx*ns,nx*ns), work(nx*ns), c xa(nx), ax(nq,nx-1), wx(nq,nx-1), theta(ns), c pi(ns,ns), c0(nx,ns), c1(nx,ns) integer iwork(nx*ns), nxp(nx-1), ibc(nx*ns) c nxm1 = nx-1 ne = nxm1 na = nx*ns c c Open input and output files. c open(unit=5, file='sgm1dl.inp') open(unit=7, file='sgm1dl.nxt') open(unit=8, file='sgm1dl.dat') c c Read in parameters, grid, boundary conditions, and if specified, c the initial guess for consumption. c read(5,*) alpha read(5,*) beta read(5,*) delta read(5,*) tau read(5,*) gamma read(5,*) ievenx read(5,*) initialize read(5,*) nbc alph1 = alpha-1.d0 delt1 = 1.d0-delta if (ievenx.eq.1) then read(5,*) xa(1) read(5,*) xa(nx) tem = (xa(nx)-xa(1))/float(nxm1) do 10 i=2,nxm1 xa(i) = xa(i-1)+tem 10 continue else do 20 i=1,nx read(5,*) xa(i) 20 continue endif do 30 i=1,ns read(5,*) theta(i) theta(i) = dexp(theta(i)) 30 continue do 40 j=1,ns do 40 i=1,ns read(5,*) pi(i,j) 40 continue if (initialize.eq.1) then do 50 j=1,ns do 50 i=1,nx read(5,*) c0(i,j) ibc(i+(j-1)*nx) = 0 50 continue else c c If initial guess for consumption not given, set c = .5*output. c do 60 j=1,ns do 60 i=1,nx ii = i+(j-1)*nx c0(i,j) = 0.5*(theta(j)*xa(i)**alpha+delt1*xa(i)) ibc(ii) = 0 60 continue endif c c Impose boundary conditions. c do 70 i=1,nbc read(5,*) ii,jj,tem ibc(ii+(jj-1)*nx) = 1 c0(ii,jj) = tem 70 continue c c Other parameters that are set inside the program c * nxp(i), i=1,nx are the # of quadrature points per x interval c * stopping criteria c * maximum number of interations c do 80 i=1,nxm1 nxp(i) = 2 80 continue crit = 0.0000001d0 maxit = 20 c c Compute quadrature abscissas and weights. c do 110 i=1,nxm1 call qgausl (nxp(i), xa(i), xa(i+1), ax(1,i), wx(1,i)) 110 continue c c The Newton-Raphson iteration, i.e., find x such that f(x)=0 where c x are coefficients of the finite element approximation. c do 130 it=1,maxit c c Initialize f and its derivative df. c do 140 i=1,na f(i) = 0.d0 do 140 j=1,na df(i,j) = 0.d0 140 continue c c Construct f(c0) = int R(x,y;c0) N_n(x,y) dx dy and its derivative on c each element n. c do 150 j=1,ns do 160 n=1,ne c c Element n is the interval [xa(n),xa(n+1)] c x1 = xa(n) x2 = xa(n+1) c c Compute f at all quadrature points on the element. c do 170 i=1,nxp(n) x = ax(i,n) u = wx(i,n) xi = 2.d0*(x-x1)/(x2-x1)-1.d0 bs1 = 0.5d0*(1.d0-xi) bs2 = 0.5d0*(1.d0+xi) c c Compute consumption using finite element approximation. c c = c0(n,j)*bs1+c0(n+1,j)*bs2 c c Compute the derivatives of the penalty function (s,ds). c y = theta(j)*x**alpha d = dmin1(y-c,0.d0) s = gamma*d*d ds = -2.d0*gamma*d c c Find next period's capital stock on grid. c xt = y+delt1*x-c ii = 1 do 180 k=1,nxm1 if (xt.ge.xa(k)) ii = k 180 continue x1t = xa(ii) x2t = xa(ii+1) dxt = x2t-x1t xi = 2.d0*(xt-x1t)/dxt-1.d0 bs1t = 0.5d0*(1.d0-xi) bs2t = 0.5d0*(1.d0+xi) dftxt = alpha*xt**alph1 c c Initialize variables for summation. c sum1 = 0.d0 dsum1 = 0.d0 do 190 k=1,na work(k) = 0.d0 190 continue c c The following loop computes the sum in the Euler equation: c c beta sum pi(j,k) c~^{-tau}(alpha*exp(ln_theta(k)) c *x~^(alpha-1)+1-delta) c do 200 k=1,ns c c Given next period's state (xt,ln_theta(k)), compute the c following variables: consumption (ct), marginal utility c (dutct), the derivative of consumption with respect to c capital (dctxt), the derivative of output with respect c to capital (dftkt), derivatives of the penalty function c next period (st,dst), and the integral of the Euler c equation (sum1). c ct = c0(ii,k)*bs1t+c0(ii+1,k)*bs2t dutct = beta*ct**(-tau) dctxt = -(c0(ii,k)-c0(ii+1,k))/dxt dftkt = theta(k)*dftxt+delt1 d = dmin1(theta(k)*xt**alpha-ct,0.d0) st = gamma*d*d dst = -2.d0*gamma*d sum1 = sum1 + pi(j,k)*(dutct*dftkt-beta*delt1*st) c c Compute the derivatives of the integral with respect to c the coefficients of current consumption (dsum1) and the c coefficients of next period consumption (dsum2). c dsum1 = dsum1 + pi(j,k)*(dutct*(tau*dctxt*dftkt/ c ct-alph1*theta(k)*dftxt/xt)-beta*delt1* c dst*(theta(k)*dftxt-dctxt)) dsum2 = pi(j,k)*(-tau*dutct*dftkt/ct-beta*delt1* c dst) m = ii+(k-1)*nx work(m) = work(m) + dsum2*bs1t work(m+1) = work(m+1) + dsum2*bs2t 200 continue c c The vector `work' contains derivatives of the residual. c dres1 = dsum1+tau*c**(-tau-1.d0)+ds m = n+(j-1)*nx work(m) = work(m) + dres1*bs1 work(m+1) = work(m+1) + dres1*bs2 c c Add the residual x quadrature weights x shape functions to f(.). c res = (sum1-c**(-tau)+s)*u f(m) = f(m) + res*bs1 f(m+1) = f(m+1) + res*bs2 c c Add derivative of the residual x quadrature weights x shape c functions to df(.). c do 220 k=1,na df(m,k) = df(m,k) + work(k)*u*bs1 df(m+1,k) = df(m+1,k) + work(k)*u*bs2 220 continue 170 continue 160 continue 150 continue c c Invert the derivative of f (i.e., df). c do 111 j=1,na do 111 i=1,na write(55,*) df(i,j) 111 continue ii=0 do 230 i=1,na if (ibc(i).eq.0) then ii = ii+1 f(ii) = f(i) do 240 j=1,na df(ii,j) = df(i,j) 240 continue do 250 j=1,na df(j,ii) = df(j,i) 250 continue endif 230 continue call dgesv(na-nbc,1,df,na,iwork,f,na,info) c c Issue warnings if necessary. c if (info.lt.0) write(*,*) 'WARNING: illegal value for', c ' argument ',-info,' in DGESV' if (info.gt.0) write(*,*) 'WARNING: singular matrix', c ' when solving df*x=f in SGM1DL' c c Do Newton update. c sum1 = 0.d0 jj = 0 do 260 j=1,ns do 260 i=1,nx ii = (j-1)*nx+i if (ibc(ii).eq.0) then jj = jj+1 c c1(i,j) = c0(i,j) - f(ii) c sum1 = sum1 + f(ii)*f(ii) c1(i,j) = c0(i,j) - f(jj) sum1 = sum1 + f(jj)*f(jj) else c1(i,j) = c0(i,j) endif 260 continue sum1 = dsqrt(sum1)/float(na-nbc) c c Check to see if the solution is converged. If so, stop. c write(*,*) sum1 if (sum1.lt.crit) go to 999 do 270 j=1,ns do 270 i=1,nx c0(i,j) = c1(i,j) 270 continue 130 continue 999 continue c c Write new input file with results. c write(7,'(9X,F12.6,12X,''alpha in (0,1) '')') alpha write(7,'(9X,F12.6,12X,''beta in (0,1) '')') beta write(7,'(9X,F12.6,12X,''delta in [0,1] '')') delta write(7,'(9X,F12.6,12X,''tau >= 0 '')') tau write(7,'(3X,F18.6,12X,''gamma >= 0 '')') gamma write(7,'(9X,I12,12X,''evenly-spaced grid (1=yes,0=no)'')') ievenx write(7,'(20X,''1'',12X,''read initial guess for c (1=yes)'')') write(7,'(9X,I12,12X,''number of boundary constraints'')') nbc write(7,*) write(7,'(9X,F12.6,12X,''lower bound for capital grid'')') xa(1) if (ievenx.eq.0) then do 280 i=2,nxm1 write(7,1) xa(i) 280 continue endif write(7,'(9X,F12.6,12X,''upper bound for capital grid'')') xa(nx) write(7,*) write(7,'(9X,F12.6,12X, c ''grid for log of the technology shock'')') dlog(theta(1)) do 281 i=2,ns write(7,1) dlog(theta(i)) 281 continue write(7,*) write(7,'(9X,F12.6,12X, c ''transition probabilities column by column'')') pi(1,1) do 290 i=2,ns write(7,1) pi(i,1) 290 continue do 292 j=2,ns do 292 i=1,ns write(7,1) pi(i,j) 292 continue write(7,*) write(7,'(9X,F12.6,12X,''consumption function '')') c1(1,1) do 300 i=2,nx write(7,1) c1(i,1) 300 continue do 305 j=2,ns do 305 i=1,nx write(7,1) c1(i,j) 305 continue write(7,*) ii = 0 do 310 j=1,ns do 310 i=1,nx if (ibc(i+(j-1)*nx).eq.1) then if (ii.eq.0) then write(7,'(1X,I3,1X,I3,1X,F12.6,12X, c ''boundary conditions '')') i,j,c1(i,j) ii = 1 else write(7,2) i,j,c1(i,j) endif endif 310 continue c c Write parameters and solution to data file for Matlab routine. c write(8,1) alpha write(8,1) beta write(8,1) delta write(8,1) tau write(8,*) gamma write(8,*) nx do 320 i=1,nx write(8,1) xa(i) 320 continue write(8,*) ns do 330 i=1,ns write(8,1) dlog(theta(i)) 330 continue do 340 j=1,ns do 340 i=1,ns write(8,1) pi(i,j) 340 continue do 350 j=1,ns do 350 i=1,nx write(8,1) c1(i,j) 350 continue write(8,*) it c c Format statements. c 1 format (9X,F12.6) 2 format (1X,I3,1X,I3,1X,F12.6) stop end