* * MAC8ADNE(T).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 EXP4Q101 IN THE PAPER * * the user is recommended to first take a look at mac33net.f * which is easier and has more explanatory comments * * 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,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ 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/ hx(1000),hw(1000) common/solut/ aold(npar),anew(npar) common/ddddd/ ein(5,5,5,5,5) common/trunc/ tr1,tr2,tr3,tr4 common/is/ ispoint2(5,5,5,5),jt2max dimension root(3000),xsol2(5) c c npar1 is number of parameters for individual policy function c npar2 is number of parameters for distribution & bond price c open(33,file='bigcoef.bi8',status='unknown') open(34,file='bigcoef.new',status='unknown') open(35,file='bigpar.bi8',status='unknown') open(36,file='cheb0.dat', status='unknown') open(37,file='prob.dat', status='unknown') open(39,file='herm.dat', status='unknown') open(40,file='herm.pro', status='unknown') open(41,file='bigagg.bi8',status='unknown') open(42,file='bigagg.new',status='unknown') open(99,file='initial.dat', status='unknown') read(36,*) (root(j),j=1,870) 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) nincome = 1 do 5 j = 1,nincome if(nincome.eq.1) then qx(1) = 0. qw(1) = 1. endif 5 continue do 666 ijk = 1,3 read(35,*) i1max,i2max,i3max,i4max,i5max,i6max, $ is1max,is2max,is3max,is4max,is5max,is6max, $ blow,bhigh,bondp,beta,gam,eps1,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) if(ijk.eq.3) then pagg(1,1) = 0.7 pagg(1,2) = 0.3 pagg(2,2) = 0.7 pagg(2,1) = 0.3 endif c c c prob(i,j) = probability of going from j to i c c nsimpson = 101 c c nsimpson should be odd c hw(1) = 1./3. hw(nsimpson) = 1./3. hx(1) = blow*2. hx(nsimpson) = bhigh*1.1 bdiff = hx(nsimpson)-hx(1) is = 1 do 910 i = 2,nsimpson-1 if(is.gt.0) then hw(i) = 4./3. is = -1 else hw(i) = 2./3. is = 1 endif write(*,*) hw(i) hx(i) = hx(1) + dble(i-1)*bdiff/dble(nsimpson-1) 910 continue c c note that the simpson points should go outside the c lower and upper bound for the bond holdings to ensure c that I get enough mass at the end points c imax = ipoint(i1max,i2max,i3max) + 18 write(*,*) imax if(imax.ne.npar) stop * * * * Construct Gridpoints * * * stbar(.,.) contains the boundaries * * stu(1,2) = 0.3772 stu(2,2) = 0.6228 if(ijk.eq.1) then stu(1,3) = 0.991 stu(2,3) = 1.047 endif if(ijk.eq.3) then stu(1,3) = 0.991 stu(2,3) = 1.047 endif if(ijk.eq.2) then stu(1,3) = 1./1.02 stu(2,3) = 1.02 endif 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 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 below we first find the appropriate roots, and then c find the corresponding unscaled numbers 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 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 jt = 0 do 42 is3 = 1,is3max do 42 is4 = 1,is4max do 42 is5 = 1,is5max do 42 is6 = 1,is6max if(is4.eq.1.and.is5.ge.is5max-1) goto 42 if(is4.eq.1.and.is6.ge.is6max-1) goto 42 if(is5.eq.1.and.is4.ge.is4max-1) goto 42 if(is5.eq.1.and.is6.ge.is6max-1) goto 42 if(is6.eq.1.and.is4.ge.is4max-1) goto 42 if(is6.eq.1.and.is5.ge.is5max-1) goto 42 if(is4.eq.2.and.is5.ge.is5max) goto 42 if(is4.eq.2.and.is6.ge.is6max) goto 42 if(is5.eq.2.and.is4.ge.is4max) goto 42 if(is5.eq.2.and.is6.ge.is6max) goto 42 if(is6.eq.2.and.is4.ge.is4max) goto 42 if(is6.eq.2.and.is5.ge.is5max) goto 42 jt = jt + 1 write(*,*) jt ispoint2(is3,is4,is5,is6) = jt 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 jt2max = jt * * * calculate the coefficients of the density at the grid points * * do 1000 k3 = 1,is3max do 1000 k4 = 1,is4max do 1000 k5 = 1,is5max do 1000 k6 = 1,is6max xg = stu(k3,3) if(k4.eq.1.and.k5.ge.is5max-1) goto 1000 if(k4.eq.1.and.k6.ge.is6max-1) goto 1000 if(k5.eq.1.and.k4.ge.is4max-1) goto 1000 if(k5.eq.1.and.k6.ge.is6max-1) goto 1000 if(k6.eq.1.and.k4.ge.is4max-1) goto 1000 if(k6.eq.1.and.k5.ge.is5max-1) goto 1000 if(k4.eq.2.and.k5.ge.is5max) goto 1000 if(k4.eq.2.and.k6.ge.is6max) goto 1000 if(k5.eq.2.and.k4.ge.is4max) goto 1000 if(k5.eq.2.and.k6.ge.is6max) goto 1000 if(k6.eq.2.and.k4.ge.is4max) goto 1000 if(k6.eq.2.and.k5.ge.is5max) goto 1000 u2 = stu(k4,4) u3 = stu(k5,5) u4 = stu(k6,6) numpar = 5 tr2 = u2 tr3 = u3 tr4 = u4 xsol2(1) = 1./sqrt(2.*3.141592654*u2) xsol2(2) = 0. xsol2(3) = -0.5/u2 xsol2(4) = 0. xsol2(5) = 0. read(99,'(5f14.6)') (xsol2(j),j=1,5) call newt(xsol2,numpar,check) ein(1,k3,k4,k5,k6) = xsol2(1) ein(2,k3,k4,k5,k6) = xsol2(2) ein(3,k3,k4,k5,k6) = xsol2(3) ein(4,k3,k4,k5,k6) = xsol2(4) ein(5,k3,k4,k5,k6) = xsol2(5) write(*,'(5f14.6)') (xsol2(j),j=1,5) 1000 continue * * * Iterating on Aggregate Policy Rules * * do 448 itera = 1,100 write(*,*) itera,448 * * * Iterating on S(beta) * do 300 iter = 1,10 call sbeta() 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 aold(i) = eps1*anew(i) + (1.-eps1)*aold(i) 310 continue c (*,'(6f13.6)') ( aold(i),i=1,npar) c write(*,'(6f13.6)') ( anew(i),i=1,npar) write(*,*) 'indiv' , accur if(accur.lt.0.05) goto 305 300 continue 305 continue * * * UPDATE THE COEFFICIENTS IN THE AGGREGATE POLICY RULES * * call update() 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 aggold(i,j) = eps2*aggnew(i,j) + (1.-eps2)*aggold(i,j) 3310 continue write(*,*) 'macro', accur if(accur.lt.0.001) goto 3305 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 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,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ 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/ hx(1000),hw(1000) common/solut/ aold(npar),anew(npar) common/istate/ k3,k4,k5,k6 common/trunc/ tr1,tr2,tr3,tr4 common/ifu/ ifunc common/den/ e1,e2,e3,e4,e5 common/ddddd/ ein(5,5,5,5,5) common/is/ ispoint2(5,5,5,5),jt2max dimension xsol(1),xsol2(5),xx(npar2,npar2),xy(500),yreg(4000,4) dimension sim1(1000),sim2(1000) 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 if(k4.eq.1.and.k5.ge.is5max-1) goto 1000 if(k4.eq.1.and.k6.ge.is6max-1) goto 1000 if(k5.eq.1.and.k4.ge.is4max-1) goto 1000 if(k5.eq.1.and.k6.ge.is6max-1) goto 1000 if(k6.eq.1.and.k4.ge.is4max-1) goto 1000 if(k6.eq.1.and.k5.ge.is5max-1) goto 1000 if(k4.eq.2.and.k5.ge.is5max) goto 1000 if(k4.eq.2.and.k6.ge.is6max) goto 1000 if(k5.eq.2.and.k4.ge.is4max) goto 1000 if(k5.eq.2.and.k6.ge.is6max) goto 1000 if(k6.eq.2.and.k4.ge.is4max) goto 1000 if(k6.eq.2.and.k5.ge.is5max) goto 1000 e1 = ein(1,k3,k4,k5,k6) e2 = ein(2,k3,k4,k5,k6) e3 = ein(3,k3,k4,k5,k6) e4 = ein(4,k3,k4,k5,k6) e5 = ein(5,k3,k4,k5,k6) jt1 = ispoint2(k3,k4,k5,k6) c write(*,*) jt1 c write(*,*) k3,k4,k5,k6 xg = stu(k3,3) * * 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 * * c c note that the mean equals zero c second = 0. third = 0. fourth = 0. hh = bdiff/dble(nsimpson-1) do 1300 k = 1,nincome do 1350 j = 1,nsimpson hin = hx(j) if(hx(j).lt.blow/xg) hin = blow/xg if(hx(j).gt.bhigh/xg) hin = bhigh/xg wpoor = hin + qx(k) + stu(1,2) indiv = 1 bpoor = bfunc(wpoor,indiv,k3,k4,k5,k6) wrich = hin + qx(k) + stu(2,2) indiv = 2 brich = bfunc(wrich,indiv,k3,k4,k5,k6) sim1(j) = bpoor sim2(j) = brich 1350 continue c c c don't forget that this only works when individual prob = 0.5 c c do 1360 j = 1,nsimpson dd = hh*dens(e1,e2,e3,e4,e5,hx(j))*hw(j)*qw(k) second = sim1(j)**2.*dd*prob(1,1) + $ sim2(j)**2.*dd*prob(1,2) + second third = sim1(j)**3.*dd*prob(1,1) + $ sim2(j)**3.*dd*prob(1,2) + third fourth = sim1(j)**4.*dd*prob(1,1) + $ sim2(j)**4.*dd*prob(1,2) + fourth 1360 continue 1300 continue c write(*,'(3f14.6)') second,third,fourth yreg(jt1,2) = sca(second,4) yreg(jt1,3) = sca(third,5) yreg(jt1,4) = sca(fourth,6) 1000 continue * * * do the regression * * do 5000 ireg = 1,4 JMAX = jt2max NPMAX = NPAR2 DO 110 J = 1,NPMAX XY(J) = 0.0 DO 110 K = 1,NPMAX XX(J,K) = 0.0 110 CONTINUE DO 130 K = 1,NPMAX DO 130 J = 1,NPMAX DO 130 I = 1,JMAX XX(J,K) = XX(J,K) + CHEB2(I,J)*CHEB2(I,K) 130 CONTINUE c write(*,'(20f4.1)') ((xx(ii,jj),ii=1,npmax),jj=1,npmax) CALL INVERT(XX,NPMAX,DET) 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) = 0. DO 145 K = 1,NPMAX AGGNEW(IREG,J) = AGGNEW(IREG,J)+XX(J,K)*XY(K) 145 CONTINUE c DO 150 I = 1,6 c WRITE(*,*) I,AGGNEW(IREG,I) c150 CONTINUE 5000 continue end 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,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ 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/ hx(1000),hw(1000) common/solut/ aold(npar),anew(npar) common/istate/ k3,k4,k5,k6 common/den/ e1,e2,e3,e4,e5 hh = bdiff/dble(nsimpson-1) price(1) = qq arich = 0. indiv = 2 do 930 i = 1,nsimpson do 930 k = 1,nincome hin = hx(i) if(hx(i).lt.blow/xg) hin = blow/xg if(hx(i).gt.bhigh/xg) hin = bhigh/xg wealth= hin + qx(k) + stu(indiv,2) arich = arich + bfunc(wealth,indiv,k3,k4,k5,k6)*hw(i)*qw(k)* $ hh*dens(e1,e2,e3,e4,e5,hx(i)) 930 continue apoor = 0. indiv = 1 do 830 i = 1,nsimpson do 830 k = 1,nincome hin = hx(i) if(hx(i).lt.blow/xg) hin = blow/xg if(hx(i).gt.bhigh/xg) hin = bhigh/xg wealth = hin + qx(k) + stu(indiv,2) apoor = apoor + bfunc(wealth,indiv,k3,k4,k5,k6)*hw(i)*qw(k)* $ hh*dens(e1,e2,e3,e4,e5,hx(i)) 830 continue usrf = arich + apoor end subroutine funcv(n,x,fvec) implicit real*8(a-h,o-z),integer(i-n) parameter(npar=50) common/rpar/ xg,blow,bhigh,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome common/hermite/ hx(1000),hw(1000) common/trunc/ tr1,tr2,tr3,tr4 dimension x(n),fvec(n),hermx(30) e1 = x(1) e2 = x(2) e3 = x(3) e4 = x(4) e5 = x(5) a1 = 0. a2 = 0. a3 = 0. a4 = 0. a0 = 0. hh = bdiff/dble(nsimpson-1) do 930 i = 1,nsimpson hin = hx(i) if(hin.lt.blow/xg) hin = blow/xg if(hin.gt.bhigh/xg) hin = bhigh/xg a0 = a0 + hh* dens(e1,e2,e3,e4,e5,hx(i))*hw(i) a1 = a1 + hh*hin *dens(e1,e2,e3,e4,e5,hx(i))*hw(i) a2 = a2 + hh*hin**2.*dens(e1,e2,e3,e4,e5,hx(i))*hw(i) a3 = a3 + hh*hin**3.*dens(e1,e2,e3,e4,e5,hx(i))*hw(i) a4 = a4 + hh*hin**4.*dens(e1,e2,e3,e4,e5,hx(i))*hw(i) 930 continue fvec(2) = a1 - 0. fvec(3) = (a2 - tr2)*10. fvec(4) = (a3 - tr3)*100. fvec(5) = (a4 - tr4)*100. fvec(1) = a0 - 1. c write(*,'(5f14.6)') (fvec(j),j=1,5) c write(*,'(5f14.6)') a0,a1,a2,a3,a4 end 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,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ 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)**2.) s4 = sca(s4,4) ieq = 3 s5 = aggfunc(s3old,s4old,s5old,s6old,ieq) s5 = scinv(s5,5) s5 = s5/(stu(lnew,3)**3.) s5 = sca(s5,5) ieq = 4 s6 = aggfunc(s3old,s4old,s5old,s6old,ieq) s6 = scinv(s6,6) s6 = s6/(stu(lnew,3)**4.) 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 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 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,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome sca = ( 2.*x - stbar(i,2) - stbar(i,1) ) $ / (stbar(i,2)-stbar(i,1)) end 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,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ is1max,is2max,is3max,is4max,is5max,is6max,imax,nincome scinv = ( stbar(i,2) + stbar(i,1) $ + x*(stbar(i,2) - stbar(i,1)) ) / 2. end * * * Define Chebyshev Polynomials * * 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 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 integer function ipoint(i1,i2,i3) implicit real*8(a-h,o-z),integer(i-n) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ 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 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,nsimpson, $ 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 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,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ 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/ hx(1000),hw(1000) 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 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,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ 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 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,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ 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 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,bdiff,bondp,beta,gam,price(10) common/ipar/ i1max,i2max,i3max,i4max,i5max,i6max,nsimpson, $ 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/ hx(1000),hw(1000) 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 real*8 function dens(e1,e2,e3,e4,e5,bb) implicit real*8(a-h,o-z),integer(i-n) dens = e1*exp(e2*bb+e3*bb**2.+e4*bb**3.+e5*bb**4.) end 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=. SUBROUTINE INVERT(A,N,D) C C C THIS SUBROUTINE COMPUTES THE INVERSE OF A MATRIX. C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(N,N) ,L(400) , M(400) COMMON // ISING D=1.D0 DO 80 K=1,N L(K)=K M(K)=K BIGA=A(K,K) DO 20 I=K,N DO 20 J=K,N IF(DABS(BIGA)-DABS(A(I,J))) 10,20,20 10 BIGA=A(I,J) L(K)=I M(K)=J 20 CONTINUE IF (DABS(BIGA).LT.(1.0E-15)) GO TO 99 J=L(K) IF(L(K)-K) 35,35,25 25 DO 30 I=1,N HOLD=-A(K,I) A(K,I)=A(J,I) 30 A(J,I)=HOLD 35 I=M(K) IF(M(K)-K) 45,45,37 37 DO 40 J=1,N HOLD=-A(J,K) A(J,K)=A(J,I) 40 A(J,I)=HOLD 45 DO 55 I=1,N IF(I-K) 50,55,50 50 A(I,K)=A(I,K)/(-A(K,K)) 55 CONTINUE DO 65 I=1,N DO 65 J=1,N IF(I-K) 57,65,57 57 IF(J-K) 60,65,60 60 A(I,J)=A(I,K)*A(K,J)+A(I,J) 65 CONTINUE DO 75 J=1,N IF(J-K) 70,75,70 70 A(K,J)=A(K,J)/A(K,K) 75 CONTINUE D=D*A(K,K) A(K,K)=1.D0/A(K,K) 80 CONTINUE K=N 100 K=K-1 IF(K) 150,150,103 103 I=L(K) IF(I-K) 120,120,105 105 DO 110 J=1,N HOLD=A(J,K) A(J,K)=-A(J,I) 110 A(J,I)=HOLD 120 J=M(K) IF(J-K) 100,100,125 125 DO 130 I=1,N HOLD=A(K,I) A(K,I)=-A(J,I) 130 A(J,I)=HOLD GO TO 100 99 CONTINUE ISING = 1 WRITE(11,991) ISING 991 FORMAT(' SINGULAR MATRIX = = > ISING = ',I4) 150 CONTINUE RETURN END 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