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 + terms for penalty function 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 <= F(k,theta)-(1-delta)*k 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]) -gam*penalty(i[t]<0) | 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 penalty = (c-theta*k^alpha)^3+abs((c-theta*k^alpha)^3) c steady state for k (without penalty imposed): 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) trial.inp: alpha,beta,delta,rho,sige,tau,gam 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) trial.out: general results c (c) trial.nxt: new input file 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, 6-09-93, ERM c parameter (na=72,nx=23,ny=2,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) integer iwork(nt), nxp(nx), nyp(ny), nb(na), c jj(nep), kk(nep), ibc(na) c open(unit=5, file='trial.inp') open(unit=7, file='trial.out') open(unit=8, file='trial.nxt') 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.00000001d0 maxit = 100 ii = 0 do 95 j = 1,ny1 do 95 i = 1,nx1 ii = ii+1 out(ii) = ya(j)*xa(i)**alpha 95 continue 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)-.22) 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 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 itask = 1 call dgefs (df, na, nnb, f, itask, ind, work, iwork, rcond) 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,*) 'RESULTS FOR TRIAL.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,'('' gam = '',E9.3,'' '')') gam write(7,*) write(7,'('' Steady State for capital = '',E9.3,'' '')') c ((1.d0-beta*delt1)/(beta*alpha))**(1.d0/(alpha-1.d0)) write(7,*) write(7,*) ' Results from dgefs.f ' if (ind .eq. -10) then write (7,*) ' error code =', ind else if (ind .lt. 0) then write (7,*) ' error code =', ind stop else write (7,*) ' number of accurate digits =', ind end if c write(7,*) write(7,*) ' Coefficients for Shape functions', c ' (consumption,investment)' 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