c TITLE: implementation of finite-element method for stochastic c growth model, 2-D c c PROBLEM: find c^l=c^l(k) that solves the functional equation: c u'(c^l)=beta sum_{j=1:n} pi(l,j) u'(c~)F'(k~,theta^j) c + terms for penalty function c c l=1,...n, where c c~ = c^j(k~) c k~ = F(k,theta^l)-c^l(k), c c^l(0) = 0 c c^l(k) <= F(k,theta^l)-(1-delta)*k c c ALGORITHM: the finite element method with piecewise bilinear shape c functions N(x) and c(x) = sum Na(x)*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 penalty = (c-theta*k^alpha)^3+abs((c-theta*k^alpha)^3) c theta = n-state Markov chain c steady state for k (without penalty imposed): c ((1-beta*(1-delta))/(beta*alpha*theta))^(1/(alpha-1)) c c TEST CASE: delta=1, tau=1, gam=0 c ==> c(k) = (1-beta*alpha)*theta*k^alpha c c INPUTS: (a) qr2b.inp: alpha,beta,delta,rho,sige,tau,gam c theta(i) i=1,ns c mcp(i,j) i=1,ns j=1,ns c c0(i,j) i=1,na, j=1,ns c xa(i),nxp(i) i=1,nx c xa(nx+1) c where c0's are initial guesses for ci's c xa's are grid pts. c nzp(i) are # of points per interval for z c OUTPUTS: (b) qr2b.out: general results c (c) qr2b.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=21,ns=2,nas=42,nx=20,nxm=5,nt=1000000) implicit double precision (a-g,o-z) double precision f(nas), df(nas,nas), c work(nt), rcond, xa(nx+1), xx(nxm,nx), wx(nxm,nx), c x(nxm), u(nxm), out(na,ns), c0(na,ns), c1(na,ns), c wrk(ns,2), mcp(ns,ns), theta(ns) integer iwork(nt), nxp(nx) c open(unit=5, file='qr2b.inp') open(unit=7, file='qr2b.out') open(unit=8, file='qr2b.nxt') pi = 3.141592654d0 nx1 = nx+1 read(5,*) alpha,beta,delta,tau,gam do 5 i=1,ns read(5,*) theta(i) 5 continue do 10 i=1,ns read(5,*) (mcp(i,j),j=1,ns) 10 continue do 20 i=1,na read(5,*) (c0(i,j),j=1,ns) 20 continue do 30 i=1,nx read(5,*) xa(i),nxp(i) 30 continue read(5,*) xa(nx1) c alph1 = alpha-1.d0 delt1 = 1.d0-delta 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 crit = 0.00000001d0 maxit = 100 do 100 it=1,maxit do 110 i=1,nas f(i) = 0.d0 do 110 j=1,nas df(i,j) = 0.d0 110 continue do 140 l=1,ns do 150 n=1,nx n1 = nxp(n) x1 = xa(n) x2 = xa(n+1) dx = x2-x1 nd1= n nd2= n+1 do 160 i=1,n1 x(i) = xx(i,n) u(i) = wx(i,n) 160 continue do 170 i=1,n1 xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx c = c0(nd1,l)*xp+c0(nd2,l)*xm d = c-theta(l)*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*theta(l)-c+delt1*x(i) dftxt = alpha*xt**alph1 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) sum1 = 0.d0 sum2 = 0.d0 sum3 = 0.d0 do 220 m=1,nas work(m) = 0.d0 220 continue x1t = xa(ii) x2t = xa(ii+1) dxt = x2t-x1t xpt = (x2t-xt)/dxt xmt = (xt-x1t)/dxt nd1t = ii nd2t = ii+1 do 190 m=1,ns ct = c0(nd1t,m)*xpt+c0(nd2t,m)*xmt dt = ct-theta(m)*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,m)-c0(nd2t,m))/dxt dftkt = theta(m)*dftxt+delt1 sum1 = sum1 + mcp(l,m)*(dutct*dftkt-beta*delt1*st) dsum = mcp(l,m)*dutct*(tau*dctxt*dftkt/ct c -alph1*theta(m)*dftxt/xt)+ c mcp(l,m)*beta*delt1*dstdc*(dctxt- c theta(m)*dftxt) sum2 = sum2 + dsum*xp sum3 = sum3 + dsum*xm dsum = -mcp(l,m)*(tau*dutct*dftkt/ct+ c beta*delt1*dstdc) wrk(m,1) = dsum*xpt wrk(m,2) = dsum*xmt 190 continue tem1 = tau*c**(-tau-1.d0)+dsdc m = (l-1)*na work(m+nd1) = sum2 + tem1*xp work(m+nd2) = sum3 + tem1*xm do 232 m = 1,ns mt = (m-1)*na work(mt+nd1t) = work(mt+nd1t) + wrk(m,1) work(mt+nd2t) = work(mt+nd2t) + wrk(m,2) 232 continue res = sum1-c**(-tau)+s tem1 = u(i)*res m = (l-1)*na f(m+nd1) = f(m+nd1) + tem1*xp f(m+nd2) = f(m+nd2) + tem1*xm c tem1 = u(i) do 250 mm=1,nas tem2 = tem1*work(mm) df(m+nd1,mm) = df(m+nd1,mm) + tem2*xp df(m+nd2,mm) = df(m+nd2,mm) + tem2*xm 250 continue 170 continue 150 continue 140 continue c call dgeco(df, nas, nas, iwork, rcond, work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when solving df*x=f' call dgesl(df, nas, nas, iwork, f, 0) sum1 = 0.d0 do 300 j=1,ns do 300 i=1,na tem1 = f((j-1)*na+i) c1(i,j) = c0(i,j) - tem1 sum1 = sum1 + tem1*tem1 300 continue sum1 = dsqrt(sum1)/float(nas) c write(*,*) ' at iteration ', it,' the residual is ', sum1 if (sum1.lt.crit) then go to 999 endif do 310 j=1,ns do 310 i=1,na c0(i,j) = c1(i,j) 310 continue 100 continue 999 continue c c print results c write(7,*) 'RESULTS FOR QR2C.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)) c write(7,*) write(7,*) ' Coefficients for Shape functions (consumption)' do 320 i=1,na write(7,802) (c0(i,j),j=1,ns) 320 continue write(7,*) write(7,*) ' Coefficients for Shape functions (investment)' do 325 i=1,na write(7,802) (theta(j)*xa(i)**alpha-c0(i,j),j=1,ns) 325 continue write(7,*) write(7,*) ' Number of iterations = ', it write(8,806) alpha,beta,delta,tau,gam do 330 i=1,ns write(8,*) theta(i) 330 continue do 340 i=1,ns write(8,*) (mcp(i,j),j=1,ns) 340 continue do 360 i=1,na write(8,802) (c0(i,j),j=1,ns) 360 continue do 370 i=1,nx write(8,804) xa(i),nxp(i) 370 continue write(8,801) xa(nx1) 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