c PROGRAM: SGM2DL.F c c DESCRIPTION: implementation of finite-element method for stochastic c growth model, 2-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 theta[t+1] = theta[t]**rho*exp(eps[t]) c eps[t] ~ N(0,sige^2) c c[t],k[t] >= 0 c i[t] >= 0 if penalty parameter gamma is c sufficiently large c c PROBLEM: find c=c(x,y) that solves the functional equation: c c c^{-tau} = beta int c~^{-tau}(alpha*exp(tanh^{-1}(y~)) c *x~^(alpha-1)+1-delta) exp(-omega^2)/sqrt(pi) domega c c where c c~ = c(x~,y~), c x~ = exp(tanh^{-1}(y))*x^alpha+(1-delta)*x-c(x,y) c y~ = tanh(rho*tanh^{-1}(y)+sqrt(2)*sige*omega) c c(0,y) = 0 c y in (-1,1) c x in [0,xmax] c c NOTE: y=tanh(ln(theta)), where theta is defined above. c c ALGORITHM: the finite element method with piecewise bilinear c shape functions N_a(x,y), ie., c(x,y) = sum N_a(x,y)*ca, c where ca'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) sgm2dl.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 rho persistence of theta real in (-1,1) c sige standard deviation real >= 0 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 ieveny option for evenly 1 (=yes) or 0 c spaced y-partition c initialize option to provide 1 (=yes) or 0 c initial guess for c c nbc number of boundary integer c constraints c x grid for capital nx x 1 real vector c w/ x(1),x(2),..x(nx) c or 2 x 1 with x(1) c and x(nx) if c ievenx = 1 c y grid for technology, ny x 1 real vector c eg. tanh(log(theta)) w/ y(1),y(2),..y(ny) c or 2 x 1 with y(1) c and y(ny) if c ieveny = 1 c c initial guess for na x 1 real vector c consumption function where na=nx*ny; ith c (if initialize=1) element of c is c function for node i c (see EXAMPLE below) c a,c_a boundary conditions if nbc>0, reads c (for some in nbc x 2 matrix c nodes a) with integer (node c #) in column 1 and c real (constrained c consumption value) c in column 2 c c EXAMPLE: suppose the mesh is given by c c y-grid c ------ c y(4) *--------*--------*--------*--------* c |16 |17 |18 |19 |20 c | | | | | c | | | | | c |e=9 |e=10 |e=11 |e=12 | c y(3) *--------*--------*--------*--------* c |11 |12 |13 |14 |15 c | | | | | c | | | | | c |e=5 |e=6 |e=7 |e=8 | c y(2) *--------*--------*--------*--------* c |6 |7 |8 |9 |10 c | | | | | c | | | | | c |e=1 |e=2 |e=3 |e=4 | c y(1) *--------*--------*--------*--------* c 1 2 3 4 5 c c x-grid: x(1) x(2) x(3) x(4) x(5) c c c where *'s mark nodes and e's are element numbers. c The *'s and e's must be marked in a particular c way: first mark in the x-direction and then in c the y-direction until e=12 (=ne) and a=20 (=na) c where a is the node number. The smallest x and c y values are at node 1. For this example, the c number of points in the x partition (nx) is 5; c the number of points in the y partition (ny) is 4. c c c Notes: c i. y=tanh(ln(theta)) -- remember to convert c theta values for input file; c ii. see sgm2dl.inp for an example input file; c iii. if grid size changes, edit this file and c change dimensions nx and ny; c iv. all other inputs are set in this file; c see comment lines below. c c OUTPUT FILES: (b) sgm2dl.nxt: new input file with solution for c(x,y) c (c) sgm2dl.dat: input file for Matlab post-processing c routine -- run sgm2dl in Matlab. c c REQUIRED (BLAS) daxpy, dcopy, ddot, dgemm, dgemv, c PROGRAMS: dger, dnrm2, drot, dscal, dswap, c dtrmm, dtrmv, dtrsm, idamax c (LAPACK) dgesv c (PRESS, ET.AL.) qgaush, 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: sgm1dl.f, sgm1dq.f, sgm2dq.f c c Ellen McGrattan, 2-20-91 c Revised, 1-25-96, ERM c * c c Set dimensions: c nx = length of xa c ny = length of ya c nep = number of quadrature points for integral in Euler equation c nq = maximum number of quadrature points per element (in x or c y direction); nq>=nxp(i) and nq>=nyp(i) all i. c parameter (nx=7,ny=4,nep=10, nq=5) implicit double precision (a-h,o-z) double precision f(nx*ny), df(nx*ny,nx*ny), work(nx*ny), c xa(nx), ya(ny), ax(nq,nx-1), wx(nq,nx-1), ay(nq,ny-1), c wy(nq,ny-1), omega(nep), we(nep), c0(nx*ny), c1(nx*ny) integer iwork(nx*ny), nxp(nx-1), nyp(ny-1), ibc(nx*ny) c nxm1 = nx-1 nym1 = ny-1 ne = nxm1*nym1 na = nx*ny c c Open input and output files. c open(unit=5, file='sgm2dl.inp') open(unit=7, file='sgm2dl.nxt') open(unit=8, file='sgm2dl.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,*) rho read(5,*) sige read(5,*) tau read(5,*) gamma read(5,*) ievenx read(5,*) ieveny 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 if (ieveny.eq.1) then read(5,*) ya(1) read(5,*) ya(ny) tem = (ya(ny)-ya(1))/float(nym1) do 30 i=2,nym1 ya(i) = ya(i-1)+tem 30 continue else do 40 i=1,ny read(5,*) ya(i) 40 continue endif if (initialize.eq.1) then do 50 i=1,na read(5,*) c0(i) ibc(i) = 0 50 continue else c c If initial guess for consumption not given, set c = .5*output. c do 60 j=1,ny do 60 i=1,nx ii = (j-1)*nx+i tech = dsqrt((1.d0+ya(j))/(1.d0-ya(j))) c0(ii) = 0.5*(tech*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,tem ibc(ii) = 1 c0(ii) = 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 * nyp(i), i=1,ny are the # of quadrature points per y interval c * stopping criteria c * maximum number of interations. c do 80 i=1,nxm1 nxp(i) = 3 80 continue do 90 i=1,nym1 nyp(i) = 3 90 continue crit = 0.0000001d0 maxit = 100 pi = 3.141592654d0 sq2 = dsqrt(2.d0) sqpii = 1.d0/dsqrt(pi) c c Compute quadrature abscissas and weights. c call qgaush (nep, omega, we) sum1 = 0.d0 do 100 i=1,nep we(i) = we(i)*sqpii sum1 = sum1 + we(i) 100 continue do 110 i=1,nep we(i) = we(i)/sum1 110 continue do 120 i=1,nxm1 call qgausl (nxp(i), xa(i), xa(i+1), ax(1,i), wx(1,i)) 120 continue do 130 i=1,nym1 call qgausl (nyp(i), ya(i), ya(i+1), ay(1,i), wy(1,i)) 130 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 140 it=1,maxit c c Initialize f and its derivative df. c do 150 i=1,na f(i) = 0.d0 do 150 j=1,na df(i,j) = 0.d0 150 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 160 n=1,ne c c Element n is the rectangle [xa(ix),xa(ix+1)] x [ya(iy),ya(iy+1)]. c ix = mod(n-1,nxm1)+1 iy = (n-1)/nxm1+1 x1 = xa(ix) x2 = xa(ix+1) dx = x2-x1 y1 = ya(iy) y2 = ya(iy+1) dy = y2-y1 c c Assign numbers to the nodes of element n. c nd1 = ix+nx*(iy-1) nd2 = nd1+1 nd4 = nd1+nx nd3 = nd4+1 c c Compute f at all quadrature points on the element. c do 170 j=1,nyp(iy) y = ay(j,iy) v = wy(j,iy) do 170 i=1,nxp(ix) x = ax(i,ix) u = wx(i,ix) xi = 2.d0*(x-x1)/(x2-x1)-1.d0 eta = 2.d0*(y-y1)/(y2-y1)-1.d0 c c Assign numbers to the basis functions defined on element n. c bs1 = 0.25d0*(1.d0-xi)*(1.d0-eta) bs2 = 0.25d0*(1.d0+xi)*(1.d0-eta) bs3 = 0.25d0*(1.d0+xi)*(1.d0+eta) bs4 = 0.25d0*(1.d0-xi)*(1.d0+eta) c c Compute consumption using finite element approximation. c c = c0(nd1)*bs1+c0(nd2)*bs2+c0(nd3)*bs3+c0(nd4)*bs4 c c Compute the derivatives of the penalty function (s,ds). c tech = 0.5d0*(dlog(1.d0+y)-dlog(1.d0-y)) gnp = dexp(tech)*x**alpha d = dmin1(gnp-c,0.d0) s = gamma*d*d ds = -2.d0*gamma*d c c Find next period's capital stock on grid. c xt = gnp+delt1*x-c dftxt = alpha*xt**(alpha-1.d0) ii = 1 do 180 k=1,nxm1 if (xt.ge.xa(k)) ii = k 180 continue x1t = xa(ii) x2t = xa(ii+1) xi = 2.d0*(xt-x1t)/(x2t-x1t)-1.d0 dxi = 2.d0/(x2t-x1t) 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 integral: c c beta int c~^{-tau}(alpha*exp(tanh^{-1}(y~)) c *x~^(alpha-1)+1-delta) exp(-eta^2)/sqrt(pi) deta. c do 200 k=1,nep c c Compute next period technology (yt). c yt = dtanh(rho*tech+sq2*sige*omega(k)) jj = 1 do 210 l=1,nym1 if (yt.ge.ya(l)) jj = l 210 continue y1t = ya(jj) y2t = ya(jj+1) eta = 2.d0*(yt-y1t)/(y2t-y1t)-1.d0 c c Assign numbers to the nodes of element of (xt,yt). c nd1t = (jj-1)*nx+ii nd2t = nd1t+1 nd4t = nd1t+nx nd3t = nd4t+1 c c Assign numbers to the basis functions defined on element c of (xt,yt) and take derivatives with respect to xt. c bs1t = 0.25d0*(1.d0-xi)*(1.d0-eta) bs2t = 0.25d0*(1.d0+xi)*(1.d0-eta) bs3t = 0.25d0*(1.d0+xi)*(1.d0+eta) bs4t = 0.25d0*(1.d0-xi)*(1.d0+eta) dbs1 =-0.25d0*(1.d0-eta)*dxi dbs3 = 0.25d0*(1.d0+eta)*dxi dbs2 =-dbs1 dbs4 =-dbs3 c c Given next period's state vector (xt,yt), 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), and the integral of the Euler c equation (sum1). c ct = c0(nd1t)*bs1t+c0(nd2t)*bs2t+ c c0(nd3t)*bs3t+c0(nd4t)*bs4t dctxt = c0(nd1t)*dbs1+c0(nd2t)*dbs2+ c c0(nd3t)*dbs3+c0(nd4t)*dbs4 dutct = beta*ct**(-tau) theta = dsqrt((1.d0+yt)/(1.d0-yt)) dftkt = theta*dftxt+delt1 d = dmin1(theta*xt**alpha-ct,0.d0) st = gamma*d*d dst = -2.d0*gamma*d sum1 = sum1 + we(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 + we(k)*(dutct*(tau*dctxt*dftkt/ct- c alph1*theta*dftxt/xt)-beta*delt1* c dst*(theta*dftxt-dctxt)) dsum2 = we(k)*(-tau*dutct*dftkt/ct-beta*delt1*dst) work(nd1t) = work(nd1t) + dsum2*bs1t work(nd2t) = work(nd2t) + dsum2*bs2t work(nd3t) = work(nd3t) + dsum2*bs3t work(nd4t) = work(nd4t) + dsum2*bs4t 200 continue c c The vector `work' contains derivatives of the residual. c dres = dsum1+tau*c**(-tau-1.d0)+ds work(nd1) = work(nd1) + dres*bs1 work(nd2) = work(nd2) + dres*bs2 work(nd3) = work(nd3) + dres*bs3 work(nd4) = work(nd4) + dres*bs4 c c Add the residual x quadrature weights x shape functions to f(.). c res = (sum1-c**(-tau)+s)*u*v f(nd1) = f(nd1) + res*bs1 f(nd2) = f(nd2) + res*bs2 f(nd3) = f(nd3) + res*bs3 f(nd4) = f(nd4) + res*bs4 c c Add derivative of the residual x quadrature weights x shape c functions to df(.). c do 220 k=1,na tem = u*v*work(k) df(nd1,k) = df(nd1,k) + tem*bs1 df(nd2,k) = df(nd2,k) + tem*bs2 df(nd3,k) = df(nd3,k) + tem*bs3 df(nd4,k) = df(nd4,k) + tem*bs4 220 continue 170 continue 160 continue c c Invert the derivative of f (i.e., df). c 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 SGM2DL' c c Do Newton update. c sum1 = 0.d0 ii = 0 do 260 i=1,na if (ibc(i).eq.0) then ii = ii+1 c1(i) = c0(i) - f(ii) sum1 = sum1 + f(ii)*f(ii) else c1(i) = c0(i) endif 260 continue sum1 = dsqrt(sum1)/float(na-nbc) c c Check to see if the solution is converged. If so, stop. c if (sum1.lt.crit) go to 999 do 270 i=1,na c0(i) = c1(i) 270 continue 140 continue 999 continue c c Write new input file with results. c write(7,'(5X,F12.6,12X,''alpha '')') alpha write(7,'(5X,F12.6,12X,''beta '')') beta write(7,'(5X,F12.6,12X,''delta '')') delta write(7,'(5X,F12.6,12X,''rho '')') rho write(7,'(5X,F12.6,12X,''sige '')') sige write(7,'(5X,F12.6,12X,''tau '')') tau write(7,'(5X,F12.6,12X,''gamma '')') gamma write(7,'(5X,I12,12X,''even spacing for x grid (1=yes)'')') ievenx write(7,'(5X,I12,12X,''even spacing for y grid (1=yes)'')') ieveny write(7,'(16X,''1'',12X,''read initial guess for c (1=yes)'')') write(7,'(5X,I12,12X,''number of boundary constraints'')') nbc write(7,*) if (ievenx.eq.1) then write(7,'(5X,F12.6,12X,''lower bound for capital '')') xa(1) write(7,'(5X,F12.6,12X,''upper bound for capital '')') xa(nx) else write(7,'(5X,F12.6,12X,''grid for capital '')') xa(1) do 280 i=2,nx write(7,1) xa(i) 280 continue endif write(7,*) if (ieveny.eq.1) then write(7,'(5X,F12.6,12X, c ''lower bound for tanh(ln(theta))'')') ya(1) write(7,'(5X,F12.6,12X, c ''upper bound for tanh(ln(theta))'')') ya(ny) else write(7,'(5X,F12.6,12X,''grid for tanh(ln(theta))'')') ya(1) do 290 i=2,ny write(7,1) ya(i) 290 continue endif write(7,*) write(7,'(5X,F12.6,12X,''consumption function '')') c1(1) do 300 i=2,na write(7,1) c1(i) 300 continue write(7,*) ii = 0 do 310 i=1,na if (ibc(i).eq.1) then if (ii.eq.0) then write(7,'(1X,I3,1X,F12.6,8X,''boundary conditions '')') c i,c1(i) else write(7,2) i,c1(i) endif ii = 1 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) rho write(8,1) sige write(8,1) tau write(8,*) nx do 320 i=1,nx write(8,1) xa(i) 320 continue write(8,*) ny do 330 i=1,ny write(8,1) ya(i) 330 continue do 340 i=1,na write(8,1) c1(i) 340 continue c c Format statements. c 1 format (5X,F12.6) 2 format (1X,I3,1X,F12.6) stop end