c PROGRAM: SGM2DQS.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} = beta int c~^{-tau}(alpha*exp(tanh^{-1}(y~)) c *x~^(alpha-1)+1-delta) 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 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)*ca, c where ca'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) sgm2dqs.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 initialize option to provide 1 (=yes) or 0 c initial guess for c c nbc number of boundary integer c conditions c x,y,c mesh and initial na x 3 real matrix c (for each guesses for if initialize = 1 c global consumption and na x 2 other- c node) wise, where na = c # of global nodes c elemental global nodes corre- ne x 8 integer c nodes sponding to each matrix, where ne = c element # of elements c (see EXAMPLE below) 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 |43 44 |45 46 |47 48 |49 50 |51 c ^ | | | | | c | *38 *39 *40 *41 *42 c | | | | | | c |e=9 |e=10 |e=11 |e=12 | c y *----*----*----*----*----*----*----*----* c |29 30 |31 32 |33 34 |35 36 |37 c p | | | | | c a *24 *25 *26 *27 *28 c r | | | | | c t |e=5 |e=6 |e=7 |e=8 | c i *----*----*----*----*----*----*----*----* c t |15 16 |17 18 |19 20 |21 22 |23 c i | | | | | c o *10 *11 *12 *13 *14 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 The *'s and e's can be marked anyway the user wants c but the labeling must be consistent with the input c file. For this mesh, the number of nodes (na) is c 51 and the number of elements (ne) is 12. There c is one type of element: an 8-node serendipity c element. When constructing local approximations, c I use the following nodal ordering: c c 4 7 3 c *----*----* c | | c 8* *6 (Local nodes 1-8) c | | c *----*----* c 1 5 2 c c To make the transformation between local and c global nodes, I require that the user type in the c global nodes corresponding to the 8 local nodes. c For example, for elements 1-4, I require that the c following lines be put in the input file (under c ``elemental nodes''): c c 1,3,17,15,2,11,16,10 c 3,5,19,17,4,12,18,11 c 5,7,21,19,6,13,20,12 c 7,9,23,21,8,14,22,13 c c Notes: c i. y=tanh(ln(theta)) -- remember to convert c theta values for input file; c ii. see sgm2dqs.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 v. The Matlab file nodes8se.m can be used to c generate the elemental nodes. c c c OUTPUT FILES: (b) sgm2dqs.nxt: new input file with solution for c(x,y) c (c) sgm2dqs.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: Jeff Eishen (notes) c Taylor and Uhlig (1990), JBES c McGrattan, Federal Reserve Bank of Minneapolis, SR#164 c c SEE ALSO: sgm1dl.f sgm1dq.f, sgm2dq.f, sgm2d.f c c Ellen McGrattan, 2-20-91 c Revised, 4-2-96, ERM c * c c Set dimensions: c na = number of nodes c ne = number of elements 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) i=1:ne. c parameter (na=73,ne=18,nep=10,nq=5) implicit double precision (a-h,o-z) double precision f(na), df(na,na), work(na), xa(na), ya(na), c ax(nq,ne), wx(nq,ne), ay(nq,ne), wy(nq,ne), c omega(nep), we(nep), c0(na), c1(na) integer iwork(na), nxp(ne), nyp(ne), ibc(na), in(ne,8) c c Open input and output files. c open(unit=5, file='sgm2dqs.inp') open(unit=7, file='sgm2dqs.nxt') open(unit=8, file='sgm2dqs.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,*) initialize read(5,*) nbc alph1 = alpha-1.d0 delt1 = 1.d0-delta xmin = 1000000.d0 ymin = 1.d0 xmax = 0.d0 ymax =-1.d0 if (initialize.eq.1) then do 10 i=1,na read(5,*) xa(i),ya(i),c0(i) xmax = dmax1(xmax,xa(i)) xmin = dmin1(xmin,xa(i)) ymax = dmax1(ymax,ya(i)) ymin = dmin1(ymin,ya(i)) 10 continue else c c If initial guess for consumption not given, set c = .5*output. c do 20 i=1,na read(5,*) xa(i),ya(i) c0(i) = 0.5d0*dsqrt((1.d0+ya(i))/(1.d0-ya(i)))*xa(i)**alpha xmax = dmax1(xmax,xa(i)) xmin = dmin1(xmin,xa(i)) ymax = dmax1(ymax,ya(i)) ymin = dmin1(ymin,ya(i)) 20 continue endif do 30 i=1,ne read(5,*) (in(i,j),j=1,8) 30 continue do 40 i=1,nbc read(5,*) ii,tem ibc(ii) = 1 c0(ii) = tem 40 continue c c Other parameters that are set inside the program c * nxp(i), i=1,ne are the # of quadrature points for x on element i c * nyp(i), i=1,ne are the # of quadrature points for y on element i c * stopping criteria c * maximum number of interations c do 70 i=1,ne nxp(i) = 3 70 continue do 80 i=1,ne nyp(i) = 3 80 continue crit = 0.0000001d0 maxit = 1000 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 90 i=1,nep we(i) = we(i)*sqpii sum1 = sum1 + we(i) 90 continue do 100 i=1,nep we(i) = we(i)/sum1 100 continue do 110 i=1,ne call qgausl(nxp(i), xa(in(i,1)), xa(in(i,2)), ax(1,i), wx(1,i)) call qgausl(nyp(i), ya(in(i,1)), ya(in(i,4)), ay(1,i), wy(1,i)) 110 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 130 it=1,maxit c c Initialize f and its derivative df. c do 140 i=1,na f(i) = 0.d0 do 140 j=1,na df(i,j) = 0.d0 140 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 150 n=1,ne c c Assign numbers to the nodes of element n. c nd1 = in(n,1) nd2 = in(n,2) nd3 = in(n,3) nd4 = in(n,4) nd5 = in(n,5) nd6 = in(n,6) nd7 = in(n,7) nd8 = in(n,8) c c Element n is the rectangle with 8 nodes. c x1 = xa(nd1) x2 = xa(nd2) dx = x2-x1 y1 = ya(nd1) y2 = ya(nd4) dy = y2-y1 c c Compute f at all quadrature points on the element. c do 170 j=1,nyp(n) y = ay(j,n) v = wy(j,n) do 170 i=1,nxp(n) x = ax(i,n) u = wx(i,n) 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*(1.d0-xi)*(1.d0-eta)*(-1.d0-xi-eta) bs2 = 0.25*(1.d0+xi)*(1.d0-eta)*(-1.d0+xi-eta) bs3 = 0.25*(1.d0+xi)*(1.d0+eta)*(-1.d0+xi+eta) bs4 = 0.25*(1.d0-xi)*(1.d0+eta)*(-1.d0-xi+eta) bs5 = 0.50*(1.d0-xi*xi)*(1.d0-eta) bs6 = 0.50*(1.d0+xi)*(1.d0-eta*eta) bs7 = 0.50*(1.d0-xi*xi)*(1.d0+eta) bs8 = 0.50*(1.d0-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 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 Compute next period's capital stock. c xt = dexp(tech)*x**alpha+delt1*x-c dftxt = alpha*xt**(alpha-1.d0) 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)) c c Find element containing point (xt,yt). c xtem = dmax1(xt,xmin) xtem = dmin1(xtem,xmax) ytem = dmax1(yt,ymin) ytem = dmin1(ytem,ymax) do 210 l=1,ne i1 = in(l,1) i2 = in(l,2) i4 = in(l,4) if ((xtem.ge.xa(i1)).and.(ytem.ge.ya(i1)).and. c (xtem.le.xa(i2)).and.(ytem.le.ya(i4))) jj = l 210 continue nd1t = in(jj,1) nd2t = in(jj,2) nd3t = in(jj,3) nd4t = in(jj,4) nd5t = in(jj,5) nd6t = in(jj,6) nd7t = in(jj,7) nd8t = in(jj,8) x1t = xa(nd1t) x2t = xa(nd2t) dxt = x2t-x1t y1t = ya(nd1t) y2t = ya(nd4t) dyt = y2t-y1t c c Assign numbers to the basis functions defined on element c of (xt,yt) and take derivatives with respect to xt. c xi = 2.d0*(xt-x1t)/(x2t-x1t)-1.d0 dxi = 2.d0/(x2t-x1t) eta = 2.d0*(yt-y1t)/(y2t-y1t)-1.d0 bs1t = 0.25*(1.d0-xi)*(1.d0-eta)*(-1.d0-xi-eta) bs2t = 0.25*(1.d0+xi)*(1.d0-eta)*(-1.d0+xi-eta) bs3t = 0.25*(1.d0+xi)*(1.d0+eta)*(-1.d0+xi+eta) bs4t = 0.25*(1.d0-xi)*(1.d0+eta)*(-1.d0-xi+eta) bs5t = 0.50*(1.d0-xi*xi)*(1.d0-eta) bs6t = 0.50*(1.d0+xi)*(1.d0-eta*eta) bs7t = 0.50*(1.d0-xi*xi)*(1.d0+eta) bs8t = 0.50*(1.d0-xi)*(1.d0-eta*eta) dbs1 = 0.25*(1.d0-eta)*(2.d0*xi+eta)*dxi dbs2 = 0.25*(1.d0-eta)*(2.d0*xi-eta)*dxi dbs3 = 0.25*(1.d0+eta)*(2.d0*xi+eta)*dxi dbs4 = 0.25*(1.d0+eta)*(2.d0*xi-eta)*dxi dbs5 =-xi*(1.d0-eta)*dxi dbs6 = 0.50*(1.d0-eta*eta)*dxi dbs7 =-xi*(1.d0+eta)*dxi dbs8 =-0.50*(1.d0-eta*eta)*dxi 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 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 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 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 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 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 220 continue 170 continue 150 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 130 continue 999 continue c c Write new input file with results. c write(7,'(5X,F12.6,30X,''alpha in (0,1) '')') alpha write(7,'(5X,F12.6,30X,''beta in (0,1) '')') beta write(7,'(5X,F12.6,30X,''delta in [0,1] '')') delta write(7,'(5X,F12.6,30X,''rho in (-1,1)'')') rho write(7,'(5X,F12.6,30X,''sige >= 0 '')') sige write(7,'(5X,F12.6,30X,''tau >= 0 '')') tau write(7,'(5X,F12.6,30X,''gamma >= 0 '')') gamma write(7,'(5X,I12,30X,''read initial guess for c (1=yes)'')') 1 write(7,'(5X,I12,30X,''number of boundary constraints'')') nbc write(7,*) write(7,'(5X,F12.6,1X,F12.6,1X,F12.6,4X, c''k, z, c(k,z) at nodes 1,2,...A'')') xa(1),ya(1),c1(1) do 280 i=2,na write(7,2) xa(i),ya(i),c1(i) 280 continue write(7,*) write(7,'(6X,8(1X,I3),9X, c''global nodes for each element'')') (in(1,j),j=1,8) do 290 i=2,ne write(7,3) (in(i,j),j=1,8) 290 continue write(7,*) ii = 0 do 300 i=1,na if (ibc(i).eq.1) then if (ii.eq.0) then write(7,'(1X,I3,1X,F12.6,30X,''boundary conditions '')') c i,c1(i) ii = 1 else write(7,4) i,c1(i) endif endif 300 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,*) na do 320 i=1,na write(8,1) xa(i) 320 continue do 330 i=1,na write(8,1) ya(i) 330 continue do 340 i=1,na write(8,1) c1(i) 340 continue c c Format statements. c 1 format (5X,F12.6) 2 format (4X,3(1X,F12.6)) 3 format (6X,8(1X,I3)) 4 format (1X,I3,1X,F12.6) stop end