c PROGRAM: SGM1DQ.F c c DESCRIPTION: implementation of finite-element method for stochastic c growth model, 1-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] drawn from a Markov chain c c[t],k[t] >= 0 c i[t] >= 0 for sufficiently large value of gamma c c PROBLEM: find c=c(x,i) that solves the functional equation: c c c^{-tau} = beta pi(i,j) c~^{-tau}(alpha*theta(j) c *x~^(alpha-1)+1-delta) c c where c c~ = c(x~,j), c x~ = theta(j)*x^alpha+(1-delta)*x-c(x,i) c pi(i,j)= probability of transiting from state i to j c c(0,i) = 0, all i c x in [0,xmax] c i,j in {1,2,...,I} c c ALGORITHM: the finite element method with piecewise bilinear c shape functions N_a(x), ie., c(x,i) = sum N_a(x)*c^i_a, c where c^i_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) sgm1dq.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 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 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 if ievenx=1, c 2 x 1 real vector c with x(1),x(nx) c otherwise, c nx x 1 real vector c with entire grid c ln(theta) grid for technology ns x 1 real vector c pi transition matrix ns*ns x 1 real vectorc (column by column) c c initial guess for nx*ns x 1 real c consumption function vector; read c(:,1) c (if initialize=1) first, c(:,2) c second, etc. c i,j,c points c(x(i),j) triplets of (integer,c that are to be fixed integer,real) c during computation c c EXAMPLE: suppose the mesh is given by c c c | e=1 | e=2 | e=3 | < element #s c *----*-----*-----*---------*----*----* c 1 2 3 4 5 6 7 < node #s c x(1) x(2) x(3) x(4) x(5) x(6) x(7) < x values c c For this example, the length of the x partition (nx) c is 7 and the number of elements (ne) is 3. c c c Notes: c i. see sgm1dq.inp for an example input file; c ii. if grid size changed, edit this file and c change dimensions nx and ny; c iii. all other inputs are set in this file; c see comment lines below. c c OUTPUT FILES: (b) sgm1dq.nxt: new input file with solution for c(x,y) c (c) sgm1dq.dat: input file for Matlab post-processing c routine -- run sgm2dq in Matlab. c c REQUIRED (BLAS) dgemm,dtrsm,idamax,dger,dscal,dswap c PROGRAMS: (LAPACK) dgesv c (PRESS) 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, sgm2dl.f, sgm2dq.f c c Ellen McGrattan, 2-20-91 c Revised, 5-19-97, ERM c * c c Set dimensions: c nx = length of x grid -- must be an odd number c ns = length of theta 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=161,ns=1, nq=5) implicit double precision (a-g,o-z) double precision f(nx*ns), df(nx*ns,nx*ns), work(nx*ns), c xa(nx), ax(nq,(nx+1)/2), wx(nq,(nx+1)/2), theta(ns), c pi(ns,ns), c0(nx,ns), c1(nx,ns) integer iwork(nx*ns), nxp((nx+1)/2), ibc(nx*ns) c na = nx*ns nxm1 = nx-1 ne = (nx+1)/2-1 c c Open input and output files. c open(unit=5, file='sgm1dq.inp') open(unit=7, file='sgm1dq.nxt') open(unit=8, file='sgm1dq.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,*) tau read(5,*) gamma read(5,*) ievenx 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,nxm1 xa(i) = xa(i-1)+tem 10 continue else do 20 i=1,nx read(5,*) xa(i) 20 continue endif do 30 i=1,ns read(5,*) theta(i) theta(i) = dexp(theta(i)) 30 continue do 40 j=1,ns do 40 i=1,ns read(5,*) pi(i,j) 40 continue if (initialize.eq.1) then do 50 j=1,ns do 50 i=1,nx read(5,*) c0(i,j) ibc(i+(j-1)*nx) = 0 50 continue else c c If initial guess for consumption not given, set c = .5*output. c do 60 j=1,ns do 60 i=1,nx ii = i+(j-1)*nx c0(i,j) = 0.5*(theta(j)*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,jj,tem ibc(ii+(jj-1)*nx) = 1 c0(ii,jj) = 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 * stopping criteria c * maximum number of interations c do 80 i=1,ne nxp(i) = 3 80 continue crit = 0.0000001d0 maxit = 20 c c Compute quadrature abscissas and weights. c do 110 i=1,ne call qgausl (nxp(i), xa(2*i-1), xa(2*i+1), ax(1,i), wx(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,j;c0) N_n(x) dx and its derivative on c each element n. c do 150 j=1,ns do 160 n=1,ne c c Element n is the interval [xa(n),xa(n+1)] c x1 = xa(2*n-1) dx = xa(2*n+1)-x1 a = (xa(2*n)-x1)/dx a1 = 1.d0-a c c Compute f at all quadrature points on the element. c do 170 i=1,nxp(n) x = ax(i,n) u = wx(i,n) xi = (x-x1)/dx bs1 = (1.d0-xi)*(1.d0-xi/a) bs2 = xi/(a*a1)*(1.d0-xi) bs3 = -a*xi/a1*(1.d0-xi/a) c c Compute consumption using finite element approximation. c c = c0(2*n-1,j)*bs1+c0(2*n,j)*bs2+c0(2*n+1,j)*bs3 c c Compute the derivatives of the penalty function (s,ds). c y = theta(j)*x**alpha d = dmin1(y-c,0.d0) s = gamma*d*d ds = -2.d0*gamma*d c c Find next period's capital stock on grid. c xt = y+delt1*x-c ii = 1 do 180 k=1,nxm1 if (xt.ge.xa(k)) ii = k 180 continue ii = ii-mod(ii+1,2) x1t = xa(ii) dxt = (xa(ii+2)-x1t) at = (xa(ii+1)-x1t)/dxt a1t = 1.d0-at xi = (xt-x1t)/dxt dxi = 1.d0/dxt bs1t = (1.d0-xi)*(1.d0-xi/at) bs2t = xi/(at*a1t)*(1.d0-xi) bs3t = -at*xi/a1t*(1.d0-xi/at) dbs1 = (-at-1.d0+2.d0*xi)*dxi/at dbs2 = (1.d0-2.d0*xi)*dxi/at/a1t dbs3 = (-at+2.d0*xi)*dxi/a1t dftxt = alpha*xt**alph1 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 sum in the Euler equation: c c beta sum pi(j,k) c~^{-tau}(alpha*theta(k) c *x~^(alpha-1)+1-delta) c do 200 k=1,ns c c Given next period's state (xt,theta(k)), 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), derivatives of the penalty function c next period (st,dst), and the integral of the Euler c equation (sum1). c ct = c0(ii,k)*bs1t+c0(ii+1,k)*bs2t+ c c0(ii+2,k)*bs3t dutct = beta*ct**(-tau) dctxt = c0(ii,k)*dbs1+c0(ii+1,k)*dbs2+ c c0(ii+2,k)*dbs3 dftkt = theta(k)*dftxt+delt1 d = dmin1(theta(k)*xt**alpha-ct,0.d0) st = gamma*d*d dst = -2.d0*gamma*d sum1 = sum1 + pi(j,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 + pi(j,k)*(dutct*(tau*dctxt*dftkt/ c ct-alph1*theta(k)*dftxt/xt)-beta*delt1* c dst*(theta(k)*dftxt-dctxt)) dsum2 = pi(j,k)*(-tau*dutct*dftkt/ct-beta*delt1* c dst) m = ii+(k-1)*nx work(m) = work(m) + dsum2*bs1t work(m+1) = work(m+1) + dsum2*bs2t work(m+2) = work(m+2) + dsum2*bs3t 200 continue c c The vector `work' contains derivatives of the residual. c dres1 = dsum1+tau*c**(-tau-1.d0)+ds m = 2*n-1+(j-1)*nx work(m) = work(m) + dres1*bs1 work(m+1) = work(m+1) + dres1*bs2 work(m+2) = work(m+2) + dres1*bs3 c c Add the residual x quadrature weights x shape functions to f(.). c res = (sum1-c**(-tau)+s)*u f(m) = f(m) + res*bs1 f(m+1) = f(m+1) + res*bs2 f(m+2) = f(m+2) + res*bs3 c c Add derivative of the residual x quadrature weights x shape c functions to df(.). c do 220 k=1,na tem1 = u*work(k) df(m,k) = df(m,k) + tem1*bs1 df(m+1,k) = df(m+1,k) + tem1*bs2 df(m+2,k) = df(m+2,k) + tem1*bs3 220 continue 170 continue 160 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 SGM1DL' c c Do Newton update. c sum1 = 0.d0 ii = 0 do 260 j=1,ns do 260 i=1,nx if (ibc(i+(j-1)*nx).eq.0) then ii = ii+1 c1(i,j) = c0(i,j) - f(ii) sum1 = sum1 + f(ii)*f(ii) else c1(i,j) = c0(i,j) 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 j=1,ns do 270 i=1,nx c0(i,j) = c1(i,j) 270 continue 130 continue 999 continue c c Write new input file with results. c write(7,'(9X,F12.6,12X,''alpha in (0,1) '')') alpha write(7,'(9X,F12.6,12X,''beta in (0,1) '')') beta write(7,'(9X,F12.6,12X,''delta in [0,1] '')') delta write(7,'(9X,F12.6,12X,''tau >= 0 '')') tau write(7,'(3X,F18.6,12X,''gamma >= 0 '')') gamma write(7,'(9X,I12,12X,''evenly-spaced grid? (1=yes)'')') ievenx write(7,'(20X,''1'',12X,''read initial guess for c (1=yes)'')') write(7,'(9X,I12,12X,''number of boundary constraints'')') nbc write(7,*) write(7,'(9X,F12.6,12X,''lower bound for capital grid'')') xa(1) if (ievenx.eq.0) then do 280 i=2,nxm1 write(7,1) xa(i) 280 continue endif write(7,'(9X,F12.6,12X,''upper bound for capital grid'')') xa(nx) write(7,*) write(7,'(9X,F12.6,12X, c ''grid for log of the technology shock'')') dlog(theta(1)) do 281 i=2,ns write(7,1) dlog(theta(i)) 281 continue write(7,*) write(7,'(9X,F12.6,12X, c ''transition probabilities column by column'')') pi(1,1) do 290 i=2,ns write(7,1) pi(i,1) 290 continue do 292 j=2,ns do 292 i=1,ns write(7,1) pi(i,j) 292 continue write(7,*) write(7,'(9X,F12.6,12X,''consumption function '')') c1(1,1) do 300 i=2,nx write(7,1) c1(i,1) 300 continue do 305 j=2,ns do 305 i=1,nx write(7,1) c1(i,j) 305 continue write(7,*) ii = 0 do 310 j=1,ns do 310 i=1,nx if (ibc(i+(j-1)*nx).eq.1) then if (ii.eq.0) then write(7,'(1X,I3,1X,I3,1X,F12.6,12X, c ''boundary conditions '')') i,j,c1(i,j) ii = 1 else write(7,2) i,j,c1(i,j) 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) tau write(8,1) gamma write(8,*) nx do 320 i=1,nx write(8,1) xa(i) 320 continue write(8,*) ns do 330 i=1,ns write(8,1) theta(i) 330 continue do 340 j=1,ns do 340 i=1,ns write(8,1) pi(i,j) 340 continue do 350 j=1,ns do 350 i=1,nx write(8,1) c1(i,j) 350 continue c c Format statements. c 1 format (9X,F12.6) 2 format (1X,I3,1X,I3,1X,F12.6) stop end