* * MAC33NET.F, PROGRAM USED TO SOLVE FOR THE HETEROGENEOUS AGENT * MODEL IN MY PAPER * * "SOLVING DYNAMIC MODELS WITH AGGREGATE SHOCKS AND * HETEROGENEOUS AGENTS" * * Journal of Macroeconomic Dynamics (forthcoming) * * THIS PROGRAM CORRESPONDS WITH EXP2Q30 IN THE PAPER * * * the notation in this program is not consistent with the * notation in the paper. In particular, the program follows * the Heaton and Lucas notation, i.e. each agent receives * a share of a pie. In the paper each agent receives income * that is expressed relative to the per capita income level * * thus the income realization in the notation of the paper * is twice the value of the income realization in the program * implicit real*8(a-h,o-z),integer(i-n) parameter(npar=50,npar2=20) common/state/ stu(1000,10),stbar(10,2),st(1000,10) common/cheb/ cheb(4000,50),cheb2(500,500) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome common/prob/ prob(2,2),qx(24),qw(24),pagg(2,2) common/coef/ aggold(4,100),aggnew(4,100) common/hermite/ hhx(1000),hhw(1000),hxr(30),hxp(30),hw(30) common/solut/ aold(npar),anew(npar) dimension root(3000), $ roothx(1000),roothw(1000) c c npar1 is number of parameters for individual policy function c npar2 is number of parameters for distribution an bond price c open(33,file='bigcoef.b33',status='unknown') open(34,file='bigcoef.new',status='unknown') open(35,file='bigpar.b33',status='unknown') open(36,file='cheb0.dat', status='unknown') open(37,file='pr.dat', status='unknown') open(39,file='herm.dat', status='unknown') open(40,file='herm.pro', status='unknown') open(41,file='bigagg.b33',status='unknown') open(42,file='bigagg.new',status='unknown') read(36,*) (root(j),j=1,870) read(39,*) (roothx(j),j=1,780) read(40,*) (roothw(j),j=1,780) nhermite = 30 do 910 i = 1,nhermite istart = (nhermite-5)*30 hhx(i) = roothx(istart+i) hhw(i) = roothw(istart+i) 910 continue read(37,*) prob(1,1),prob(2,2),prob(1,2),prob(2,1), $ pagg(1,1),pagg(2,2),pagg(1,2),pagg(2,1) c c c In an earlier program I allowed for an additional income c shock that could take on nincome values. In this program c nincome = 1, the value of the additional shock is equal to c zero with probability equal to 1. c nincome = 1 do 5 j = 1,nincome qx(1) = 0. qw(1) = 1. 5 continue c c c prob(i,j) = probability of going from j to i c c do 666 ijk = 1,2 prob(1,1) = 0.5 prob(1,2) = 0.5 prob(2,1) = 0.5 prob(2,2) = 0.5 read(35,*) i1max,i2max,i3max,i4max,i5max,i6max, $ is1max,is2max,is3max,is4max,is5max,is6max, $ blow,bhigh,bondp,beta,gam,eps,eps2, $ bl4,bu4,bl5,bu5,bl6,bu6 read(33,*) (aold(j),j=1,npar) read(41,*) ((aggold(i,j),j=1,npar2),i=1,4) 4 continue imax = ipoint(i1max,i2max,i3max) + 18 if(imax.ne.npar) stop * * * * Construct Gridpoints * * * stbar(.,.) contains the boundaries * * stu(1,2) = 0.3772 stu(2,2) = 0.6228 stu(1,3) = 0.991 stu(2,3) = 1.047 stbar(1,1) = stu(1,2) + qx(1) + blow stbar(1,2) = stu(2,2) + qx(nincome) + bhigh stbar(4,1) = bl4 stbar(4,2) = bu4 stbar(5,1) = bl5 stbar(5,2) = bu5 stbar(6,1) = bl6 stbar(6,2) = bu6 c c c note that implicitly there is a scaling for the discrete c valued state variables (i=2,3). c c Since you never interpolate these values, you don't care c what this scaling is [it is not sca(.)] c c c below we first find the appropriate roots, and then c find the corresponding unscaled numbers c c do 10 is1 = 1,is1max istart = (is1max-2)*30 st(is1,1) = root(istart+is1) 10 continue do 20 is1 = 1,is1max bb = st(is1,1) stu(is1,1) = scinv(bb,1) 20 continue do 11 is2 = 1,is2max istart = (is2max-2)*30 st(is2,2) = root(istart+is2) 11 continue do 12 is3 = 1,is3max istart = (is3max-2)*30 st(is3,3) = root(istart+is3) 12 continue do 13 is4 = 1,is4max istart = (is4max-2)*30 st(is4,4) = root(istart+is4) 13 continue do 21 is4 = 1,is4max bb = st(is4,4) stu(is4,4) = scinv(bb,4) 21 continue do 14 is5 = 1,is5max istart = (is5max-2)*30 st(is5,5) = root(istart+is5) 14 continue do 22 is5 = 1,is5max bb = st(is5,5) stu(is5,5) = scinv(bb,5) 22 continue do 15 is6 = 1,is6max istart = (is6max-2)*30 st(is6,6) = root(istart+is6) 15 continue do 23 is6 = 1,is6max bb = st(is6,6) stu(is6,6) = scinv(bb,6) 23 continue c c c Construct the chebyshev polynomials at the grid points c (Since these have to be calculated many times, it is c worthwhile to save them). c Note that no complete basis is used. c do 40 is6 = 1,is6max do 40 is5 = 1,is5max do 40 is4 = 1,is4max do 40 is3 = 1,is3max do 40 is2 = 1,is2max do 40 is1 = 1,is1max jt = ispoint(is1,is2,is3,is4,is5,is6) ss1 = st(is1,1) ss2 = st(is2,2) ss3 = st(is3,3) ss4 = st(is4,4) ss5 = st(is5,5) ss6 = st(is6,6) do 41 j1 = 1,i1max+1 do 41 j2 = 1,i2max+1 do 41 j3 = 1,i3max+1 i1 = j1 - 1 i2 = j2 - 1 i3 = j3 - 1 j = ipoint(i1,i2,i3) cheb(jt,j) = hh(ss1,ss2,ss3,i1,i2,i3) 41 continue jjmax = ipoint(i1max,i2max,i3max) cheb(jt,jjmax+1) = ss4 cheb(jt,jjmax+2) = ss5 cheb(jt,jjmax+3) = ss6 cheb(jt,jjmax+4) = h(ss4,2) cheb(jt,jjmax+5) = h(ss5,2) cheb(jt,jjmax+6) = h(ss6,2) cheb(jt,jjmax+7) = ss4*ss5 cheb(jt,jjmax+8) = ss4*ss6 cheb(jt,jjmax+9) = ss5*ss6 cheb(jt,jjmax+10) = ss4*ss3 cheb(jt,jjmax+11) = ss5*ss3 cheb(jt,jjmax+12) = ss6*ss3 cheb(jt,jjmax+13) = h(ss4,2)*ss3 cheb(jt,jjmax+14) = h(ss5,2)*ss3 cheb(jt,jjmax+15) = h(ss6,2)*ss3 cheb(jt,jjmax+16) = ss4*ss5*ss3 cheb(jt,jjmax+17) = ss4*ss6*ss3 cheb(jt,jjmax+18) = ss5*ss6*ss3 40 continue do 42 is6 = 1,is6max do 42 is5 = 1,is5max do 42 is4 = 1,is4max do 42 is3 = 1,is3max jt = ispoint2(is3,is4,is5,is6) ss3 = st(is3,3) ss4 = st(is4,4) ss5 = st(is5,5) ss6 = st(is6,6) cheb2(jt,1) = 1. cheb2(jt,2) = ss3 cheb2(jt,3) = ss4 cheb2(jt,4) = ss5 cheb2(jt,5) = ss6 cheb2(jt,6) = ss3*ss4 cheb2(jt,7) = ss3*ss5 cheb2(jt,8) = ss3*ss6 cheb2(jt,9) = h(ss4,2) cheb2(jt,10) = h(ss5,2) cheb2(jt,11) = h(ss6,2) cheb2(jt,12) = ss4*ss5 cheb2(jt,13) = ss4*ss6 cheb2(jt,14) = ss5*ss6 cheb2(jt,15) = h(ss4,2)*ss3 cheb2(jt,16) = h(ss5,2)*ss3 cheb2(jt,17) = h(ss6,2)*ss3 cheb2(jt,18) = ss4*ss5*ss3 cheb2(jt,19) = ss4*ss6*ss3 cheb2(jt,20) = ss5*ss6*ss3 43 continue 42 continue * * * Iterating on Aggregate Policy Rules * * do 448 itera = 1,100 * * * Iterating on S(beta) * do 300 iter = 1,1000 call sbeta() * * checking for accuracy * accur = 0. do 310 i = 1,npar if(abs(aold(i)).gt.0.001) then accur = accur + abs( (aold(i)-anew(i))/aold(i) ) else accur = accur + abs( (aold(i)-anew(i)) ) endif 310 continue c write(*,'(6f13.6)') ( aold(i),i=1,npar) c write(*,'(6f13.6)') ( anew(i),i=1,npar) write(*,*) 'indiv' , accur,iter if(accur.lt.0.05) goto 305 * * use successive approximation to update the coefficients * do 311 i = 1,npar aold(i) = eps*anew(i) + (1.-eps)*aold(i) 311 continue 300 continue 305 continue * * CALCULATE THE AGGREGATE COEFFICIENTS TAKING THE INDIVIDUAL * POLICY RULE AS GIVEN * call update() * * checking for accuracy * accur = 0. do 3310 i = 1,4 do 3310 j = 1,npar2 if(abs(aggold(i,j)).gt.0.001) then accur = accur + abs( (aggold(i,j)-aggnew(i,j))/aggold(i,j) ) else accur = accur + abs( (aggold(i,j)-aggnew(i,j)) ) endif 3310 continue write(*,*) 'macro', accur,ijk if(accur.lt.0.001) goto 3305 * * use successive approximation to update the coefficients * do 3311 i = 1,4 do 3311 j = 1,npar2 aggold(i,j) = eps2*aggnew(i,j) + (1.-eps2)*aggold(i,j) 3311 continue 3300 continue 448 continue 3305 continue write(34,*) (aold(j),j=1,imax) write(42,*) ((aggold(i,j),j=1,npar2),i=1,4) 666 continue stop end C C SUBROUTINE TO CALCULATE THE COEFFICIENTS FOR THE AGGREGATE C POLICY RULES TAKING THE INDIVIDUAL POLICY RULES AS GIVEN C subroutine update() implicit real*8(a-h,o-z),integer(i-n) parameter(npar=50,npar2=20) common/state/ stu(1000,10),stbar(10,2),st(1000,10) common/cheb/ cheb(4000,50),cheb2(500,500) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome common/prob/ prob(2,2),qx(24),qw(24),pagg(2,2) common/coef/ aggold(4,100),aggnew(4,100) common/hermite/ hhx(1000),hhw(1000),hxr(30),hxp(30),hw(30) common/solut/ aold(npar),anew(npar) common/istate/ k3,k4,k5,k6 common/trunc/ tr1,tr2 common/ifu/ ifunc dimension xsol(1),xsol2(2),xx(500,500),xy(500),yreg(4000,4) external funcv,usrf * * k3: current aggregate shock * k4: mean of the wealth distribution of the poor * k5: var of the wealth distribution of the poor * k6: var of the wealth distribution of the rich * do 1000 k3 = 1,is3max do 1000 k4 = 1,is4max do 1000 k5 = 1,is5max do 1000 k6 = 1,is6max jt1 = ispoint2(k3,k4,k5,k6) xg = stu(k3,3) up = stu(k4,4) varp = stu(k5,5) varr = stu(k6,6) ur = - up c c the algorithm has to find a normal which when all the c mass below blow/xg and above bhigh/xg has an average c equal to tr1 and a variance equal to tr2 c This subroutine that calculates the error for a given c mean and variance is "funcv", "Newt" is an equation solver, c and xsol2 is the vector with the initial guesses and the c solution c numpar = 2 tr1 = up tr2 = varp xsol2(1) = up xsol2(2) = varp call newt(xsol2,numpar,check) up = xsol2(1) sigp = sqrt(xsol2(2)) tr1 = ur tr2 = varr xsol2(1) = ur xsol2(2) = varr call newt(xsol2,numpar,check) ur = xsol2(1) sigr = sqrt(xsol2(2)) c c c rescale the hermite roots and weights c c Recall that the hermite stuff is for using a density exp(-x^2) c c if you take the transformation y = sig*sqrt(2)*x + u, c then the density of y is that of the normal (except for c the division by sqrt(pi). (Don't forget to divide by the c Jacobian to understand this, i.e. divide by sig*sqrt(2) c c do 920 i = 1,nhermite hxr(i) = sigr*sqrt(2.)*hhx(i) + ur if(hxr(i).lt.blow/xg) hxr(i) = blow/xg if(hxr(i).gt.bhigh/xg) hxr(i) = bhigh/xg hxp(i) = sigp*sqrt(2.)*hhx(i) + up if(hxp(i).lt.blow/xg) hxp(i) = blow/xg if(hxp(i).gt.bhigh/xg) hxp(i) = bhigh/xg hw(i) = hhw(i)/sqrt(3.141592654) 920 continue * * * calculate the equilibrium bond price * * qlow = 0.95 qhigh = 1. accrtbis = 0.0001 price(1) = rtsec(usrf,qlow,qhigh,accrtbis) yreg(jt1,1) = price(1) * * * calculate the moments * * avepoor = 0. varpoor = 0. varrich = 0. do 1300 k = 1,nincome do 1300 j = 1,nhermite wpoor = hxp(j) + qx(k) + stu(1,2) indiv = 1 bpoor = bfunc(wpoor,indiv,k3,k4,k5,k6) wrich = hxr(j) + qx(k) + stu(2,2) indiv = 2 brich = bfunc(wrich,indiv,k3,k4,k5,k6) w1 = bpoor w2 = brich avepoor = w1*qw(k)*hw(j)*prob(1,1) + $ w2*qw(k)*hw(j)*prob(1,2) + avepoor varpoor = w1**2.*qw(k)*hw(j)*prob(1,1) + $ w2**2.*qw(k)*hw(j)*prob(1,2) + varpoor varrich = w1**2.*qw(k)*hw(j)*prob(2,1) + $ w2**2.*qw(k)*hw(j)*prob(2,2) + varrich 1300 continue z1 = avepoor z2 = varpoor-avepoor**2. z3 = varrich-avepoor**2. c c c note that these formula's implicitly assume that there are c as many poor as rich (this will be true in the steady state c if prob(i,1) + prob(i,2) = 1) c c in general the formula should be something like c c wealth(1,1)*prob(1,1)*(number of poor last period / number of poor this period) c + c wealth(1,2)*prob(1,2)*(number of rich last period / number of poor this period) c yreg(jt1,2) = sca(z1,4) yreg(jt1,3) = sca(z2,5) yreg(jt1,4) = sca(z3,6) 1000 continue * * do the regression (the orthogonality property is imposed) * do 5000 ireg = 1,4 JMAX = is3max*is4max*is5max*is6max NPMAX = NPAR2 DO 110 J = 1,NPMAX XY(J) = 0.0 XX(J,J) = 0.0 110 CONTINUE DO 130 J = 1,NPMAX DO 130 I = 1,JMAX XX(J,J) = XX(J,J) + CHEB2(I,J)*CHEB2(I,J) 130 CONTINUE DO 135 J = 1,NPMAX XX(J,J) = 1./XX(J,J) 135 CONTINUE DO 140 J = 1,NPMAX DO 140 I = 1,JMAX XY(J) = XY(J) + CHEB2(I,J)*YREG(I,IREG) 140 CONTINUE DO 145 J = 1,NPMAX AGGNEW(IREG,J) = XX(J,J)*XY(J) 145 CONTINUE 5000 continue end c c this function calculates the aggregate demand of bonds c for a given bondprice qq. c real*8 function usrf(qq) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=50) common/state/ stu(1000,10),stbar(10,2),st(1000,10) common/cheb/ cheb(4000,50),cheb2(500,500) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome common/prob/ prob(2,2),qx(24),qw(24),pagg(2,2) common/coef/ aggold(4,100),aggnew(4,100) common/hermite/ hhx(1000),hhw(1000),hxr(30),hxp(30),hw(30) common/solut/ aold(npar),anew(npar) common/istate/ k3,k4,k5,k6 price(1) = qq arich = 0. indiv = 2 do 930 i = 1,nhermite do 930 k = 1,nincome wealth= hxr(i)+ qx(k) + stu(indiv,2) arich = arich + bfunc(wealth,indiv,k3,k4,k5,k6)*hw(i)*qw(k) 930 continue apoor = 0. indiv = 1 do 830 i = 1,nhermite do 830 k = 1,nincome wealth= hxp(i) + qx(k) + stu(indiv,2) apoor = apoor + bfunc(wealth,indiv,k3,k4,k5,k6)*hw(i)*qw(k) 830 continue usrf = arich + apoor end c c the algorithm has to find a normal which when all the c mass below blow/xg and above bhigh/xg has an average c equal to tr1 and a variance equal to tr2 c This subroutine calculates the error when the mean is c equal to x(1) and the variance is equal to x(2) c subroutine funcv(n,x,fvec) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=50) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome common/hermite/ hhx(1000),hhw(1000),hxr(30),hxp(30),hw(30) common/trunc/ tr1,tr2 dimension x(n),fvec(n),hermx(30) ave = x(1) sig = sqrt(x(2)) a1 = 0. a2 = 0. do 930 i = 1,nhermite hermx(i) = sig*sqrt(2.)*hhx(i) + ave if(hermx(i).lt.blow/xg) hermx(i) = blow/xg if(hermx(i).gt.bhigh/xg) hermx(i) = bhigh/xg a1 = a1 + hermx(i) * hhw(i)/sqrt(3.141592654) a2 = a2 + hermx(i)**2. * hhw(i)/sqrt(3.141592654) 930 continue a2 = a2 - a1*a1 fvec(1) = a1 - tr1 fvec(2) = a2 - tr2 end c c subroutine that calculates the coefficients of the individual c policy rule taking the coefficients of the aggregate policy c rule as given c subroutine sbeta() implicit real*8(a-h,o-z),integer(i-n) parameter(npar=50) dimension y(20000),xx(npar,npar),xy(4000) common/state/ stu(1000,10),stbar(10,2),st(1000,10) common/cheb/ cheb(4000,50),cheb2(500,500) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome common/prob/ prob(2,2),qx(24),qw(24),pagg(2,2) common/coef/ aggold(4,100),aggnew(4,100) common/solut/ aold(npar),anew(npar) c c c all variables are scaled relative to the aggregate endowment c c the euler equation is written as c c c(t)**gam = beta * E [ c(t+1)**gam * (agg(t+1)/agg(t)**gam ] c c c (agg(t+1)/agg(t)**gam = stu(.,3) c ktel = 0 do 10 is1 = 1,is1max do 10 is2 = 1,is2max do 10 is3 = 1,is3max do 10 is4 = 1,is4max do 10 is5 = 1,is5max do 10 is6 = 1,is6max ktel = ktel + 1 s1old = st(is1,1) s2old = st(is2,2) s3old = st(is3,3) s4old = st(is4,4) s5old = st(is5,5) s6old = st(is6,6) ieq = 1 bondp = aggfunc(s3old,s4old,s5old,s6old,ieq) jt = ispoint(is1,is2,is3,is4,is5,is6) * * calculate next period's bond holdings * exp1 = polg(npar,aold,is1,is2,is3,is4,is5,is6) cons1 = ( beta*exp1 / bondp )**(1./gam) bond1 = ( stu(is1,1) - cons1 ) / bondp if(bond1.lt.blow) then bond1 = blow endif if(bond1.gt.bhigh) then bond1 = bhigh endif * * calculate conditional expectation * expec1 = 0. do 30 jnew = 1,2 do 30 knew = 1,nincome do 30 lnew = 1,2 salary = ( stu(jnew,2) + qx(knew) ) wealth = bond1/stu(lnew,3) + salary s1 = sca(wealth,1) s2 = st(jnew,2) s3 = st(lnew,3) ieq = 2 s4 = aggfunc(s3old,s4old,s5old,s6old,ieq) s4 = scinv(s4,4) s4 = s4/stu(lnew,3) s4 = sca(s4,4) ieq = 3 s5 = aggfunc(s3old,s4old,s5old,s6old,ieq) s5 = scinv(s5,5) s5 = s5/(stu(lnew,3)*stu(lnew,3)) s5 = sca(s5,5) ieq = 4 s6 = aggfunc(s3old,s4old,s5old,s6old,ieq) s6 = scinv(s6,6) s6 = s6/(stu(lnew,3)*stu(lnew,3)) s6 = sca(s6,6) ieq = 1 bondp = aggfunc(s3,s4,s5,s6,ieq) exp2 = pol(npar,aold,s1,s2,s3,s4,s5,s6) cons2 = ( beta*exp2 / bondp )**(1./gam) bond2 = ( wealth - cons2 ) / bondp if(bond2.lt.blow) then bond2 = blow cons2 = ( wealth - bond2*bondp ) endif if(bond2.gt.bhigh) then bond2 = bhigh cons2 = ( wealth - bond2*bondp ) endif expec1 = expec1 $ + qw(knew)*pagg(lnew,is3)*prob(jnew,is2)*cons2**gam $ * stu(lnew,3)**gam 30 continue y(jt) = expec1 10 continue * * do the regression * JMAX = IS1MAX*IS2MAX*IS3MAX*IS4MAX*IS5MAX*IS6MAX DO 110 J = 1,NPAR XY(J) = 0.0 XX(J,J) = 0.0 110 CONTINUE DO 130 J = 1,NPAR DO 130 I = 1,JMAX XX(J,J) = XX(J,J) + CHEB(I,J)*CHEB(I,J) 130 CONTINUE DO 135 J = 1,NPAR XX(J,J) = 1./XX(J,J) 135 CONTINUE DO 140 J = 1,NPAR DO 140 I = 1,JMAX XY(J) = XY(J) + CHEB(I,J)*Y(I) 140 CONTINUE DO 145 J = 1,NPAR ANEW(J) = XX(J,J)*XY(J) 145 CONTINUE c DO 150 I = 1,5 c WRITE(*,*) I,ANEW(I) c150 CONTINUE c write(*,'(6f13.6)') (anew(i),i=1,npar) end c c since the program uses chebyshev nodes all state variables c must be between -1 and 1. This function does the scaling. c real*8 function sca(x,i) implicit real*8(a-h,o-z),integer(i-n) common/state/ stu(1000,10),stbar(10,2),st(1000,10) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome sca = ( 2.*x - stbar(i,2) - stbar(i,1) ) $ / (stbar(i,2)-stbar(i,1)) end c c since the program uses chebyshev nodes all state variables c must be between -1 and 1. This function gives the inverse, i.e. c for a number between -1 and 1, it gives the original value. c real*8 function scinv(x,i) implicit real*8(a-h,o-z),integer(i-n) common/state/ stu(1000,10),stbar(10,2),st(1000,10) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome scinv = ( stbar(i,2) + stbar(i,1) $ + x*(stbar(i,2) - stbar(i,1)) ) / 2. end * * This function defines the univariate Chebyshev polynomial term * (i-th order) * real*8 function h(x,i) implicit real*8(a-h,o-z),integer(i-n) dimension a(200) a(1) = 1. a(2) = x if(i.le.1) then h = a(i+1) goto 99 endif if(i.gt.1) then do 10 j = 2,i a(j+1) = 2.*x*a(j) - a(j-1) 10 continue h = a(i+1) endif 99 continue end * * This function multi(=3)variate chebysev term * real*8 function hh(s1,s2,s3,i1,i2,i3) implicit real*8(a-h,o-z),integer(i-n) hh = h(s1,i1)*h(s2,i2)*h(s3,i3) end * * this function indicates the location in the vector of * coefficients for the indicated order of the three * state variables. (this is only used for the individual problem) * integer function ipoint(i1,i2,i3) implicit real*8(a-h,o-z),integer(i-n) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome ipoint = 1 + i2 $ + i3*(i2max+1) $ + i1*(i3max+1)*(i2max+1) c c i* goes from 0 to i*max c c with this convention increasing i1max adds parameters at the c end of the parameter vector c end * * each element of the grid is given a number. this function * indicates that number. This is for the grid used to solve * the individual problem. * integer function ispoint(is1,is2,is3,is4,is5,is6) implicit real*8(a-h,o-z),integer(i-n) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome ispoint = is1 $ + (is2-1)*(is1max) $ + (is3-1)*(is1max*is2max) $ + (is4-1)*(is1max*is2max*is3max) $ + (is5-1)*(is1max*is2max*is3max*is4max) $ + (is6-1)*(is1max*is2max*is3max*is4max*is5max) c c is* goes from 1 to is*max c end * * each element of the grid is given a number. this function * indicates that number. This is for the grid used to solve * the aggregate problem. * integer function ispoint2(is3,is4,is5,is6) implicit real*8(a-h,o-z),integer(i-n) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome ispoint2 = is3 $ + (is4-1)*(is3max) $ + (is5-1)*(is3max*is4max) $ + (is6-1)*(is3max*is4max*is5max) c c is* goes from 1 to is*max c end * * this function calculates the value of the aggregate * policy function when the inputs for the four state * variables are equal to ss3, ss4, ss5, and ss6. * Recall that these are the scaled inputs * * * ss3: current aggregate shock * ss4: mean of the wealth distribution of the poor * ss5: variance of the wealth distribution of the poor * ss6: variance of the wealth distribution of the rich * \ * * ireg = 1: bond price * ireg = 2: next period's mean of the poor * ireg = 3: next period's variance of the poor * ireg = 4: next period's variance of the rich * real*8 function aggfunc(ss3,ss4,ss5,ss6,ireg) Implicit real*8(a-h,o-z),integer(i-n) parameter(npar2=20) common/state/ stu(1000,10),stbar(10,2),st(1000,10) common/cheb/ cheb(4000,50),cheb2(500,500) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome common/prob/ prob(2,2),qx(24),qw(24),pagg(2,2) common/coef/ aggold(4,100),aggnew(4,100) common/hermite/ hhx(1000),hhw(1000),hxr(30),hxp(30),hw(30) aggfunc = 0. do 44 i = 1,npar2 if(i.eq.1) term = 1 if(i.eq.2) term = ss3 if(i.eq.3) term = ss4 if(i.eq.4) term = ss5 if(i.eq.5) term = ss6 if(i.eq.6) term = ss4*ss3 if(i.eq.7) term = ss5*ss3 if(i.eq.8) term = ss6*ss3 if(i.eq.9) term = h(ss4,2) if(i.eq.10) term = h(ss5,2) if(i.eq.11) term = h(ss6,2) if(i.eq.12) term = ss4*ss5 if(i.eq.13) term = ss4*ss6 if(i.eq.14) term = ss5*ss6 if(i.eq.15) term = h(ss4,2)*ss3 if(i.eq.16) term = h(ss5,2)*ss3 if(i.eq.17) term = h(ss6,2)*ss3 if(i.eq.18) term = ss4*ss5*ss3 if(i.eq.19) term = ss4*ss6*ss3 if(i.eq.20) term = ss5*ss6*ss3 aggfunc = aggfunc + aggold(ireg,i)*term 44 continue end * * this function calculates the value of the individual * policy function (i.e., the parameterized expectation) * when the inputs for the four state * variables are equal to s1, s2, s3, s4, s5, and s6. * Recall that these are the scaled inputs * * * s1: wealth * s2: income realization * s3: current aggregate shock * s4: mean of the wealth distribution of the poor * s5: variance of the wealth distribution of the poor * s6: variance of the wealth distribution of the rich * real*8 function pol(npar,coef,s1,s2,s3,s4,s5,s6) implicit real*8(a-h,o-z),integer(i-n) common/state/ stu(1000,10),stbar(10,2),st(1000,10) common/cheb/ cheb(4000,50),cheb2(500,500) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome dimension coef(npar) pol = 0. do 20 j1 = 1,i1max+1 do 20 j2 = 1,i2max+1 do 20 j3 = 1,i3max+1 i1 = j1 - 1 i2 = j2 - 1 i3 = j3 - 1 jj = ipoint(i1,i2,i3) pol = pol + coef(jj)*hh(s1,s2,s3,i1,i2,i3) 20 continue jjmax = ipoint(i1max,i2max,i3max) pol = pol + coef(jjmax+1)*s4 pol = pol + coef(jjmax+2)*s5 pol = pol + coef(jjmax+3)*s6 pol = pol + coef(jjmax+4 )*h(s4,2) pol = pol + coef(jjmax+5 )*h(s5,2) pol = pol + coef(jjmax+6 )*h(s6,2) pol = pol + coef(jjmax+7 )*s4*s5 pol = pol + coef(jjmax+8 )*s4*s6 pol = pol + coef(jjmax+9 )*s5*s6 pol = pol + coef(jjmax+10)*s4*s3 pol = pol + coef(jjmax+11)*s5*s3 pol = pol + coef(jjmax+12)*s6*s3 pol = pol + coef(jjmax+13)*h(s4,2)*s3 pol = pol + coef(jjmax+14)*h(s5,2)*s3 pol = pol + coef(jjmax+15)*h(s6,2)*s3 pol = pol + coef(jjmax+16)*s4*s5*s3 pol = pol + coef(jjmax+17)*s4*s6*s3 pol = pol + coef(jjmax+18)*s5*s6*s3 end * * * this function is the same as pol(.) but does it only * at the grid points. Since all the chebyshev terms * at the grid point values are saved, using this function * is faster than using pol(.) * real*8 function polg(npar,coef,is1,is2,is3,is4,is5,is6) implicit real*8(a-h,o-z),integer(i-n) common/state/ stu(1000,10),stbar(10,2),st(1000,10) common/cheb/ cheb(4000,50),cheb2(500,500) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome dimension coef(npar) jt = ispoint(is1,is2,is3,is4,is5,is6) polg = 0. do 20 jj = 1,npar polg = polg + coef(jj)*cheb(jt,jj) 20 continue end * * * this function calculates the demand for bonds at the aggregate * grid point. This function is used to calculate, among other * things, the average demand for bonds. Therefore the wealth * values used vary but all the other inputs are equal to ghe * values at the grid points. * * real*8 function bfunc(wealth,is2,is3,is4,is5,is6) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=50) common/state/ stu(1000,10),stbar(10,2),st(1000,10) common/cheb/ cheb(4000,50),cheb2(500,500) common/rpar/ xg,blow,bhigh,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nhermite, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome common/prob/ prob(2,2),qx(24),qw(24),pagg(2,2) common/coef/ aggold(4,100),aggnew(4,100) common/hermite/ hhx(1000),hhw(1000),hxr(30),hxp(30),hw(30) common/solut/ aold(npar),anew(npar) s1 = sca(wealth,1) s2 = st(is2,2) s3 = st(is3,3) s4 = st(is4,4) s5 = st(is5,5) s6 = st(is6,6) bondp = price(1) exp2 = pol(npar,aold,s1,s2,s3,s4,s5,s6) cons2 = ( beta*exp2 / bondp )**(1./gam) bfunc = ( wealth - cons2 ) / bondp if(bfunc.lt.blow) then bfunc = blow endif if(bfunc.gt.bhigh) then bfunc = bhigh endif c write(*,*) wealth,bfunc,is2,is3 end c c BELOW ARE TWO SUBROUTINES TO NUMERICALLY SOLVE EQUATIONS c NEWT and RTSEC (Both are from numerical recipes). c SUBROUTINE newt(x,n,check) INTEGER n,nn,NP,MAXITS LOGICAL check DOUBLE PRECISION x(n),fvec,TOLF,TOLMIN,TOLX,STPMX PARAMETER (NP=40,MAXITS=200,TOLF=1.d-8,TOLMIN=1.d-12,TOLX=3.d-16, *STPMX=10000.d0) COMMON /newtv/ fvec(NP),nn SAVE /newtv/ CU USES fdjac,fmin,lnsrch,lubksb,ludcmp INTEGER i,its,j,indx(NP) DOUBLE PRECISION d,den,f,fold,stpmax,sum,temp,test,fjac(NP,NP),g *(NP),p(NP), *xold(NP),fmin EXTERNAL fmin nn=n f=fmin(x) test=0.d0 do 11 i=1,n if(abs(fvec(i)).gt.test)test=abs(fvec(i)) 11 continue if(test.lt..01d0*TOLF)return sum=0.d0 do 12 i=1,n sum=sum+x(i)**2 12 continue stpmax=STPMX*max(sqrt(sum), dble(n)) do 21 its=1,MAXITS call fdjac(n,x,fvec,NP,fjac) do 14 i=1,n sum=0.d0 do 13 j=1,n sum=sum+fjac(j,i)*fvec(j) 13 continue g(i)=sum 14 continue do 15 i=1,n xold(i)=x(i) 15 continue fold=f do 16 i=1,n p(i)=-fvec(i) 16 continue call ludcmp(fjac,n,NP,indx,d) call lubksb(fjac,n,NP,indx,p) call lnsrch(n,xold,fold,g,p,x,f,stpmax,check,fmin) test=0.d0 do 17 i=1,n if(abs(fvec(i)).gt.test)test=abs(fvec(i)) 17 continue if(test.lt.TOLF)then check=.false. return endif if(check)then test=0.d0 den=max(f,.5d0*n) do 18 i=1,n temp=abs(g(i))*max(abs(x(i)),1.d0)/den if(temp.gt.test)test=temp 18 continue if(test.lt.TOLMIN)then check=.true. else check=.false. endif return endif test=0.d0 do 19 i=1,n temp=(abs(x(i)-xold(i)))/max(abs(x(i)),1.d0) if(temp.gt.test)test=temp 19 continue if(test.lt.TOLX)return 21 continue pause 'MAXITS exceeded in newt' END C (C) Copr. 1986-92 Numerical Recipes Software D.Q=. SUBROUTINE lubksb(a,n,np,indx,b) INTEGER n,np,indx(n) DOUBLE PRECISION a(np,np),b(n) INTEGER i,ii,j,ll DOUBLE PRECISION sum ii=0 do 12 i=1,n ll=indx(i) sum=b(ll) b(ll)=b(i) if (ii.ne.0)then do 11 j=ii,i-1 sum=sum-a(i,j)*b(j) 11 continue else if (sum.ne.0.d0) then ii=i endif b(i)=sum 12 continue do 14 i=n,1,-1 sum=b(i) do 13 j=i+1,n sum=sum-a(i,j)*b(j) 13 continue b(i)=sum/a(i,i) 14 continue return END C (C) Copr. 1986-92 Numerical Recipes Software D.Q=. SUBROUTINE ludcmp(a,n,np,indx,d) INTEGER n,np,indx(n),NMAX DOUBLE PRECISION d,a(np,np),TINY PARAMETER (NMAX=500,TINY=1.0d-20) INTEGER i,imax,j,k DOUBLE PRECISION aamax,dum,sum,vv(NMAX) d=1.d0 do 12 i=1,n aamax=0.d0 do 11 j=1,n if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) 11 continue if (aamax.eq.0.d0) pause 'singular matrix in ludcmp' vv(i)=1.d0/aamax 12 continue do 19 j=1,n do 14 i=1,j-1 sum=a(i,j) do 13 k=1,i-1 sum=sum-a(i,k)*a(k,j) 13 continue a(i,j)=sum 14 continue aamax=0.d0 do 16 i=j,n sum=a(i,j) do 15 k=1,j-1 sum=sum-a(i,k)*a(k,j) 15 continue a(i,j)=sum dum=vv(i)*abs(sum) if (dum.ge.aamax) then imax=i aamax=dum endif 16 continue if (j.ne.imax)then do 17 k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum 17 continue d=-d vv(imax)=vv(j) endif indx(j)=imax if(a(j,j).eq.0.d0)a(j,j)=TINY if(j.ne.n)then dum=1.d0/a(j,j) do 18 i=j+1,n a(i,j)=a(i,j)*dum 18 continue endif 19 continue return END C (C) Copr. 1986-92 Numerical Recipes Software D.Q=. SUBROUTINE fdjac(n,x,fvec,np,df) INTEGER n,np,NMAX DOUBLE PRECISION df(np,np),fvec(n),x(n),EPS PARAMETER (NMAX=40,EPS=1.d-8) CU USES funcv INTEGER i,j DOUBLE PRECISION h,temp,f(NMAX) do 12 j=1,n temp=x(j) h=EPS*abs(temp) if(h.eq.0.d0)h=EPS x(j)=temp+h h=x(j)-temp call funcv(n,x,f) x(j)=temp do 11 i=1,n df(i,j)=(f(i)-fvec(i))/h 11 continue 12 continue return END C (C) Copr. 1986-92 Numerical Recipes Software D.Q=. FUNCTION fmin(x) INTEGER n,NP DOUBLE PRECISION fmin,x(*),fvec PARAMETER (NP=40) COMMON /newtv/ fvec(NP),n SAVE /newtv/ CU USES funcv INTEGER i DOUBLE PRECISION sum call funcv(n,x,fvec) sum=0.d0 do 11 i=1,n sum=sum+fvec(i)**2 11 continue fmin=0.5d0*sum return END C (C) Copr. 1986-92 Numerical Recipes Software D.Q=. SUBROUTINE lnsrch(n,xold,fold,g,p,x,f,stpmax,check,func) INTEGER n LOGICAL check DOUBLE PRECISION f,fold,stpmax,g(n),p(n),x(n),xold(n),func,ALF *,TOLX PARAMETER (ALF=1.d-4,TOLX=3.d-16) EXTERNAL func CU USES func INTEGER i DOUBLE PRECISION a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2 *,slope,sum,temp, *test,tmplam check=.false. sum=0.d0 do 11 i=1,n sum=sum+p(i)*p(i) 11 continue sum=sqrt(sum) if(sum.gt.stpmax)then do 12 i=1,n p(i)=p(i)*stpmax/sum 12 continue endif slope=0.d0 do 13 i=1,n slope=slope+g(i)*p(i) 13 continue test=0.d0 do 14 i=1,n temp=abs(p(i))/max(abs(xold(i)),1.d0) if(temp.gt.test)test=temp 14 continue alamin=TOLX/test alam=1.d0 1 continue do 15 i=1,n x(i)=xold(i)+alam*p(i) 15 continue f=func(x) if(alam.lt.alamin)then do 16 i=1,n x(i)=xold(i) 16 continue check=.true. return else if(f.le.fold+ALF*alam*slope)then return else if(alam.eq.1.d0)then tmplam=-slope/(2.d0*(f-fold-slope)) else rhs1=f-fold-alam*slope rhs2=f2-fold2-alam2*slope a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2) b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/(alam-alam2) if(a.eq.0.d0)then tmplam=-slope/(2.d0*b) else disc=b*b-3.d0*a*slope tmplam=(-b+sqrt(disc))/(3.d0*a) endif if(tmplam.gt..5d0*alam)tmplam=.5d0*alam endif endif alam2=alam f2=f fold2=fold alam=max(tmplam,.1d0*alam) goto 1 END C (C) Copr. 1986-92 Numerical Recipes Software D.Q=. FUNCTION RTSEC(FUNC,X1,X2,XACC) IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N) PARAMETER (MAXIT=30) FL=FUNC(X1) F=FUNC(X2) IF(ABS(FL).LT.ABS(F))THEN RTSEC=X1 XL=X2 SWAP=FL FL=F F=SWAP ELSE XL=X1 RTSEC=X2 ENDIF DO 11 J=1,MAXIT DX=(XL-RTSEC)*F/(F-FL) XL=RTSEC FL=F RTSEC=RTSEC+DX F=FUNC(RTSEC) IF(ABS(DX).LT.XACC.OR.F.EQ.0.)RETURN 11 CONTINUE PAUSE 'RTSEC exceed maximum iterations' END