Cc Open program Cc convergence criterion with series cc Setup of the model: cc in file psis.for the parametrized expectation; cc in der.for the derivatives w.r.t. parameters; cc in ginvert.for the procedure to calculate the state cc variables from the variables(t-1) and the parametrized cc expectations cc also in ginvert.for the expression be predicted in cc the par. expectation (phi) C Subroutines for generating random numbers, C and matrix algebra are included at the end. C It needs the input files: C openpar.dat - (contains parameter values) C openin.dat - (contains initial conditions for the iterations C on beta, and for the capital stock) C C names of arrays and parameters: c x and xlog - the vector of variables and of their logs as in MM c nx - is the total number of x variables c nxini - is the number of x variables that have to be initialized c (they have to be the first ones in the z vector) c phi - the function of x(t+1) that has to be predicted (as in MM) c psi - the param. expectation c thed - vector of AR(1) shocks (as u in MM) c (the shock part is still unchanged, so thed is a single skock c and ro and sig are single parameters) c bet(.) - the parameters in the param. expectation C sbet(.) - result of non-linear regression c neq - the number of par. exp. expressions c np(.) - the number of parameters in each param. exp. expression c alpha - a vector of parameters given in the model c nal - number of given parameters alpha C del - discount factor in utility C depr - undepreciated portion of capital C gam - relative risk aversion C ro- AR(1) parameter of log(thed(t)) C sig - standard deviation of innovation in log(thed(t)) C When solving the model for different parameters keep in mind that it C is safest to change parameters in small steps. See discussion of C homotopy in den Haan & Marcet, JBES. program pea implicit real*8(A-H,O-Z) integer t, time, nx, nxini, neq, np(10), nal, hom real z, rstart real xmean(15), xstd(15), xcorr(15,15) real b(10,10) real psiold(9000,10) character*100 alnotes common x(9000,12), xlog(9000,12), nx, nxini, $ thed(9000), thedlog(9000), phi(9000,10), psi(9000,10), $ bet(10,10), sbet(10,10), np, neq, $ alpha(20), nal open(8, file='openpar.dat' , status='old' ) read(8,*) time, nfor read(8,*) neq read(8,*) (np(j), j = 1,neq) read(8,*) nx, nxini read(8,*) (x(1,j), j = 1,nxini) read(8,*) eps, acc, accse, accnlr,iseed open(5, file='alpha.dat' , status='old' ) read(5,*) nal read(5,*) (alpha(j), j = 1,nal) read(5,*) hom if (hom.eq.1) then read(5,*) ialp, alpend, step endif read(5,*) ro, sig read(5,*) alnotes if (hom.ne.1) goto 24 alpzero=alpha(ialp) alvar = (alpend - alpzero) check = alvar*step if(check.lt.0) step = -step if (alvar.eq.0) hom=0 nshom = abs(alvar)/step if (nshom.gt.100) then write(*,*) 'gap between initial and final value for homotopy $ very large (',nshom,' steps), are you sure?(yes/no)' read(*,*) answer if (answer.ne."yes") stop endif 24 continue c c the parameters should be in different rows for different equations c in file openin.dat c open(2, file='openin.dat', status='old') do 26 j = 1, neq ii=np(j) read(2,*) (bet(i,j), i=1,ii) 26 continue do 27 i = 1, nxini xlog(1,i) = log(x(1,i)) 27 continue C creation of exogenous income process C z=rstart(iseed) call exproc(time, thed, thedlog, ro, sig) C C iterations on S to find fixed point. First loop sets the initial C condition for the first NLLS regression; after iter=1, the initial C condition is the previous sbet. C Do 10 j=1,neq ii=np(j) Do 10 i=1, ii sbet(i,j) = bet(i,j) 10 continue C start iterations C difse is the distance between the simul.series with bet and sbet C divi is the divisor C homotopy step start 12 continue divi = (time-150+1)*neq difse = 0. Do 30 iter = 1, 30 difseold = difse difse = 0. C C subroutine simul calculates solutions to the endogenous series for C a given bet, and the sum of the squares of (psi - psiold) C open(7, file='openout.dat' , status='unknown' ) call simul(time) C calculation of S(bet)=value of NLR when bet is used in generating C simulations C call sbeta(time,accnlr) C Calculate psi with the new sbeta on the old simulated series Do 211 j=1,neq ii=np(j) Do 211 i=1, ii b(i,j) = sbet(i,j) 211 continue Do 311 t = 150, time do 312 j=1, neq 312 psiold(t,j) = psi(t,j) include 'psis.for' Do 311 j = 1, neq difse = difse + % (psi(t,j) - psiold(t,j))*(psi(t,j) - psiold(t,j)) 311 continue C C Check convergence criterion C difse = difse / divi write(*,*) '¦¦b-Sb¦¦ = ' , difse write(7,*) '¦¦b-Sb¦¦ = ' , difse div = eps * max(0.00001,difse) didi = ((difseold - difse) / div) write(*,*) '-%diff= ' , didi write(7,*) '-%diff= ' , didi write(*,*) if (difse.lt.accse) then write(*,*) 'fixed point attained' goto 80 endif if (didi.ge.0. and. didi.lt.acc) then write(*,*) 'sum of squares decreasing too slowly' goto 80 endif C C updating of beta C Do 47 j=1, neq ii = np(j) Do 47 i=1, ii bet(i,j) = bet(i,j)*(1-eps) + sbet(i,j)*eps 47 continue write(*,*) write(7,*) iter do 50 j=1,neq ii=np(j) write(7,300) (bet(i,j), i=1,ii ) write(7,300) (sbet(i,j), i=1,ii ) 50 continue 30 continue write(7,'(a)') ' iterations on beta did not converge' write(7,'(a)') ' -----------------------------------' 80 continue C C the fixed point is written in the jbesin.dat file, to be used as initial C condition for the next computation. Note how the old 'jbesin.dat' file C will be lost unless it was saved before. C rewind 2 do 31, j=1, neq ii=np(j) write(2,*) ( bet(i,j), i=1,ii ) 31 continue C Homotopy step if (hom.ne.1) goto 85 w1 = (alpha(ialp)-alpzero)/alvar if (w1.ge.1) then write(*,*) 'homotopy concluded from...', alpzero write(*,*) 'to...', alpha(ialp) goto 85 endif alpha(ialp) = alpha(ialp) + step write(7,*) '>>>>>>>>>>>>>>>' write(7,'(a6,i2,a2,f10.5)') 'alpha(', ialp,') =', alpha(ialp) write(*,'(a6,i2,a2,f10.5)') 'alpha(', ialp,') =', alpha(ialp) goto 12 85 continue C C resulting simulations of series are written to an output file C with sample statistics write(7,*) 'simul. series ' do 118 i = 1, nx xmean(i) = 0 xstd(i) = 0 xcorr(i,i) = 1 do 119 j=1, i-1 119 xcorr(i,j) = 0 118 continue do 120 t=2, time do 128 i=1, nx 128 xmean(i) = xmean(i) + x(t,i) write(7,300) (x(t,i), i =1, nx) 120 continue do 131 i = 1, nx 131 xmean(i) = xmean(i)/(time-1) write(7,*) 'Average' write(7,300) (xmean(i), i=1, nx) do 150 t=2, time-1 do 150 i=1, nx xstd(i) = xstd(i) + (x(t,i)-xmean(i))*(x(t,i)-xmean(i)) do 159 j=1, i-1 159 xcorr(i,j) = xcorr(i,j) + £ (x(t,i) - xmean(i))*(x(t,j) - xmean(j)) 150 continue do 160 i=1, nx 160 xstd(i) = sqrt(xstd(i)/(time-1)) write(7,*) 'st. dev.' write(7,300) (xstd(i), i=1, nx) do 165 i=1, nx do 165 j=1, i-1 if (xstd(i)*xstd(j).gt.0) then xcorr(i,j) = xcorr(i,j)/((time-1)*xstd(i)*xstd(j)) xcorr(j,i) = xcorr(i,j) endif 165 continue write(7,*) 'corr. mat' do 170 i=1, nx write(7,300) (xcorr(i,j), j=1, nx) 170 continue 300 format(1x,6(f12.7,1x)) stop end C **************************************************************** C this subroutine simulates the series for a given beta C subroutine simul(time) implicit real*8(A-H,O-Z) integer t, time real rhs real b(10,10) common x(9000,12), xlog(9000,12), nx, nxini, $ thed(9000), thedlog(9000), phi(9000,10), psi(9000,10), $ bet(10,10), sbet(10,10), np(10), neq, $ alpha(20), nal Do 11 t = 2, time Do 21 j=1, neq ii=np(j) do 21 i=1, ii b(i,j)=bet(i,j) 21 continue include 'psis.for' include 'ginvert.for' xxx = (x(t,1)-x(t,2))/x(t,1) if (t.lt.90 .or. t.gt.3950) then write(7,*) t write(7,300) x(t,1), x(t,2), psi(t,1), psi(t,2) write(7,300) x(t,3), thed(t), xxx endif 300 format(1x,6(f12.8,1x)) 11 continue return end C ************************************************************* C Calculation of S(beta) C Subroutine sbeta(time,accnlr) implicit real*8(A-H,O-Z) integer t, time common x(9000,12), xlog(9000,12), nx, nxini, $ thed(9000), thedlog(9000), phi(9000,10), psi(9000,10), $ bet(10,10), sbet(10,10), np(10), neq, $ alpha(20), nal real*8 dep, der(10), coef(10), xx(10,10), xy(10) real*8 b(10,10) 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 See appendix of den Haan & Marcet for a detailed description of this C algorithm. C If accnlr is set to a large number, the subroutine only performs C one iteration in the NLLS algorithm; this makes the algorithm faster C without loosing any precision near the fixed point. C Do 20 je=1, neq ii=np(je) Do 21 iter=1,100 do 10 i= 1,ii xy(i) = 0. do 10 j=i,ii xx(i,j) = 0. 10 continue Do 30 t= 152, time-1 do 22 i=1, ii b(i,je) = sbet(i,je) 22 continue include 'psis.for' include 'der.for' dep = phi(t,je) - psi(t,je) do 33 i=1, ii 33 dep = dep + der(i) * sbet(i,je) do 35 i= 1, ii xy(i) = xy(i) + der(i) * dep do 35 j= i ,ii 35 xx(i,j) = xx(i,j) + der(i) * der(j) 30 continue do 40 i=2,ii do 40 j=1,i-1 40 xx(i,j) = xx(j,i) c do 174 i=1, ii c do 174 j=1,ii c write(*,*) 'elem di xx ', i, ',', j, ' ', xx(i,j) c 174 continue call invert(xx,ii,d) call mult(xx,xy,coef,ii) diffe = 0. do 50 i = 1,ii div = abs(sbet(i,je)) div = max(div, 0.0001) diffe = diffe + abs( (sbet(i,je)-coef(i))/div ) 50 sbet(i,je) = coef(i) if (diffe.lt.accnlr) go to 100 21 continue write(*,'(a)') ' NLLS not solved' write(*,'(a)') ' -----4-2---------' stop 100 continue 20 continue return end C *************************************************************** C this subroutine creates the exogenous income process; it throws away C the first 50 elements C subroutine exproc(time, thed, thedlog, ro, sig ) implicit real*8(A-H,O-Z) integer time, t real*8 thed(9000), thedlog(9000) real z, rnor thed(1) = 1 thedlog(1) = 0 do 10 t=2,time z = rnor() thedlog(t) = ro*thedlog(t-1) + dble(z)*sig 10 thed(t) = exp( thedlog(t) ) return end C ********************************************************** C 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) REAL*8 A(10,*) , L(10) , M(10) 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 , NROW) C C C THIS SUBROUTINE CALCULATES THE PRODUCT G OF A MATRIX C E AND A VECTOR F. C C IMPLICIT REAL*8(A-H,O-Z) REAL*8 E(10,*) , F(10) , G(10) DO 832 I = 1,NROW G(I) = 0.0 832 CONTINUE DO 834 I = 1,NROW DO 835 K = 1,NROW G(I) = E(I,K)*F(K) + G(I) 835 CONTINUE 834 CONTINUE RETURN END C C ***************************************** C ***************************************** C REAL FUNCTION RNOR() REAL AA,B,C,C1,C2,PC,X,Y,XN,V(65),RSTART,U(17),S,T,UN INTEGER J,IA,IB,IC,II,JJ,ID,III,JJJ SAVE U,II,JJ C DATA AA,B,C/12.37586,.4878992,12.67706/ DATA C1,C2,PC,XN/.9689279,1.301198,.1958303E-1,2.776994/ DATA V/ .3409450, .4573146, .5397793, .6062427, .6631691 +, .7136975, .7596125, .8020356, .8417227, .8792102, .9148948 +, .9490791, .9820005, 1.0138492, 1.0447810, 1.0749254, 1.1043917 +,1.1332738, 1.1616530, 1.1896010, 1.2171815, 1.2444516, 1.2714635 +,1.2982650, 1.3249008, 1.3514125, 1.3778399, 1.4042211, 1.4305929 +,1.4569915, 1.4834526, 1.5100121, 1.5367061, 1.5635712, 1.5906454 +,1.6179680, 1.6455802, 1.6735255, 1.7018503, 1.7306045, 1.7598422 +,1.7896223, 1.8200099, 1.8510770, 1.8829044, 1.9155830, 1.9492166 +,1.9839239, 2.0198430, 2.0571356, 2.0959930, 2.1366450, 2.1793713 +,2.2245175, 2.2725185, 2.3239338, 2.3795007, 2.4402218, 2.5075117 +,2.5834658, 2.6713916, 2.7769943, 2.7769943, 2.7769943, 2.7769943/ C Load data array in case user forgets to initialize. C This array is the result of calling UNI 100000 times C with seed 305. DATA U/ *0.8668672834288, 0.3697986366357, 0.8008968294805, *0.4173889774680, 0.8254561579836, 0.9640965269077, *0.4508667414265, 0.6451309529668, 0.1645456024730, *0.2787901807898, 0.06761531340295, 0.9663226330820, *0.01963343943798, 0.02947398211399, 0.1636231515294, *0.3976343250467, 0.2631008574685/ C DATA II,JJ / 17, 5 / C C***FIRST EXECUTABLE STATEMENT RNOR C C Fast part... C C C Basic generator is Fibonacci C UN = U(II)-U(JJ) IF(UN.LT.0.0) UN = UN+1. U(II) = UN C U(II) and UN are uniform on [0,1) C VNI is uniform on [-1,1) VNI = UN + UN -1. II = II-1 IF(II.EQ.0)II = 17 JJ = JJ-1 IF(JJ.EQ.0)JJ = 17 C INT(UN(II)*128) in range [0,127], J is in range [1,64] J = MOD(INT(U(II)*128),64)+1 C Pick sign as VNI is positive or negative RNOR = VNI*V(J+1) IF(ABS(RNOR).LE.V(J))RETURN C C Slow part; AA is a*f(0) X = (ABS(RNOR)-V(J))/(V(J+1)-V(J)) C Y is uniform on [0,1) Y = U(II)-U(JJ) IF(Y.LT.0.0) Y = Y+1. U(II) = Y II = II-1 IF(II.EQ.0)II = 17 JJ = JJ-1 IF(JJ.EQ.0)JJ = 17 C S = X+Y IF(S.GT.C2)GO TO 11 IF(S.LE.C1)RETURN IF(Y.GT.C-AA*EXP(-.5*(B-B*X)**2))GO TO 11 IF(EXP(-.5*V(J+1)**2)+Y*PC/V(J+1).LE.EXP(-.5*RNOR**2))RETURN C C Tail part; .3601016 is 1./XN C Y is uniform on [0,1) 22 Y = U(II)-U(JJ) IF(Y.LE.0.0) Y = Y+1. U(II) = Y II = II-1 IF(II.EQ.0)II = 17 JJ = JJ-1 IF(JJ.EQ.0)JJ = 17 C X = 0.3601016*LOG(Y) C Y is uniform on [0,1) Y = U(II)-U(JJ) IF(Y.LE.0.0) Y = Y+1. U(II) = Y II = II-1 IF(II.EQ.0)II = 17 JJ = JJ-1 IF(JJ.EQ.0)JJ = 17 IF( -2.*LOG(Y).LE.X**2 )GO TO 22 RNOR = SIGN(XN-X,RNOR) RETURN 11 RNOR = SIGN(B-B*X,RNOR) RETURN C C C Fill ENTRY RSTART(ISEED) IF(ISEED.NE.0) THEN C C Set up ... C Generate random bit pattern in array based on given seed C II = 17 JJ = 5 IA = MOD(ABS(ISEED),32707) IB = 1111 IC = 1947 DO 2 III = 1,17 S = 0.0 T = .50 C Do for each of the bits of mantissa of word C Loop over 64 bits, enough for all known machines C in single precision DO 3 JJJ = 1,64 ID = IC-IA IF(ID.GE.0)GOTO 4 ID = ID+32707 S = S+T 4 IA = IB IB = IC IC = ID 3 T = .5*T 2 U(III) = S ENDIF C Return floating echo of ISEED RSTART=ISEED RETURN END C C **********************************************************