c TITLE: implementation of finite-element method for stochastic c growth model, 2-D c c PROBLEM: find c=c(k,theta) that solves the functional equation: c u'(c) = beta int u'(c~)F'(k~,theta~)g(eps)deps c c where c c~ = c(k~,theta~), c k~ = F(k,theta)-c(k,theta), c theta~ = h(theta,epsilon), c c(0,theta) = 0 c c ALGORITHM: the finite element method with piecewise bilinear shape c functions N(x,y) and c(x,y) = sum Na(x,y)*ca, ca constant c a c MODEL: the discrete-time nonstochastic growth model with c agents choosing c(t) to maximize c c E[ sum beta^t u(c[t]) | k_0] c t c subj. to k[t+1] = F(k[t],theta[t]) - c[t] c c FUNCTIONS: u(c) = c^(1-tau)/(1-tau) c F(k,theta) = theta*k^alpha + (1-delta)*k c h(theta,e) = theta^rho*exp(e) c steady state for k: c ((1-beta*(1-delta))/(beta*alpha))^(1/(alpha-1)) c c TEST CASE: delta=1, tau=1 c ==> c(k,theta) = (1-beta*alpha)*theta*k^alpha c c INPUTS: (a) qr2.inp: alpha,beta,delta,rho,sige,tau c c0(i),ibc(i) i=1,na c xa(i),nxp(i) i=1,nx c xa(nx+1) c ya(i),nyp(i) i=1,ny c ya(ny+1) c where c0's are initial guesses for ci's c ibc's are 0 or 1 if pt. is essential boundary cond. c xa's,ya's are grid pts. c nzp(i) are # of points per interval for z c OUTPUTS: (b) qr2.out: general results c (c) qr2.nxt: new input file c (d) qr2.dat: input file for Matlab c c REFERENCE: Jeff Eishen (notes) c Taylor and Uhlig (1990), JBES c c SEE ALSO: grow2db.f in ~erm/fem/grow directory c c Ellen McGrattan, 2-20-91 c Revised, 10-22-93, ERM c parameter (na=189,nx=26,ny=6,nxm=5,nym=5,nep=10,nt=1000000) implicit double precision (a-g,o-z) double precision f(na), df(na,na), ft(na), dft(na,na), c work(nt), rcond, xa(nx+1), ya(ny+1), xx(nxm,nx), c wx(nxm,nx), yy(nym,ny), wy(nym,ny), x(nxm), u(nxm), c y(nym), v(nym), yyt(nep), ep(nep), out(na), c we(nep), wgt(nep), c0(na), c1(na), wrk(nep,4), c param(7),f1(na),f2(na),cp(na),cm(na),del(na) integer iwork(nt), nxp(nx), nyp(ny), nb(na), c jj(nep), kk(nep), ibc(na) c open(unit=5, file='qr2b.inp') open(unit=7, file='qr2b.out') open(unit=8, file='qr2b.nxt') open(unit=9, file='qr2b.dat') pi = 3.141592654d0 nx1 = nx+1 ny1 = ny+1 read(5,*) alpha,beta,delta,rho,sige,tau,gam do 10 i=1,na read(5,*) c0(i),ibc(i) 10 continue n = 1 do 20 i=1,na if (ibc(i).eq.0) then nb(n) = i n = n + 1 endif 20 continue nnb = n - 1 do 30 i=1,nx read(5,*) xa(i),nxp(i) 30 continue read(5,*) xa(nx1) do 40 i=1,ny read(5,*) ya(i),nyp(i) 40 continue read(5,*) ya(ny1) c alph1 = alpha-1.d0 delt1 = 1.d0-delta call qgausl (nep, -2.88d0*sige, 2.88d0*sige, ep, we) tem1 = 1.d0/(2.d0*sige*sige) tem2 = 1.d0/dsqrt(pi/tem1) sum1 = 0.d0 do 50 i=1,nep wgt(i) = dexp(-ep(i)*ep(i)*tem1)*tem2*we(i) sum1 = sum1 + wgt(i) 50 continue do 60 i=1,nep wgt(i) = wgt(i)/sum1 60 continue do 80 i=1,nx call qgausl (nxp(i), xa(i), xa(i+1), x, u) do 80 j=1,nxp(i) xx(j,i) = x(j) wx(j,i) = u(j) 80 continue do 90 i=1,ny call qgausl (nyp(i), ya(i), ya(i+1), y, v) do 90 j=1,nyp(i) yy(j,i) = y(j) wy(j,i) = v(j) 90 continue ne=nx*ny crit = 0.0000001d0 do 95 j = 1,ny1 do 95 i = 1,nx1 ii = ii+1 out(ii) = ya(j)*xa(i)**alpha 95 continue maxit = 100 do 100 it=1,maxit do 110 i=1,na f(i) = 0.d0 do 110 j=1,na df(i,j) = 0.d0 110 continue do 150 n=1,ne ix = mod(n-1,nx)+1 iy = int((n-1)/nx)+1 n1 = nxp(ix) n2 = nyp(iy) x1 = xa(ix) x2 = xa(ix+1) dx = x2-x1 y1 = ya(iy) y2 = ya(iy+1) dy = y2-y1 nd1= ix+nx1*(iy-1) nd2= nd1+1 nd4= nd1+nx1 nd3= nd4+1 do 160 i=1,n1 x(i) = xx(i,ix) u(i) = wx(i,ix) 160 continue do 170 j=1,n2 y(j) = yy(j,iy) v(j) = wy(j,iy) yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 170 i=1,n1 xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx c = c0(nd1)*xp*yp+c0(nd2)*xm*yp+ c c0(nd3)*xm*ym+c0(nd4)*xp*ym d = c-y(j)*x(i)**alpha c term in objective: -gam*(d*d*d+abs(d*d*d)) s = gam*3.d0*d*(d+dsign(d,0.d0)) dsdc = gam*6.d0*(d+dsign(d,0.d0)) xt = x(i)**alpha*y(j)-c+delt1*x(i) dftxt = alpha*xt**alph1 ii = 1 do 180 m=1,nx j1 = int( dmax1(xt-xa(m),0.d0)/(xt-xa(m))+.2 ) ii = ii + j1 180 continue ii = max(1,ii-1) do 190 l=1,nep yyt(l) = y(j)**rho*dexp(ep(l)) jj(l) = 1 do 200 m=1,ny j1 = int( dmax1(yyt(l)-ya(m),0.d0)/ c (yyt(l)-ya(m))+.2) jj(l) = jj(l) + j1 200 continue jj(l) = max(1,jj(l)-1) 190 continue sum1 = 0.d0 sum2 = 0.d0 sum3 = 0.d0 sum4 = 0.d0 sum5 = 0.d0 do 220 l=1,na work(l) = 0.d0 220 continue do 230 l=1,nep yt = yyt(l) x1t = xa(ii) x2t = xa(ii+1) dxt = x2t-x1t y1t = ya(jj(l)) y2t = ya(jj(l)+1) dyt = y2t-y1t xpt = (x2t-xt)/dxt xmt = (xt-x1t)/dxt ypt = (y2t-yt)/dyt ymt = (yt-y1t)/dyt nd1t = (jj(l)-1)*nx1+ii nd2t = nd1t+1 nd4t = nd1t+nx1 nd3t = nd4t+1 kk(l) = nd1t c c compute c * next period consumption: ct c * beta*du(c[t+1])/dc[t+1]: dutct c * dc[t+1]/dk[t+1]: dctxt c * dF[t+1]/dk[t+1]: dftkt c * int u'(c~)*F'(k~)f(e)de: sum1 c * d(sum1)/d(c_a): dsum c * d(sum1)/d(c_b): dsum (2nd instance) c where c_a,c_b are unknowns in c,ct respectively c ct = c0(nd1t)*xpt*ypt+c0(nd2t)*xmt*ypt+ c c0(nd3t)*xmt*ymt+c0(nd4t)*xpt*ymt dt = ct-yt*xt**alpha st = gam*3.d0*dt*(dt+dsign(dt,0.d0)) dstdc = gam*6.d0*(dt+dsign(dt,0.d0)) dutct = beta*ct**(-tau) dctxt = (-(c0(nd1t)-c0(nd2t))*ypt c +(c0(nd3t)-c0(nd4t))*ymt)/dxt dftkt = yt*dftxt+delt1 sum1 = sum1 + wgt(l)*(dutct*dftkt-beta*delt1*st) dsum = wgt(l)*dutct*(tau*dctxt*dftkt/ct c -alph1*yt*dftxt/xt)+ c wgt(l)*beta*delt1*dstdc*(dctxt-yt*dftxt) sum2 = sum2 + dsum*xp*yp sum3 = sum3 + dsum*xm*yp sum4 = sum4 + dsum*xm*ym sum5 = sum5 + dsum*xp*ym dsum = -wgt(l)*(tau*dutct*dftkt/ct+ c beta*delt1*dstdc) wrk(l,1) = dsum*xpt*ypt wrk(l,2) = dsum*xmt*ypt wrk(l,3) = dsum*xmt*ymt wrk(l,4) = dsum*xpt*ymt 230 continue tem1 = tau*c**(-tau-1.d0)+dsdc work(nd1) = sum2 + tem1*xp*yp work(nd2) = sum3 + tem1*xm*yp work(nd3) = sum4 + tem1*xm*ym work(nd4) = sum5 + tem1*xp*ym do 232 l = 1,nep nd1t = kk(l) nd2t = nd1t+1 nd4t = nd1t+nx1 nd3t = nd4t+1 work(nd1t) = work(nd1t) + wrk(l,1) work(nd2t) = work(nd2t) + wrk(l,2) work(nd3t) = work(nd3t) + wrk(l,3) work(nd4t) = work(nd4t) + wrk(l,4) 232 continue res = sum1-c**(-tau)+s tem1 = u(i)*v(j)*res f(nd1) = f(nd1) + tem1*xp*yp f(nd2) = f(nd2) + tem1*xm*yp f(nd3) = f(nd3) + tem1*xm*ym f(nd4) = f(nd4) + tem1*xp*ym c c tem1 = u(i)*v(j) c do 250 l=1,na c tem2 = tem1*work(l) c df(nd1,l) = df(nd1,l) + tem2*xp*yp c df(nd2,l) = df(nd2,l) + tem2*xm*yp c df(nd3,l) = df(nd3,l) + tem2*xm*ym c df(nd4,l) = df(nd4,l) + tem2*xp*ym c250 continue 170 continue 150 continue param(1) = alpha param(2) = beta param(3) = delta param(4) = rho param(5) = sige param(6) = tau param(7) = gam do 372 kt = 1,na del(kt) = dabs(c0(kt)*0.0001d0) do 373 mm=1,na cp(mm) = c0(mm) cm(mm) = c0(mm) 373 continue cp(kt) = cp(kt)+del(kt) cm(kt) = cm(kt)-del(kt) call qr2der(cp,param,xa,nxp,ya,nyp,f1) call qr2der(cm,param,xa,nxp,ya,nyp,f2) do 374 mt=1,na df(mt,kt) = .5d0*(f1(mt)-f2(mt))/del(kt) 374 continue 372 continue c do 444 j=1,na c write(77,801) f(j) c do 444 i=1,na c write(78,801) df(i,j) c444 continue c stop c ii=0 do 260 i=1,na if (ibc(i).eq.0) then ii = ii+1 f(ii) = f(i) do 270 j=1,na df(ii,j) = df(i,j) 270 continue do 280 j=1,na df(j,ii) = df(j,i) 280 continue endif 260 continue c call dgeco(df, na, nnb, iwork, rcond, work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when solving df*x=f' call dgesl(df, na, nnb, iwork, f, 0) sum1 = 0.d0 do 300 i=1,nnb c1(nb(i)) = c0(nb(i)) - f(i) sum1 = sum1+f(i)*f(i) 300 continue sum1 = dsqrt(sum1)/float(nnb) write(*,*) ' at iteration ', it,' the residual is ', sum1 if (sum1.lt.crit) then go to 999 endif do 310 i=1,nnb c0(nb(i)) = c1(nb(i)) 310 continue 100 continue 999 continue c c print results c write(7,*) 'GROW-2D-B RESULTS' write(7,*) write(7,*) ' Parameters ' write(7,'('' beta = '',E9.3,'' '')') beta write(7,'('' delta = '',E9.3,'' '')') delta write(7,'('' tau = '',E9.3,'' '')') tau write(7,'('' rho = '',E9.3,'' '')') rho write(7,'('' sige = '',E9.3,'' '')') sige write(7,'('' alpha = '',E9.3,'' '')') alpha write(7,*) write(7,'('' Steady State for capital = '',E9.3,'' '')') c ((1.d0-beta*delt1)/(beta*alpha))**(1.d0/(alpha-1.d0)) c write(7,*) write(7,*) ' Coefficients for Shape functions ' do 320 i=1,na write(7,802) c0(i),out(i)-c0(i) 320 continue write(7,*) write(7,*) ' Number of iterations = ', it write(8,806) alpha,beta,delta,rho,sige,tau,gam n = 1 do 340 i=1,na iwork(i) = 1 340 continue do 350 i=nnb,1,-1 iwork(nb(i)) = 0 350 continue do 360 i=1,na write(8,804) c0(i),iwork(i) 360 continue do 370 i=1,nx write(8,804) xa(i),nxp(i) 370 continue write(8,801) xa(nx1) do 380 i=1,ny write(8,804) ya(i),nyp(i) 380 continue write(8,801) ya(ny1) c 800 format (1X, I3, F12.6) 801 format (4X, E12.6) 802 format (2(4X, E12.6)) 803 format (3(4X, E12.6)) 804 format (4X, E12.6, 4X, I3) 805 format (5(2X, E8.3)) 806 format (7(1X, E10.3)) stop end