c PROGRAM: SGM2DQ.F c c DESCRIPTION: implementation of finite-element method for stochastic c growth model, 2-D, quadratic 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} -s = beta int { c~^{-tau}[alpha c *exp(tanh^{-1}(y~))*x~^(alpha-1) c +1-delta- s~*(1-delta)] c *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 s = gamma* min(exp(tanh^{-1}(y))*x^alpha-c,0)^2 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 quadratic c shape functions N_a(x,y), ie., c(x,y) = sum N_a(x,y)*c_a, c where c_a'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) sgm2dq.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 conditions c x grid for capital nx x 1 real vector c with x at all nodes c on x partition 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)) with y at all nodes c on y-partition 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 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 *----*----*----*----*----*----*----*----* c |55 56 |57 58 |59 60 |61 62 |63 c ^ | | | | | c | *46 *47 *48 *49 *50 *51 *52 *53 *54 c | | | | | | c |e=9 |e=10 |e=11 |e=12 | c y *----*----*----*----*----*----*----*----* c |37 38 |39 40 |41 42 |43 44 |45 c p | | | | | c a *28 *29 *30 *31 *32 *33 *34 *35 *36 c r | | | | | c t |e=5 |e=6 |e=7 |e=8 | c i *----*----*----*----*----*----*----*----* c t |19 20 |21 22 |23 24 |25 26 |27 c i | | | | | c o *10 *11 *12 *13 *14 *15 *16 *17 *18 c n | | | | | c |e=1 |e=2 |e=3 |e=4 | c *----*----*----*----*----*----*----*----* c 1 2 3 4 5 6 7 8 9 c c x partition ----> c c where * marks a node and e's are element numbers. c For this mesh, the number of elements (ne) is 12, c the number of x values (nx) is 9, the number of c y values (ny) is 7, and the number of nodes (na) c is 63. There is one type of element: a 9-node c element. The x-grid in this case is the x values c at nodes 1-9. The y-grid is the y values at c nodes 1,10,19,28,37,46, and 55. If the consumption c function is provided by the user, then the ith c value read in should be consumption at node i. c c c Notes: c i. y=tanh(ln(theta)) -- remember to convert c theta values for input file; c ii. see sgm2dq.inp for an example input file; c iii. if grid size changed, 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) sgm2dq.nxt: new input file with solution for c(x,y) c (c) sgm2dq.dat: input file for Matlab post-processing c routine -- run sgm2d 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, sgm2dqs.f, sgm2d.f c c Ellen McGrattan, 2-20-91 c Revised, 4-2-96, ERM c * c c Set dimensions: c nx = length of x grid -- must be an odd number c ny = length of y grid -- must be an odd number 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=21,ny=11,nep=10,nq=5) implicit double precision (a-g,o-z) double precision f(nx*ny), df(nx*ny,nx*ny), work(nx*ny), c xa(nx), ya(ny), ax(nq,(nx+1)/2), wx(nq,(nx+1)/2), c ay(nq,(ny+1)/2), wy(nq,(ny+1)/2), omega(nep), c we(nep), c0(nx*ny), c1(nx*ny) integer iwork(nx*ny), nxp((nx+1)/2), nyp((ny+1)/2), ibc(nx*ny) c na = nx*ny nex = (nx+1)/2-1 ney = (ny+1)/2-1 ne = nex*ney c c Open input and output files. c open(unit=5, file='sgm2dq.inp') open(unit=7, file='sgm2dq.nxt') open(unit=8, file='sgm2dq.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,nx-1 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(ny-1) do 30 i=2,ny-1 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 tech = 0.5d0*(dlog(1.d0+ya(j))-dlog(1.d0-ya(j))) ii = i+(j-1)*nx c0(ii) = 0.5*(dexp(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,nex nxp(i) = 3 80 continue do 90 i=1,ney nyp(i) = 3 90 continue crit = 0.0000001d0 maxit = 20 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,nex call qgausl (nxp(i), xa(2*i-1), xa(2*i+1), ax(1,i), wx(1,i)) 120 continue do 130 i=1,ney call qgausl (nyp(i), ya(2*i-1), ya(2*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+2)] x [ya(iy),ya(iy+2)]. c iex = mod(n-1,nex)+1 ix = 2*iex-1 iey = (n-1)/nex+1 iy = 2*iey-1 x1 = xa(ix) x2 = xa(ix+2) y1 = ya(iy) y2 = ya(iy+2) c c Assign numbers to the nodes of element n. c nd1 = ix+nx*(iy-1) nd5 = nd1+1 nd2 = nd1+2 nd4 = nd1+2*nx nd7 = nd4+1 nd3 = nd4+2 nd8 = nd1+nx nd9 = nd8+1 nd6 = nd8+2 c c Compute f at all quadrature points on the element. c do 170 j=1,nyp(iey) y = ay(j,iey) v = wy(j,iey) do 170 i=1,nxp(iex) x = ax(i,iex) u = wx(i,iex) 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.25*xi*(xi-1.d0)*eta*(eta-1.d0) bs2 = 0.25*xi*(xi+1.d0)*eta*(eta-1.d0) bs3 = 0.25*xi*(xi+1.d0)*eta*(eta+1.d0) bs4 = 0.25*xi*(xi-1.d0)*eta*(eta+1.d0) bs5 = 0.50*(1.d0-xi*xi)*eta*(eta-1.d0) bs6 = 0.50*xi*(xi+1.d0)*(1.d0-eta*eta) bs7 = 0.50*(1.d0-xi*xi)*eta*(eta+1.d0) bs8 = 0.50*xi*(xi-1.d0)*(1.d0-eta*eta) bs9 = (1.d0-xi*xi)*(1.d0-eta*eta) c c Compute consumption using finite element approximation. c c = c0(nd1)*bs1+c0(nd2)*bs2+c0(nd3)*bs3+ c c0(nd4)*bs4+c0(nd5)*bs5+c0(nd6)*bs6+ c c0(nd7)*bs7+c0(nd8)*bs8+c0(nd9)*bs9 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,nx-1 if (xt.ge.xa(k)) ii = k 180 continue ii = ii-mod(ii+1,2) x1t = xa(ii) x2t = xa(ii+2) 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(-omega^2)/sqrt(pi) d omega. 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,ny-1 if (yt.ge.ya(l)) jj = l 210 continue jj = jj-mod(jj+1,2) y1t = ya(jj) y2t = ya(jj+2) eta = 2.d0*(yt-y1t)/(y2t-y1t)-1.d0 c c Assign numbers to the nodes of element of (xt,yt). c nd1t = ii+nx*(jj-1) nd5t = nd1t+1 nd2t = nd1t+2 nd4t = nd1t+2*nx nd7t = nd4t+1 nd3t = nd4t+2 nd8t = nd1t+nx nd9t = nd8t+1 nd6t = nd8t+2 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.25*xi*(xi-1.d0)*eta*(eta-1.d0) bs2t = 0.25*xi*(xi+1.d0)*eta*(eta-1.d0) bs3t = 0.25*xi*(xi+1.d0)*eta*(eta+1.d0) bs4t = 0.25*xi*(xi-1.d0)*eta*(eta+1.d0) bs5t = 0.50*(1.d0-xi*xi)*eta*(eta-1.d0) bs6t = 0.50*xi*(xi+1.d0)*(1.d0-eta*eta) bs7t = 0.50*(1.d0-xi*xi)*eta*(eta+1.d0) bs8t = 0.50*xi*(xi-1.d0)*(1.d0-eta*eta) bs9t = (1.d0-xi*xi)*(1.d0-eta*eta) dbs1 = bs1t*dxi*(1.d0/xi+1.d0/(xi-1.d0)) dbs2 = bs2t*dxi*(1.d0/xi+1.d0/(xi+1.d0)) dbs3 = bs3t*dxi*(1.d0/xi+1.d0/(xi+1.d0)) dbs4 = bs4t*dxi*(1.d0/xi+1.d0/(xi-1.d0)) dbs5 = bs5t*dxi*2.d0*xi/(xi*xi-1.d0) dbs6 = bs6t*dxi*(1.d0/xi+1.d0/(xi+1.d0)) dbs7 = bs7t*dxi*2.d0*xi/(xi*xi-1.d0) dbs8 = bs8t*dxi*(1.d0/xi+1.d0/(xi-1.d0)) dbs9 = bs9t*dxi*2.d0*xi/(xi*xi-1.d0) 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+c0(nd3t)*bs3t+ c c0(nd4t)*bs4t+c0(nd5t)*bs5t+c0(nd6t)*bs6t+ c c0(nd7t)*bs7t+c0(nd8t)*bs8t+c0(nd9t)*bs9t dutct = beta*ct**(-tau) dctxt = c0(nd1t)*dbs1+c0(nd2t)*dbs2+c0(nd3t)*dbs3+ c c0(nd4t)*dbs4+c0(nd5t)*dbs5+c0(nd6t)*dbs6+ c c0(nd7t)*dbs7+c0(nd8t)*dbs8+c0(nd9t)*dbs9 techt = 0.5d0*(dlog(1.d0+yt)-dlog(1.d0-yt)) dftkt = dexp(techt)*dftxt+delt1 d = dmin1(dexp(techt)*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*dexp(techt)*dftxt/xt)-beta*delt1* c dst*(dexp(techt)*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 work(nd5t) = work(nd5t) + dsum2*bs5t work(nd6t) = work(nd6t) + dsum2*bs6t work(nd7t) = work(nd7t) + dsum2*bs7t work(nd8t) = work(nd8t) + dsum2*bs8t work(nd9t) = work(nd9t) + dsum2*bs9t 200 continue c c The vector `work' contains derivatives of the residual. c dres1 = dsum1+tau*c**(-tau-1.d0)+ds work(nd1) = work(nd1) + dres1*bs1 work(nd2) = work(nd2) + dres1*bs2 work(nd3) = work(nd3) + dres1*bs3 work(nd4) = work(nd4) + dres1*bs4 work(nd5) = work(nd5) + dres1*bs5 work(nd6) = work(nd6) + dres1*bs6 work(nd7) = work(nd7) + dres1*bs7 work(nd8) = work(nd8) + dres1*bs8 work(nd9) = work(nd9) + dres1*bs9 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 f(nd5) = f(nd5) + res*bs5 f(nd6) = f(nd6) + res*bs6 f(nd7) = f(nd7) + res*bs7 f(nd8) = f(nd8) + res*bs8 f(nd9) = f(nd9) + res*bs9 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 df(nd5,k) = df(nd5,k) + tem*bs5 df(nd6,k) = df(nd6,k) + tem*bs6 df(nd7,k) = df(nd7,k) + tem*bs7 df(nd8,k) = df(nd8,k) + tem*bs8 df(nd9,k) = df(nd9,k) + tem*bs9 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 SGM2DQ' 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) 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 in (0,1) '')') alpha write(7,'(5X,F12.6,12X,''beta in (0,1) '')') beta write(7,'(5X,F12.6,12X,''delta in [0,1] '')') delta write(7,'(5X,F12.6,12X,''rho in (-1,1)'')') rho write(7,'(5X,F12.6,12X,''sige >= 0 '')') sige write(7,'(5X,F12.6,12X,''tau >= 0 '')') tau write(7,'(5X,F12.6,12X,''gamma >= 0 '')') 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,''lower bound for tanh(ln(theta))'')') c ya(1) write(7,'(5X,F12.6,12X,''upper bound for tanh(ln(theta))'')') c 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) ii = 1 else write(7,2) i,c1(i) endif 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,1) gamma 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 write(8,*) it c c Format statements. c 1 format (5X,F12.6) 2 format (1X,I3,1X,F12.6) stop end