c c XXPE6.F, PEA PROGRAM USING PERCENTILES c c DECEMBER 1994 c c This program only runs with nper = 10. Running it with c a different number of agents would require some small c changes in generating the unemployed agents, setting the c percentiles and the dimensions c c **************************************************************** implicit real*8(A-H,O-Z) common y(10,51000),bj(100),w1(51000),y1(51000),b1(51000), $ wdist(51000,10),agg(51000), $ q(51000),ya(51000),bet(15),sbet(15), $ cj(100),fru(51000),c1(51000), $ umar(51000),temp(51000),bold(100),bound,nper,np,num,i dimension ca(51000),ww(1000),iun(51000),xinter(10) external func 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 open(31,file='pe.in',status='old') open(32,file='pe.new',status='unknown') open(33,file='pe.par',status='old') c c c example for pe.par c c c c 0.001 -40. 0.0001 0.96 .7 15 50000 10 -333921 500 200 5 1 3 7 10 c c corresponding pe.in c c 1.142463663467267E-02 -5.489783219680950E-02 -.18114419565413 c -1.437471088104099E-02 .1146736033184872 -7.643164257621127E-04 c -3.295828620506037E-03 -7.002972348455280E-03 -2.430356898619899E-02 c .1937518453431341 .1625948316166882 4.607839556308447E-02 c 4.521787375124224E-03 -.104368380698151 .1017658256177863 c c generating the exogenous income process c note that the storage of y(10,51000) requires a lot of space c in principle we do not have to do this since by setting the c iseed within the SBETA do-loop we also can get the same c random numbers but this will slow down the iterations c read(33,*) frac,bound,acc,del,eps,np,nobs,nper,idum,numiter, $ iterq,num,iz1,iz2,iz3,iz4 ikeep = idum iz1 = 1 iz2 = 3 iz3 = 7 iz4 = 10 read(31,*) (bet(j),j=1,np) do 209 jj = 1,10 xinter(jj) = 0.1*dble(jj-1) 209 continue 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 if(iun(i).eq.2) then xxx1 = ran1(idum) call locate(xinter,nper,xxx1,junem) y(junem,i) = yu else xxx1 = ran1(idum) call locate(xinter,nper,xxx1,junem1) y(junem1,i) = yu 125 continue xxx2 = ran1(idum) call locate(xinter,nper,xxx2,junem2) if(junem2.eq.junem1) goto 125 y(junem2,i) = yu endif 120 continue do 1110 i = 1,nobs ya(i) =0.01*( ya(i) + 110. ) agg(i) = 0.0 do 1125 j = 1,nper agg(i) = agg(i) + y(j,i) 1125 continue agg(i) = 0.01*agg(i)/dble(nper) 1110 continue c c do 180 j = 1,np sbet(j) = bet(j) 180 continue do 300 iter = 1,numiter c c c generating the series c c do 205 j = 1,nper bj(j) = 0.0 205 continue do 210 i = 2,nobs c c c calculating information about the wealth distribution c c do 140 j= 1,nper ww(j) = y(j,i)+bj(j) 140 continue call sort(nper,ww) wdist(i,1) = ( 0.01*( ww(iz1))) wdist(i,2) = ( 0.01*( ww(iz2))) wdist(i,3) = ( 0.01*( ww(iz3))) wdist(i,4) = ( 0.01*( ww(iz4))) b1(i) = bj(1) w1(i) = ( 0.01*(y(1,i) + bj(1) )) xxx = bet(1)*exp( bet(2)*ya(i) $ + bet(4)*wdist(i,3) + bet(5)*agg(i) $ + bet(7)*wdist(i,1) + bet(8)*wdist(i,2) $ + bet(9)*wdist(i,4) + $ bet(13)*agg(i)*ya(i)) do 214 j = 1,nper ww(j) = 0.01*(y(j,i)+bj(j) ) bbb = bj(j)/100. temp(j) = bet(3)*ww(j) + bet(6)*ww(j)*y(j,i) + $ bet(10)*bbb + bet(11)*bbb**2. + bet(12)*bbb**3. + $ bet(14)*bbb*agg(i) + bet(15)*y(j,i)*agg(i)/100. temp(j) = xxx*exp(temp(j)) bold(j) = bj(j) 214 continue c c q(i) = rtsec(func,fl,fu,facc) do 220 j = 1,nper cj(j) = q(i)/temp(j) bj(j) = ( y(j,i) + bold(j) - cj(j) ) / q(i) if(bj(j).lt.bound) then bj(j) = bound cj(j) = y(j,i) + bold(j) - bound*q(i) endif 220 continue umar(i) = del/cj(1) c1(i) = cj(1) 210 continue c c calculating sbeta c c write(*,*) ' start sbeta ' call sbeta(nobs) c c c updating beta c c write(*,*) write(*,*) iter write(*,500) (bet(i), i=1,5 ) write(*,500) (bet(i), i=6,10) write(*,500) (bet(i), i=11,15) diff = 0. do 310 i=1,np diff = diff + abs( ( bet(i)-sbet(i) ) / bet(i) ) 310 bet(i) = bet(i)*(1-eps) + sbet(i)*eps write(*,500) (sbet(i), i=1,5 ) write(*,500) (sbet(i), i=6,10) write(*,500) (sbet(i), i=11,15) write(*,'(a7,f10.5)') 'diff= ' , diff if ( diff.lt.acc ) $ go to 400 300 continue write(*,'(a)') ' iterations on beta did not converge' write(*,'(a)') ' -----------------------------------' 400 continue write(32,*) (bet(j),j=1,np) do 610 i = 501,1000 write(12,'(i6,6f11.5)') i,c1(i),y(1,i),b1(i),ya(i),q(i) 610 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 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 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) 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. 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) ) write(*,*) write(*,*) write(*,510) del,-log(del)*100.,bound write(*,*) write(*,516) np,nobs,nper write(*,*) write(*,515) frac,ikeep,iterq write(*,*) write(*,*) write(*,511) cmean,cstan write(*,512) ymean,ystan write(*,513) camean,castan write(*,514) rmean,rstan write(*,517) bmean,bstan 500 format(' bet :',6f13.6) 502 format(' not solved bond price iteration at obs.: ',i6) 510 format(' del = ',f6.3,' discount rate (%) = ',f6.3, $ ' bound = ',f7.2) 511 format(' agent 1s dlog consumption: ',2f10.6) 517 format(' agent 1s bond holdings: ',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(' np = ',i3,' nobs = ',i5,' nper =',i5) 515 format(' frac = ',f8.6,' idum = ',i8,' iterq = ',i5) stop end C C C SUBROUTINE INVERT(A,N,D) C IMPLICIT REAL*8(A-H,O-Z) REAL*8 A(N,N) , L(20) , M(20) common /ss/ ISING D=1. 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(ABS(BIGA)-ABS(A(I,J))) 10,20,20 10 BIGA=A(I,J) L(K)=I M(K)=J 20 CONTINUE IF (ABS(BIGA).LT.(1.0E-7)) 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./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 150 CONTINUE RETURN END C ************************************ SUBROUTINE MULT(E , F , G , NROWE , NROWF , NCOLF) C C C THIS SUBROUTINE CALCULATES THE PRODUCT G OF TWO MATRICES C E AND F. C C IMPLICIT REAL*8(A-H,O-Z) REAL*8 E(NROWE,NROWF) , F(NROWF,NCOLF) , G(NROWE,NCOLF) DO 831 J = 1,NCOLF DO 832 I = 1,NROWE G(I,J) = 0.0 832 CONTINUE 831 CONTINUE DO 833 J = 1,NCOLF DO 834 I = 1,NROWE DO 835 K = 1,NROWF G(I,J) = E(I,K)*F(K,J) + G(I,J) 835 CONTINUE 834 CONTINUE 833 CONTINUE RETURN END 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 C ************************************************************* C C Calculation of S(beta) C C ************************************************************* Subroutine sbeta(time) implicit real*8(A-H,O-Z),integer(i-n) integer t,time common y(10,51000),bj(100),w1(51000),y1(51000),b1(51000), $ wdist(51000,10),agg(51000), $ q(51000),ya(51000),bet(15),sbet(15), $ cj(100),fru(51000),c1(51000), $ umar(51000),temp(51000),bold(100),bound,nper,np,num,i real*8 dep, ind(15), coef(15), xx(15,15), xy(15) c c a NLLS regression is done by running an OLS regression on the derivati- c ves of pred. 'ind' contains these derivatives, and 'dep' c has the dependent variables in the corresponding OLS regression. c c c num is the number of iterations. very often you can speed up c convergence by setting num = 1. The inaccuracy this creates c disappears in subsequent iterations on S(beta) c Do 20 iter=1,num cwrite(*,*) (sbet(kk),kk=1,np) do 10 i= 1,np xy(i) = 0. do 10 j=i,np 10 xx(i,j) = 0. Do 30 t= 150, time-1 ind(1) = exp(sbet(2)*ya(t) + & sbet(3)*w1(t) + $ sbet(4)*wdist(t,3) + $ sbet(5)*agg(t) + $ sbet(6)*w1(t)*y(1,t) + $ sbet(7)*wdist(t,1) + $ sbet(8)*wdist(t,2) + $ sbet(9)*wdist(t,4) + $ sbet(10)*(b1(t)/100.) + $ sbet(11)*(b1(t)/100.)**2. + $ sbet(12)*(b1(t)/100.)**3. + $ sbet(13)*ya(t)*agg(t) + $ sbet(14)*b1(t)*agg(t)/100.+ $ sbet(15)*y(1,t)*agg(t)/100.) prov = ind(1) * sbet(1) ind(2) = prov * ya(t) ind(3) = prov * w1(t) ind(4) = prov * wdist(t,3) ind(5) = prov * agg(t) ind(6) = prov * w1(t)*y(1,t) ind(7) = prov * wdist(t,1) ind(8) = prov * wdist(t,2) ind(9) = prov * wdist(t,4) ind(10) = prov * b1(t)/100. ind(11) = prov * (b1(t)/100.)**2. ind(12) = prov * (b1(t)/100.)**3. ind(13) = prov * ya(t)*agg(t) ind(14) = prov * b1(t)*agg(t)/100. ind(15) = prov * y(1,t)*agg(t)/100. dep = umar(t+1) do 33 j = 2,np dep = dep + ind(j)*sbet(j) 33 continue do 35 i= 1,np xy(i) = xy(i) + ind(i)*dep do 35 j=i,np 35 xx(i,j) = xx(i,j) + ind(i) * ind(j) 30 continue do 40 i=2,np do 40 j=1,i-1 40 xx(i,j) = xx(j,i) call invert(xx,np,d) call mult(xx,xy,coef,np,np,1) diff = 0. do 50 i=1,np diff = diff + abs( (sbet(i)-coef(i))/sbet(i) ) sbet(i) = coef(i) 50 continue if (diff.lt.(.000001)) go to 100 20 continue write(*,'(a)') ' NLLS not solved' write(*,'(a)') ' ---------------' 100 continue 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 func(qq) implicit real*8(a-h,o-z),integer(i-n) common y(10,51000),bj(100),w1(51000),y1(51000),b1(51000), $ wdist(51000,10),agg(51000), $ q(51000),ya(51000),bet(15),sbet(15), $ cj(100),fru(51000),c1(51000), $ umar(51000),temp(51000),bold(100),bound,nper,np,num,i do 220 j = 1,nper cj(j) = qq/temp(j) bj(j) = ( y(j,i) + bold(j) - cj(j) ) / qq if(bj(j).lt.bound) then bj(j) = bound endif 220 continue c c c check aggregate demand conditions c c btot = 0.0 do 235 j = 1,nper btot = btot + bj(j) 235 continue func = btot 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 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