c This program solves dynamic economic models c by a parametrized expectations algorithm, c as described in Marcet and Lorenzoni (1998). c c Variant for EXAMPLE 3, with additional subroutines c to solve numerically for labor supply. c c It needs the program files: c - ginvert.for (the steps to invert the system g(...)=0) c - psis.for (the parametrized expressions psi for the expectation c terms) c - der.for (derivatives w.r.t. parameters of the psi) c c c It needs the input files: c - par.dat (contains parameter values controlling the algorithm) c - alpha.dat (contains parameter values for the model) c - pini.dat (contains initial conditions for the iterations c on beta) c c It creates the output files: c - openout.dat (containing information on algorithm steps c and summary statistics on the simulated series) c - series.dat (contains simulated series) c c c Names of arrays and parameters: c c z(t,i) and zlog(t,i) - Variable i and of its natural logarithm c c nz - is the total number of z variables c nzini - is the number of z variables that have to be initialized c (they have to be the first ones in the z vector) c (they may exceed the number of state variables in c the model for computational purposes) c thetalog(t,j) - normal AR(1) vector of shocks c theta(t,j) - exp(thetalog) c nthe - number of shock variables c c phi(t,i)- the expression inside the expectation that has c to be predicted c psi(t,i) - the parametrized expression used to c approximate the terms in expectation c bet(k) and gbet(k) - the parameters in the param. expectation c neq - the number of psi expressions c np(j) - the number of parameters in the j-th psi expression c c alpha(k) - a vector of parameters given in the model c nal - number of given parameters alpha c am(j,k), cm(j,k), dm(j) - arrays containing c parameters for the vector AR(1) process thetalog c c time - number of periods in simulation c nfor - maximum lead of variables in the par.expectation c eps - parameter controlling the adjustment of beta c along the algorithm iterations c acc - acceptance level for the convergence of the algorithm c accnlr - acc. level for the convergence of the nonlinear regression c iseed - integer for the random numer subroutines c maxita - maximum number of iterations for the algorithm c maxitr - maximum number of iterations for the NL regression c c Variables to control the homotopy steps: c hom - integer variable controlling the homotopy steps: c 0= no homotopy c 1= change an element (ialp) of vector alpha c 2= change an element (jam,kam) of matrix A c 3= change an element (jam,kam) of matrix C c alpend - target value for the parameter to be changed c step - step length c c c Variable: Maximum value: c time 10000 c nz 15 c nzini nz c nal 20 c neq 10 c np 10 c program pea implicit real*8(A-H,O-Z) integer t, time, nz, nzini, nthe, neq, & np(10), nal, nfor, hom, maxita, maxitr, & jam, kam real q, rstart real*8 zmean(15), zstd(15), zcorr(15,15) real*8 am(10,10), cm(10,10) character*3 answer character*100 alnotes common z(10000,15), zlog(10000,15), nz, nzini, $ theta(10000,10), thetalog(10000,10), nthe, $ phi(10000,10), psi(10000,10), $ bet(10,10), gbet(10,10), np, neq, $ alpha(20), nal, nfor CCCCCCCCCCCCCCCCCCCCC C C Additional line to use functions f and derf C external f, derf C CCCCCCCCCCCCCCCCCCCCC open(8, file='par.dat' , status='old' ) read(8,*) time, nfor read(8,*) neq read(8,*) (np(j), j = 1,neq) read(8,*) nz, nzini read(8,*) ( z(1,j), j = 1,nzini) read(8,*) eps, acc, accnlr, iseed read(8,*) maxita, maxitr open(5, file='alpha.dat' , status='old' ) read(5,*) nal read(5,*) (alpha(j), j = 1,nal) read(5,*) nthe do 20 j=1, nthe 20 read(5,*) (am(j,k), k=1, nthe) do 25 j=1, nthe 25 read(5,*) (cm(j,k), k=1, nthe) read(5,*) hom if (hom.eq.1) then read(5,*) ialp, alpend, step else if (hom.eq.2.or. hom.eq.3) then read(5,*) jam, kam, alpend, step endif read(5,*) alnotes C C Compute the number of steps for the homotopy and C check if it is too large. C if (hom.eq.1) then alpzero = alpha(ialp) else if (hom.eq.2) then alpzero = am(jam,kam) else if (hom.eq.3) then alpzero = cm(jam,kam) else goto 24 endif 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(*,*) 'The gap between starting and ending value for $ the homotopy is very large (',nshom,' steps), $ continue? (yes/no)' read(*,*) answer if (answer.ne."yes") stop endif 24 continue c Read initial parameters c The parameters are in different rows for different equations c in the file pini.dat. open(2, file='pini.dat', status='old') do 26 j = 1, neq ii=np(j) read(2,*) (bet(i,j), i=1,ii) 26 continue C Initial logs do 28 i = 1, nzini if ( z(1,i).le.0) goto 28 zlog(1,i) = log( z(1,i)) 28 continue C C Creation of exogenous markov process C q=rstart(iseed) call exproc(time, theta, thetalog, nthe, am, cm) C C Iterations on G 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 Gbet. C Do 10 j=1,neq ii=np(j) Do 10 i=1, ii gbet(i,j) = bet(i,j) 10 continue 12 continue C C Here the algorithm restarts, along the homotopy C Do 30 iter = 1, maxita C C Subroutine simul calculates solutions of the endogenous series for C a given bet C open(7, file='openout.dat' , status='unknown' ) call simul(time) C C Calculation of G(bet)=coefficients of NLR C with bet used in generating simulations C call gbeta(time,accnlr,maxitr) C C Updating of beta and checking convergence criterion C write(*,*) write(7,*) 'iter.', iter diff = 0. do 50 j=1,neq ii=np(j) write(7,300) (bet(i,j), i=1,ii ) write(7,300) (gbet(i,j), i=1,ii ) do 55 i=1,ii div=max(0.0001,gbet(i,j)) diff = diff + abs( ( bet(i,j)-gbet(i,j) ) / div ) bet(i,j) = bet(i,j)*(1-eps) + gbet(i,j)*eps 55 continue 50 continue write(*,'(a7,f10.5)') 'diff= ' , diff if ( diff.lt.acc ) go to 80 30 continue write(7,'(a)') ' iterations on beta did not converge' write(7,'(a)') ' -----------------------------------' 80 continue C Homotopy step if (hom.eq.1) then w1 = (alpha(ialp)-alpzero)/alvar else if (hom.eq.2) then w1 = (am(jam,kam)-alpzero)/alvar else if (hom.eq.3) then w1 = (cm(jam,kam)-alpzero)/alvar else goto 85 endif if (w1.ge.1) then write(*,*) 'homotopy concluded from...', alpzero write(*,*) 'to...', alpend goto 85 endif if (hom.eq.1) then alpha(ialp) = alpha(ialp) + step else if (hom.eq.2) then am(jam,kam) = am(jam,kam) + step else if (hom.eq.3) then cm(jam,kam) = cm(jam,kam) + step endif write(7,*) '>>>>>>>>>>>>>>>' if (hom.eq.1) then write(7,'(a6,i2,a2,f10.5)') 'alpha(', ialp,') =', alpha(ialp) write(*,'(a6,i2,a2,f10.5)') 'alpha(', ialp,') =', alpha(ialp) else if (hom.eq.2) then write(7,'(a6,i2,a2,i2,a2,f10.5)') 'am(', jam,',',kam,') =' $ , am(jam,kam) write(*,'(a6,i2,a2,i2,a2,f10.5)') 'am(', jam,',',kam,') =' $ , am(jam,kam) else if (hom.eq.3) then write(7,'(a6,i2,a2,i2,a2,f10.5)') 'cm(', jam,',',kam,') =' $ , cm(jam,kam) write(*,'(a6,i2,a2,i2,a2,f10.5)') 'cm(', jam,',',kam,') =' $ , cm(jam,kam) endif goto 12 85 continue C C The fixed point is written in the pini.dat file, to be used as initial C condition for the next computation. Note how the old 'pini.dat' file C will be lost unless it was saved before. C rewind 2 do 35, j=1, neq ii=np(j) write(2,*) ( bet(i,j), i=1,ii ) 35 continue C C File alpha.dat is updated C if (hom.eq.0) goto 40 rewind 5 write(5,'(i5)') nal write(5,300) (alpha(j), j = 1,nal) write(5,'(i5)') nthe do 320 j=1, nthe 320 write(5,300) (am(j,k), k=1, nthe) do 325 j=1, nthe 325 write(5,300) (cm(j,k), k=1, nthe) write(5,'(i5)') hom if (hom.eq.1) then write(5,'(i5,f12.8,f12.8)') ialp, alpend, step else if (hom.eq.2.or. hom.eq.3) then write(5,'(i5,i5,f12.8,f12.8)') jam, kam, alpend, step endif write(5,45) alnotes 40 continue 45 format('''',a70,'''') C C Resulting simulations of series are written C to the output file series.dat C open(4, file='series.dat', status='unknown') C C Summary statistics are computed and written to the output file C openout.dat C call summar(z,zmean,zstd,zcorr,nz,2,time-1) do 120 t=2, time-1 write(4,300) ( z(t,j), j=1,nz), (theta(t,k), k=1, nthe) 120 continue write(7,*) 'Mean' write(7,300) (zmean(i), i=1, nz) write(7,*) 'Standard deviations' write(7,300) (zstd(i), i=1, nz) write(7,*) 'Correlation matrix' do 170 i=1, nz write(7,300) (zcorr(i,j), j=1, nz) 170 continue 300 format(1x,6(f12.8,1x)) stop end C***************************** C Subroutines C C***************************** C C Subroutine to simulate the series for a given beta C subroutine simul(time) implicit real*8(A-H,O-Z) integer t, time real b(10,10) common z(10000,15), zlog(10000,15), nz, nzini, $ theta(10000,10), thetalog(10000,10), nthe, $ phi(10000,10), psi(10000,10), $ bet(10,10), gbet(10,10), np(10), neq, $ alpha(20), nal, nfor CCCCCCCCCCCCCCCCCCCCCCCCCCC C C Additional lines to use functions f and derf C external f, derf whold = z(1,2) C CCCCCCCCCCCCCCCCCCCCCCCCCC 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' c if (t.lt.90) then c write(7,*) t c write(7,300) ( z(t,i), i=1, nz) c write(7,300) (psi(t,j), j=1, neq) c write(7,300) (theta(t,j), j=1, nthe) c endif 300 format(1x,6(f12.8,1x)) 11 continue return end C ************************************************************* C Subroutine for the calculation of G(beta) (NLR) C Subroutine gbeta(time,accnlr,maxitr) implicit real*8(A-H,O-Z) integer t, time, maxitr common z(10000,15), zlog(10000,15), nz, nzini, $ theta(10000,10), thetalog(10000,10), nthe, $ phi(10000,10), psi(10000,10), $ bet(10,10), gbet(10,10), np(10), neq, $ alpha(20), nal, nfor 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 the psi. der.for contains these derivatives. 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, maxitr do 10 i= 1,ii xy(i) = 0. do 10 j=i,ii xx(i,j) = 0. 10 continue Do 30 t= 152, time-nfor do 22 i=1, ii b(i,je) = gbet(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) * gbet(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) call lsol(xx,xy,coef,ii,je) diff = 0. do 50 i = 1,ii div = abs(gbet(i,je)) div = max(div, 0.0001) diff = diff + abs( (gbet(i,je)-coef(i))/div ) 50 gbet(i,je) = coef(i) if (diff.lt.accnlr) go to 100 21 continue write(*,'(a)') ' NLLS not solved' write(*,'(a)') ' ----------------' stop 100 continue 20 continue return end C *************************************************************** C This subroutine creates the exogenous markov process; C it throws away the first 50 elements C subroutine exproc(time, theta, thetalog, nthe, am, cm ) implicit real*8(A-H,O-Z) integer time, t real*8 theta(10000,10), thetalog(10000,10) real*8 am(10,*), cm(10,*) real q(10), rnor do 5 j=1, nthe theta(1,j) = 1 do 5 t=1, time 5 thetalog(1,j) = 0 do 10 t=2, time do 15 j=1, nthe 15 q(j) = rnor() do 30 j=1, nthe do 25 k=1, nthe 25 thetalog(t,j) = am(j,k) * thetalog(t-1,k) + thetalog(t,j) do 20 k=1, nthe 20 thetalog(t,j) = cm(j,k) * dble(q(k)) + thetalog(t,j) 30 theta(t,j) = exp( thetalog(t,j) ) 10 continue return end C C************************ C This subroutines computes summary statistics of the z C for the periods ti to tii C subroutine summar(x,xmean,xstd,xcorr,nx,ti,tii) implicit real*8(A-H,O-Z) integer ti, tii, t real*8 x(10000,*), xmean(*), xstd(*), & xcorr(15,*), xcov(15,15) div = tii - ti + 1 do 118 i = 1, nx 118 xmean(i) = 0 do 120 t= ti, tii do 128 i=1, nx 128 xmean(i) = xmean(i) + x(t,i) 120 continue do 131 i = 1, nx 131 xmean(i) = xmean(i)/div do 140 i=1, nx do 140 j=1, i 140 xcov(i,j) = 0 do 150 t = ti, tii do 150 i=1, nx do 159 j=1, i 159 xcov(i,j) = xcov(i,j) + & (x(t,i) - xmean(i))*(x(t,j) - xmean(j)) 150 continue do 155 i=1, nx do 155 j=1, i 155 xcov(i,j) = xcov(i,j)/div do 160 i=1, nx 160 xstd(i) = sqrt(xcov(i,i)) do 165 i=1, nx do 165 j=1, i if (xstd(i)*xstd(j).eq.0) then xcorr(i,j) = 0 goto 164 endif xcorr(i,j) = xcov(i,j)/(xstd(i)*xstd(j)) 164 xcorr(j,i) = xcorr(i,j) 165 continue return end C C ********************************************************** C SUBROUTINE LSOL(A,B,Y,N,npe) C This subroutine solves the linear system Ay=b C where A is a symmetric positive definite matrix, C by means of Cholesky factorization. C npe is the index of the psi for which we are performing C the nonlinear regression. This subroutines keeps track C of it, to signal this information in case of troubles C solving the linear system (i.e. inverting XX) IMPLICIT REAL*8(A-H,O-Z) real*8 a(10,*), p(10), y(*), b(*) integer n, npe C Nondiagonal part of triangular matrix T is stored in the C lower triangular part of A, the elements on the diagonal C of T are stored in p. T is such that TT'=A do 20 j=1, n do 20 i=j, n sum = a(j,i) do 30 l=1, j-1 sum = sum - a(j,l) * a(i,l) 30 continue if (i.eq.j) then if (sum.lt.0) then write(*,*) 'Failure inverting XX in LS' stop endif p(j) = sqrt(sum) else a(i,j) = sum/p(j) endif 20 continue C TT'x=b is solved in two steps: C 1. T'x is obtained solving Tq=b C 2. x is obtained solving T'x=q C (q is stored in x) do 50 i=1, n sum = b(i) do 60 j=1, i-1, 1 sum = sum - y(j)*a(i,j) 60 continue if (p(i).eq.0) then write(*,*) 'cannot invert XX in LS ', npe stop endif y(i) = sum/p(i) 50 continue do 70 i=n, 1, -1 sum = y(i) do 80 j=n, i+1, -1 sum = sum - y(j)*a(j,i) 80 continue y(i) = sum/p(i) 70 continue return end 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 cccccccccccccc c c Additional Subroutine to solve numerically for labor supply c cccccccccccccc subroutine fderzero(f,derf,u1,u2,wh,prec,isol,alpha,fac2) implicit real*8 (A-H,O-Z) real*8 alpha(20) external f, derf u1=.0 u2=1. prec = .01 c isol remains zero if convergence is not reached isol = 0 do 10 it=1, 100 fval = f(wh,alpha,fac2) wo = wh - alpha(9)*fval/derf(wh,alpha,fac2) if (wo.lt.u1) then wh = (u1+wh)/2 else if (wo.gt.u2) then wh = (wh+u2)/2 else wh = wo endif if( abs(fval) .lt. prec ) then isol = it go to 20 endif 10 continue 20 continue return end ccccccccccccccccccc c c Additional functions for the numerical subroutine c for labor supply c ccccccccccccccccccc double precision function f(y,alpha,fac2) implicit real*8(A-H,O-Z) real*8 alpha(20) f = ((1-y)**alpha(4))-fac2*(y**(-alpha(3))) end double precision function derf(y,alpha,fac2) implicit real*8(A-H,O-Z) real*8 alpha(20) derf = - alpha(4)*((1-y)**(alpha(4)-1)) + & alpha(3)*fac2*(y**(-alpha(3)-1)) end