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 c REQUIRED (LINPACK) blas, dgeco, dgefa, dgesl c PROGRAMS: (PRESS, ET.AL.) qgausl 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=18,nx=5,ny=2,nxm=5,nym=5,nep=10,nt=10000) 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), c we(nep), wgt(nep), c0(na), c1(na), wrk(nep,4) integer iwork(nt), nxp(nx), nyp(ny), nb(na), c jj(nep), kk(nep), ibc(na) c open(unit=5, file='qr2.inp') open(unit=7, file='qr2.out') open(unit=8, file='qr2.nxt') pi = 3.141592654d0 nx1 = nx+1 ny1 = ny+1 read(5,*) alpha,beta,delta,rho,sige,tau 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.00001d0 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 = (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 xt = x(i)**alpha*y(j)-c+delt1*x(i) dftxt = alpha*xt**(alpha-1.d0) ii = 1 do 180 m=1,nx j1 = idint( 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 = idint( 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 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 dsum = wgt(l)*dutct*(tau*dctxt*dftkt/ct c -alph1*yt*dftxt/xt) sum2 = sum2 + dsum*xp*yp sum3 = sum3 + dsum*xm*yp sum4 = sum4 + dsum*xm*ym sum5 = sum5 + dsum*xp*ym dsum = -tau*dutct*dftkt*wgt(l)/ct 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) 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) 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 tem1 = u(i)*v(j) do 250 l=1,na tem2 = tem1*work(l) df(nd1,l) = df(nd1,l) + tem2*xp*yp df(nd2,l) = df(nd2,l) + tem2*xm*yp df(nd3,l) = df(nd3,l) + tem2*xm*ym df(nd4,l) = df(nd4,l) + tem2*xp*ym 250 continue 170 continue 150 continue 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) c 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,*) 'RESULTS FOR QR2.F' write(7,*) write(7,*) ' Parameters ' write(7,'('' alpha = '',E9.3,'' '')') alpha write(7,'('' beta = '',E9.3,'' '')') beta write(7,'('' delta = '',E9.3,'' '')') delta write(7,'('' rho = '',E9.3,'' '')') rho write(7,'('' sige = '',E9.3,'' '')') sige write(7,'('' tau = '',E9.3,'' '')') tau 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,801) c0(i) 320 continue write(7,*) write(7,*) ' Number of iterations = ', it write(8,806) alpha,beta,delta,rho,sige,tau 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