c c c param1.f, program to solve the standard neo-classical c growth model using parameterized expectations c c version for 211B winter 1993 c c C this program calculates solution to the simple growth model C ( first-order) C C december 16 1989 C C C It needs the input files: C aipar.dat - (contains parameter values) C ain.dat - (contains initial conditions for the iterations C on beta, and for the capital stock C c an example for the two files is the following: C **************************************************************** c aipar.dat: c 2500 c .95 .001 .95 -1. .33 .95 c 1. 0.0001 -50 c time c ro sig del gam alpha depr c eps acc idum c idum is the initial seed for the random number generator c it has to be an integer LESS than zero c ain.dat: c 1.92466704 -.57492310 -.51443190 c 5.62297797555693 C **************************************************************** implicit real*8(A-H,O-Z) integer t, time common c(5000),xk(5000),umar(5000),xklog(5000),thedlog(5000), $ thed(5000), del, gam, eps, bet(10), sbet(10), depr,alpha, np c c c np is the number of parameters in the polynomial c c np = 3 c c c reading the input data c c open(8, file='aipar.dat' , status='old' ) open(2, file='ain.dat', status='old') read(8,*) time read(8,*) ro, sig, del, gam, alpha, depr, eps, acc, idum read(2,*) ( bet(i),i=1,np ), xk(1) xklog(1) = log(xk(1) ) write(*,*) time write(*,*) (bet(i),i=1,np) C C C creation of exogenous productivity shock C C call exproc(time, thed, thedlog, ro, sig,idum) c c c initializing of sbeta c c Do 10 i=1,np 10 sbet(i) = bet(i) C the "do 30" is the important loop in this program. It does the C iterations on S to find the 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. Do 30 iter = 1, 300 C C C I: subroutine simul calculates solutions to the endogenous C series for a given beta C call simul(time) C C II: calculation of S(beta) C call sbeta(time) C C C III: updating of beta and convergence criterion C C C diff is the difference between beta and S(beta); C eps is the parameter that determines the weight of C the old beta and S(beta) in determining the new beta C C for this model eps can easily be equal to 1, for more complicated C models eps has to be closer to 0. C write(*,*) write(*,*) iter write(*,300) (bet(i), i=1,np ) diff = 0. do 50 i=1,np diff = diff + abs( ( bet(i)-sbet(i) ) / bet(i) ) 50 bet(i) = bet(i)*(1-eps) + sbet(i)*eps write(*,300) (sbet(i), i=1,np ) write(*,'(a7,f10.5)') 'diff= ' , diff if ( diff.lt.acc ) $ go to 80 30 continue write(*,'(a)') ' iterations on beta did not converge' write(*,'(a)') ' -----------------------------------' 80 continue C C the fixed point is written in the ain.new file, to be used as initial C condition for the next computation. C open( 3, file = 'ain.new', status= 'unknown' ) write(3,300) (bet(i), i=1,np ) C C C calculation of some properties of the solution: means and C variances of the consumptions C C cmean = 0. do 90 t = 50, time-1 90 cmean = cmean + c(t) cmean = cmean / dble( time-50 ) cvar = 0. do 100 t = 50, time-1 100 cvar = cvar + ( c(t) - cmean )*( c(t) - cmean ) cvar = cvar / dble( time-50 ) xkmean = 0. do 110 t = 50, time-1 110 xkmean = xkmean + xk(t) xkvar = 0.0 xkmean = xkmean / dble( time-50 ) do 120 t = 50, time -1 120 xkvar = xkvar + (xk(t)-xkmean)**2. xkvar = xkvar / dble( time-50 ) 20 continue write(*,*) cmean,xkmean write(*,*) cvar,xkvar 300 format(1x,6(f12.8,1x)) write(3,*) xkmean stop end C **************************************************************** C C this subroutine simulates the series for a given beta C C **************************************************************** subroutine simul(time) implicit real*8(A-H,O-Z) integer t, time common c(5000),xk(5000),umar(5000),xklog(5000),thedlog(5000), $ thed(5000), del, gam, eps, bet(10), sbet(10), depr,alpha, np gaminv = 1/gam Do 10 t = 2,time C C 'rhs' stands for the right hand side in the Euler equation. C NOte that rhs is also the marginal utility of consumption C rhs = del * bet(1) * exp( xklog(t-1)*bet(2) + $ thedlog(t)*bet(3) ) xkalp = xk(t-1)**alpha c(t) = rhs**gaminv prod = thed(t) * xkalp + depr*xk(t-1) if ( c(t).gt.prod ) then write(*,*) 'consumption too large' stop endif xk(t) = prod - c(t) xklog(t) = log( xk(t) ) C C 'umar' contains the stuff inside the conditional expectation that C we will try to predict in sbeta C umar(t) = rhs * ( thed(t) * alpha * xkalp / xk(t-1) $ + depr ) C write(*,*) c(t) , xk(t), xklog(t), thed(t), thedlog(t) 10 continue 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 c(5000),xk(5000),umar(5000),xklog(5000),thedlog(5000), $ thed(5000), del, gam, eps, bet(10), sbet(10), depr,alpha, np real*8 dep, ind(3), coef(3), xx(3,3), xy(3) 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 the dimension of ind, coef, xx and xy has to be set consistent c with np 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 num = 20 Do 20 iter=1,num 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( xklog(t-1)*sbet(2) + $ thedlog(t)*sbet(3) ) prov = ind(1) * sbet(1) ind(2) = prov * xklog(t-1) ind(3) = prov * thedlog(t) dep = umar(t+1) + ind(2)*sbet(2) + ind(3)*sbet(3) 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) ) 50 sbet(i) = coef(i) if (diff.lt.(.000001)) go to 100 20 continue write(*,'(a)') ' NLLS not solved' write(*,'(a)') ' ---------------' stop 100 continue end C C *************************************************************** C C this subroutine creates the exogenous income process C C *************************************************************** subroutine exproc(time, thed, thedlog, ro, sig, idum) implicit real*8(A-H,O-Z) integer time,t dimension thed(5000), thedlog(5000) thed(1) = 0. do 10 t=2,time z = gasdev(idum) thedlog(t) = ro*thedlog(t-1) + z*sig 10 thed(t) = exp( thedlog(t) ) 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(N,N) , 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 , 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