C **************************************************************** c c c simulation program for xx5.f c c January 2, 1995 C **************************************************************** implicit real*8(A-H,O-Z),integer(i-n) common/series/ y(100,26000),bj(100),bjn(100), $ q(26000),ya(26000),ww(100), $ cj(100),polc(2,2,3,25,11,1),polb(2,2,3,25,11,1), $ polq(2,2,3,25,11,1), $ prede(100),s1(2),s2(20),i3(100),s4(30),s5(20),i6(10), $ s3(100),s6(100) common/par/ bound,nper,pagg,pbu,pgu,ye,yu,del,gam, $ m1,m2,m3,m4,m5,m6,kk,m1t,m2t,m3t,m4t,m5t,m6t common/ser/ jmax(10),yy(100,10) dimension ex(2,2),accmax(3), $ intacc(3,6),c1(26000),b1(26000),ca(26000),stan(26000) dimension m2j(100),i3j(100),s4j(100),s5j(100),i6j(100), $ prob(2,2,2,3),xinter(10),iun(26000),vmean(100),var(100) external func open(33,file='todd.par',status='old') open(34,file='todd.in',status='old') read(33,*) bound,acc,del,nper,idum,numiter,gam if(nper.gt.100) write(*,*) 'nper too big, adjust dimensions' if(nper.gt.100) stop nobs = 2500 idum = -123194 fl = 0.9 fu = 1.1 facc = 0.00001 ye = 100. yu = 50. pb0 = 0.5 pb1 = 0.4 pb2 = 0.1 pg0 = 0.65 pg1 = 0.3 pg2 = 0.05 pagg = 0.9375 c c c c prob(i,j,k,l) is the probability that for next period c c m1 = j, m2 = k and m3 = l given that m1 = i today c c c prob(1,1,1,1) = pagg*pg0*0. prob(1,1,2,1) = pagg*pg0*1. prob(1,1,1,2) = pagg*pg1*0.1 prob(1,1,2,2) = pagg*pg1*0.9 prob(1,1,1,3) = pagg*pg2*0.2 prob(1,1,2,3) = pagg*pg2*0.8 prob(1,2,1,1) = (1.-pagg)*pb0*0. prob(1,2,2,1) = (1.-pagg)*pb0*1. prob(1,2,1,2) = (1.-pagg)*pb1*0.1 prob(1,2,2,2) = (1.-pagg)*pb1*0.9 prob(1,2,1,3) = (1.-pagg)*pb2*0.2 prob(1,2,2,3) = (1.-pagg)*pb2*0.8 prob(2,1,1,1) = (1.-pagg)*pg0*0. prob(2,1,2,1) = (1.-pagg)*pg0*1. prob(2,1,1,2) = (1.-pagg)*pg1*0.1 prob(2,1,2,2) = (1.-pagg)*pg1*0.9 prob(2,1,1,3) = (1.-pagg)*pg2*0.2 prob(2,1,2,3) = (1.-pagg)*pg2*0.8 prob(2,2,1,1) = pagg*pb0*0. prob(2,2,2,1) = pagg*pb0*1. prob(2,2,1,2) = pagg*pb1*0.1 prob(2,2,2,2) = pagg*pb1*0.9 prob(2,2,1,3) = pagg*pb2*0.2 prob(2,2,2,3) = pagg*pb2*0.8 c c c setting up grid c c s1: aggregate state c s2: individual income status c i3: number unemployed agents c s4: individual bond level c s5: pseudo sd c m1t = 2 m2t = 2 m3t = 3 m4t =25 m5t =11 s1(1) = 10 s1(2) = -10 s2(1) = yu s2(2) = ye i3(1) = 0 i3(2) = 1 i3(3) = 2 do 22 jj = 1,m4t s4(jj) = bound+2.*abs(bound)*dble(jj-1)/dble(m4t-1) write(*,*) jj,s4(jj) 22 continue do 23 jj = 1,m5t s5(jj) = dble(jj)*5. 23 continue c c c initializing the stuff for the cons and bond interpolation c c jmax(1) = m4t jmax(2) = m5t do 206 ii = 1,m4t yy(ii,1) = s4(ii) 206 continue do 207 ii = 1,m5t yy(ii,2) = s5(ii) 207 continue c c c do 209 jj = 1,10 xinter(jj) = 0.1*dble(jj-1) 209 continue c c c initializing the policy functions c c do 440 m1 = 1,m1t do 440 m2 = 1,m2t do 440 m3 = 1,m3t do 440 m4 = 1,m4t do 440 m5 = 1,m5t read(34,*) ij1,ij2,ij3,ij4,ij5, $ polc(m1,m2,m3,m4,m5,1), $ polb(m1,m2,m3,m4,m5,1), $ polq(m1,m2,m3,m4,m5,1) 440 continue write(*,*) 'done reading' c c generate exogenous process c ya(1) = 10 do 100 i = 2,nobs xxx = ran1(idum) if(xxx.le.pagg) then ya(i) = ya(i-1) else ya(i) = -ya(i-1) endif 100 continue c c c !!! iun(i) is not the NUMBER of unemployed c but the value of the third state variable c is iun=1, then # unemployed = 0 etc. c do 110 i = 1,nobs xxx = ran1(idum) if(ya(i).gt.0) then if(xxx.le.0.65) iun(i) = 1 if(xxx.gt.0.65.and.xxx.le.0.95) iun(i) = 2 if(xxx.gt.0.95) iun(i) = 3 else if(xxx.le.0.5) iun(i) = 1 if(xxx.gt.0.5.and.xxx.le.0.9) iun(i) = 2 if(xxx.gt.0.9) iun(i) = 3 endif 110 continue do 120 i = 1,nobs do 130 j = 1,nper y(j,i) = ye 130 continue if(iun(i).eq.1) goto 120 do 123 j = 1,nper/10 1234 continue xxx1 = ran1(idum) xxx1 = xxx1*dble(nper) jun = int(xxx1) + 1 if(y(jun,i).eq.yu) goto 1234 y(jun,i) = yu 123 continue if(iun(i).eq.2) goto 119 do 124 j = 1,nper/10 12345 continue xxx1 = ran1(idum) xxx1 = xxx1*dble(nper) jun = int(xxx1) + 1 if(y(jun,i).eq.yu) goto 12345 y(jun,i) = yu 124 continue 119 continue 120 continue do 140 j = 1,nper/2 bj(j) = 10 140 continue do 141 j = 1+nper/2,nper bj(j) = -10. 141 continue do 510 i = 1,nobs if(ya(i).ge.0.) then m1 = 1 else m1 = 2 endif do 1390 j = 1,nper vmean(j) = -bj(j)/dble(nper-1) 1390 continue do 1395 j = 1,nper var(j) = 0.0 do 1395 jk = 1,nper if(j.ne.jk) var(j) = var(j) + (bj(jk)-vmean(j))**2. 1395 continue do 1397 j = 1,nper var(j) = (var(j)/dble(nper-1))**0.5 1397 continue stan(i) = 0.0 do 1399 j = 1,nper stan(i) = stan(i) + bj(j)**2. 1399 continue stan(i) = ( stan(i)/dble(nper) )**0.5 do 420 j = 1,nper ww(j) = bj(j) + y(j,i) s5j(j) = var(j) c c define the state variables for the individual c if(y(j,i).eq.yu) then m2j(j) = 1 else m2j(j) = 2 endif bjn(j) = pol(m1,m2j(j),iun(i),bj(j),s5j(j),2) 420 continue c c c use this to calculate the conditional expectation c c c c c first calculate next period's 5th state variable c c do 2390 j = 1,nper vmean(j) = -bjn(j)/dble(nper-1) 2390 continue do 2395 j = 1,nper var(j) = 0.0 do 2395 jk = 1,nper if(j.ne.jk) var(j) = var(j) + (bjn(jk)-vmean(j))**2. 2395 continue do 2397 j = 1,nper var(j) = (var(j)/dble(nper-1))**0.5 2397 continue do 421 j = 1,nper prede(j) = 0.0 do 20 k1 = 1,m1t do 20 k2 = 1,m2t do 20 k3 = 1,m3t aaa = ( pol(k1,k2,k3,bjn(j),var(j),1) )**gam prede(j) = prede(j)+prob(m1,k1,k2,k3)*aaa 20 continue prede(j) = del*prede(j) 421 continue equil = rtsec(func,fl,fu,facc) c c c use KT-conditions and budget constraint to calculate c consumption and bond-hodings c c cc = ( prede(1) / equil ) ** (1./gam) bjn(1) = ( y(1,i)+bj(1)-cc ) / equil if(bjn(1).lt.bound) then bjn(1) = bound cc = y(1,i) + bj(1) - equil*bjn(1) endif c1(i) = cc b1(i) = bjn(1) bj(1) = bjn(1) do 430 j = 2,nper cc = ( prede(j) / equil ) ** (1./gam) bjn(j) = ( y(j,i)+bj(j)-cc ) / equil if(bjn(j).lt.bound) bjn(j) = bound bj(j) = bjn(j) 430 continue q(i) = equil 510 continue c c c CALCULATING SOME SAMPLE STATISTICS c c cmean = 0.0 cstan = 0.0 camean = 0.0 castan = 0.0 rmean = 0.0 rstan = 0.0 ymean = 0.0 ystan = 0.0 bmean = 0.0 bstan = 0.0 stanmean = 0.0 stanstan = 0.0 do 710 i = 2,nobs ca(i) = 0.0 do 711 j = 1,nper ca(i) = ca(i) + y(j,i) 711 continue ca(i) = ca(i)/dble(nper) 710 continue do 720 i = 101,nobs stanmean = stanmean + stan(i) cmean = cmean + log(c1(i)/c1(i-1)) camean = camean + log(ca(i)/ca(i-1)) rmean = rmean - log(q(i)) ymean = ymean + log(y(1,i)/y(1,i-1)) bmean = bmean + b1(i) 720 continue cmean = cmean/dble(nobs-100) camean = camean/dble(nobs-100) rmean = rmean/dble(nobs-100) ymean = ymean/dble(nobs-100) bmean = bmean/dble(nobs-100) stanmean = stanmean/dble(nobs-100) do 730 i = 101,nobs cstan = cstan + (log(c1(i)/c1(i-1)) - cmean)**2. castan = castan + (log(ca(i)/ca(i-1)) - camean)**2. rstan = rstan + (-log(q(i))-rmean)**2. ystan = ystan + (log(y(1,i)/y(1,i-1))-ymean)**2. bstan = bstan + (b1(i) - bmean)**2. stanstan = stanstan + (stan(i)-stanmean)**2. 730 continue cstan = sqrt(cstan/dble(nobs-100) ) castan = sqrt(castan/dble(nobs-100) ) rstan = sqrt(rstan/dble(nobs-100) ) ystan = sqrt(ystan/dble(nobs-100) ) bstan = sqrt(bstan/dble(nobs-100) ) stanstan = sqrt(stanstan/dble(nobs-100) ) write(*,*) write(*,*) write(*,509) del,-log(del)*100.,bound,gam write(*,*) write(*,516) ikeep,nobs,nper write(*,*) write(*,515) m1t,m2t,m3t,m4t,m5t,m6t write(*,508) s2(1),i3(1),s4(1),s5(1),i6(1) write(*,507) s2(m2t),i3(m3t),s4(m4t),s5(m5t),i6(m6t) write(*,*) write(*,*) write(*,511) cmean,cstan write(*,512) ymean,ystan write(*,513) camean,castan write(*,514) rmean,rstan write(*,517) bmean,bstan write(*,518) stanmean,stanstan 500 format(' bet :',6f13.6) 502 format(' not solved bond price iteration at obs.: ',i6) 509 format(' del = ',f6.3,' discount rate (%) = ',f6.3, $ ' bound = ',f7.2,' gam = ',f6.2) 517 format(' agent 1s bond holdings: ',2f10.6) 518 format(' cross-sectional stan dev: ',2f10.6) 511 format(' agent 1s dlog consumption: ',2f10.6) 512 format(' agent 1s dlog income: ',2f10.6) 513 format(' dlog aggregate consumption: ',2f10.6) 514 format(' interest rate: ',2f10.6) 516 format(' idum = ',i10,' nobs = ',i5,' nper =',i5) 515 format(' gridpoint: ',i8,2i6,i7,i9,i7) 507 format(' minimum grid points: ',f6.2,i4,2f9.2,i4) 508 format(' maximum grid points: ',f6.2,i4,2f9.2,i4) stop end real*8 function func(qq) implicit real*8(a-h,o-z),integer(i-n) common/series/ y(100,26000),bj(100),bjn(100), $ q(26000),ya(26000),ww(100), $ cj(100),polc(2,2,3,25,11,1),polb(2,2,3,25,11,1), $ polq(2,2,3,25,11,1), $ prede(100),s1(2),s2(20),i3(100),s4(30),s5(20),i6(10), $ s3(100),s6(100) common/par/ bound,nper,pagg,pbu,pgu,ye,yu,del,gam, $ m1,m2,m3,m4,m5,m6,kk,m1t,m2t,m3t,m4t,m5t,m6t dimension save(100) do 220 j = 1,nper cj(j) = (prede(j)/qq)**(1./gam) save(j) = ( ww(j) - cj(j) ) / qq if(save(j).lt.bound) then save(j) = bound endif 220 continue c c c check aggregate demand conditions c c btot = 0.0 do 235 j = 1,nper btot = btot + save(j) 235 continue func = btot end real*8 function pol(m1p,m2p,m3p,s4p,s5p,kkk) implicit real*8(a-h,o-z),integer(i-n) common/series/ y(100,26000),bj(100),bjn(100), $ q(26000),ya(26000),ww(100), $ cj(100),polc(2,2,3,25,11,1),polb(2,2,3,25,11,1), $ polq(2,2,3,25,11,1), $ prede(100),s1(2),s2(20),i3(100),s4(30),s5(20),i6(10), $ s3(100),s6(100) common/par/ bound,nper,pagg,pbu,pgu,ye,yu,del,gam, $ m1,m2,m3,m4,m5,m6,kk,m1t,m2t,m3t,m4t,m5t,m6t common/ser/ jmax(10),yy(100,10) dimension z3(2,2,2),z2(2,2),z1(2), $ xx(4),t(4),ll(4,2),z4(2,2,2,2) call locate(s4,m4t,s4p,j1) call locate(s5,m5t,s5p,j2) ll(1,1) = j1 ll(2,1) = j2 xx(1) = s4p xx(2) = s5p do 10 ii = 1,2 if(ll(ii,1).eq.0) then ll(ii,1) = 1 endif if(ll(ii,1).eq.jmax(ii)) then ll(ii,1) = ll(ii,1) -1 endif ll(ii,2) = ll(ii,1) + 1 t(ii) = (xx(ii)-yy(ll(ii,1),ii)) $ / (yy(ll(ii,2),ii)-yy(ll(ii,1),ii)) 10 continue do 20 l1 = 1,2 do 20 l2 = 1,2 if(kkk.eq.1) z4(l1,l2,1,1) $ = polc(m1p,m2p,m3p,ll(1,l1),ll(2,l2),1 ) if(kkk.eq.2) z4(l1,l2,1,1) $ = polb(m1p,m2p,m3p,ll(1,l1),ll(2,l2),1 ) if(kkk.eq.3) z4(l1,l2,1,1) $ = polq(m1p,m2p,m3p,ll(1,l1),ll(2,l2),1 ) 20 continue do 30 l2 = 1,2 z3(l2,1,1) = (1.-t(1))*z4(1,l2,1,1) $ + t(1) *z4(2,l2,1,1) 30 continue pol = (1.-t(2))*z3(1,1,1) $ + t(2) *z3(2,1,1) if(kkk.eq.2) then if(pol.lt.bound) pol = bound endif end subroutine locate(xx,n,x,j) implicit real*8(a-h,o-z),integer(i-n) dimension xx(n) jl = 0 ju = n+1 10 if(ju-jl.gt.1) then jm=(ju+jl)/2 if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then jl=jm else ju=jm endif goto 10 endif j=jl return end real*8 function perm(j,k,nper) implicit real*8(a-h,o-z),integer(i-n) if(j.eq.0) then perm = 1 goto 100 endif s2 = 1 s3 = 1 if(j.ge.k) then nhigh = j nlow = k else nhigh = k nlow = j endif do 10 i = nhigh+1,nper s3 = s3*dble(i) 10 continue do 20 i = 1,nlow s2 = s2*dble(i) 20 continue perm = s3/s2 100 continue end C C C C ***************************************** C C THE FOLLOWING TWO SUBROUTINES GENERATE RANDOM NUMBERS C C ***************************************** C REAL*8 FUNCTION GASDEV(IDUM) c c retruns a normally distributed deviate with zero mean and unit c variance, using RAN1(IDUM) as the source of uniform deviates IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N) save iset,gset DATA ISET/0/ IF (ISET.EQ.0) THEN 1 V1=2.*RAN1(IDUM)-1. V2=2.*RAN1(IDUM)-1. R=V1**2.+V2**2. IF(R.GE.1..OR.R.EQ.0.)GO TO 1 FAC=SQRT(-2.*LOG(R)/R) GSET=V1*FAC GASDEV=V2*FAC ISET=1 ELSE GASDEV=GSET ISET=0 ENDIF RETURN END REAL*8 FUNCTION RAN1(IDUM) c c returns a uniform random deviae between 0.0 and 1.0. c set IDUM to any negative value to initialize or reinitialize c the sequence IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N) DIMENSION R(97) PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=1./M1) PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=1./M2) PARAMETER (M3=243000,IA3=4561,IC3=51349) save iff,ix1,ix2,ix3,r DATA IFF /0/ IF(IDUM.LT.0.OR.IFF.EQ.0) THEN IFF = 1 IX1 = MOD(IC1-IDUM,M1) IX1 = MOD(IA1*IX1+IC1,M1) IX2 = MOD(IX1,M2) IX1 = MOD(IA1*IX1+IC1,M1) IX3 = MOD(IX1,M3) DO 11 J=1,97 IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) R(J)=(DBLE(IX1)*DBLE(IX2)*RM2)*RM1 11 CONTINUE IDUM = 1 ENDIF IX1=MOD(IA1*IX1+IC1,M1) IX2=MOD(IA2*IX2+IC2,M2) IX3=MOD(IA3*IX3+IC3,M3) J=1+(97*IX3)/M3 IF(J.GT.97.OR.J.LT.1)PAUSE RAN1=R(J) R(J)=(DBLE(IX1)+DBLE(IX2)*RM2)*RM1 RETURN END SUBROUTINE SORT(N,RA) implicit real*8(a-h,o-z) DIMENSION RA(N) L=N/2+1 IR=N 10 CONTINUE IF(L.GT.1)THEN L=L-1 RRA=RA(L) ELSE RRA=RA(IR) RA(IR)=RA(1) IR=IR-1 IF(IR.EQ.1)THEN RA(1)=RRA RETURN ENDIF ENDIF I=L J=L+L 20 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(RA(J).LT.RA(J+1))J=J+1 ENDIF IF(RRA.LT.RA(J))THEN RA(I)=RA(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF RA(I)=RRA GO TO 10 END REAL*8 FUNCTION RTBIS(FUNC,X1,X2,XACC) IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N) PARAMETER (JMAX=100) EXTERNAL FUNC FMID=FUNC(X2) F=FUNC(X1) IF(F*FMID.GE.0.) PAUSE 'Root must be bracketed for bisection.' IF(F.LT.0.)THEN RTBIS=X1 DX=X2-X1 ELSE RTBIS=X2 DX=X1-X2 ENDIF DO 11 J=1,JMAX DX=DX*.5 XMID=RTBIS+DX FMID=FUNC(XMID) IF(FMID.LE.0.)RTBIS=XMID IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN 11 CONTINUE PAUSE 'too many bisections' END REAL*8 FUNCTION RTSEC(FUNC,X1,X2,XACC) IMPLICIT REAL*8(A-H,O-Z),INTEGER(I-N) PARAMETER (MAXIT=30) EXTERNAL FUNC 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