C C PEAPROI1.F, Fortran program to solve the neoclassical growth model C C - Projection method & Iterative solution procedure C C NOVEMBER 9 1998 C C Written by Christian Haefke, University of California, San Diego C Support by Wouter den Haan's NSF grant is gratefully acknowledged. C C Documentation of the program is given at the end of the program C C ***************************************************************** PROGRAM PEAPROI1.F IMPLICIT NONE C Meta Parameters C ---------------- INTEGER pok,pot,knod,tnod,simtime,actime PARAMETER (pok=3,pot=3) PARAMETER (knod=7,tnod=7) PARAMETER (simtime=10000,actime=50) c c pok: highest order of (log) capital stock c pot: highest order of (log) technology stock c knod: number of capital grid points c tnod: number of productivity grid points c simtime: number of simulation periods if simulation is performed c actime: number of points at which to perform an accuracy test c INTEGER dimgrid,psisize,sim PARAMETER (dimgrid=knod*tnod,psisize=(pok+1)*(pot+1)) PARAMETER (sim=max(simtime,actime)) C Model Parameters C ---------------- DOUBLE PRECISION alpha,delta,dfactor,nu,tmin,tmax,kmin,kmax,sigma DOUBLE PRECISION rho,psitol,lrate C For the solution step C --------------------- INTEGER r,psiread INTEGER psiter,maxiter,ki,ti DOUBLE PRECISION kg(knod),tg(tnod),kgrid(dimgrid),tgrid(dimgrid) DOUBLE PRECISION tsig,lnk1,lnt1,lnk2,yvec(dimgrid),yvecr DOUBLE PRECISION xmat(dimgrid,psisize),xxmat(psisize,psisize) DOUBLE PRECISION xr(psisize,dimgrid),xmatt(psisize,dimgrid) DOUBLE PRECISION xk(dimgrid,psisize),xt(dimgrid,psisize) DOUBLE PRECISION delpsi,norm DOUBLE PRECISION psi0(psisize),psi(psisize),psin(psisize) DOUBLE PRECISION savfun,confun,numi,scalup C For the simulation step C ----------------------- INTEGER simyes,idum DOUBLE PRECISION ks(sim),ts(sim+1),epsi(sim) DOUBLE PRECISION kmins,kmaxs,kmeans,tmins,tmaxs,gasdev C For the accuracy step C --------------------- INTEGER accur INTEGER tobs,kobs,obs DOUBLE PRECISION lnk1v,lnt1v,lnk2v,incr DOUBLE PRECISION appc,actc,tgfix,kgfix,lnt1fix,lnk1fix DOUBLE PRECISION diffcr,diffck,diffct C The Common Blocks C ----------------- COMMON /consfcn/ nu,tmin,tmax,kmin,kmax COMMON /savfcn/ alpha,delta,dfactor COMMON /aexpfcn/ rho,sigma C ****************************************************************** C C 1. Initialize the parameters C ============================ C 1.1 Initialize model parameters C ------------------------------- alpha = 0.33D0 dfactor = 0.95D0 nu = 2.0D0 delta = .95D0 sigma = .01D0 rho = 0.95D0 C 1.2 Initialize algorithm parameters and algorithm options C --------------------------------------------------------- kmin = 0.8D0 kmax = 1.2D0 tmin = -3.5D0 tmax = 3.5D0 maxiter = 1000 psitol = 0.0001D0 lrate = .4D0 idum = -500 simyes = 1 accur = 1 C 1.3 Compute the steady state C ---------------------------- ks(1) = (1.0D0-dfactor+delta*dfactor)/(alpha*dfactor) ks(1) = ks(1)**(1.0D0/(alpha-1.0D0)) kmin = log(ks(1)*kmin) kmax = log(ks(1)*kmax) tsig = (1.0D0-rho*rho)**.5D0 tsig = sigma/tsig tmin = tmin*tsig tmax = tmax*tsig C 2 Some Preparatory Work C ======================= C 2.1 Create the grid C ------------------- call chebnode(knod,kg) call chebnode(tnod,tg) DO 10 ki = 1,knod DO 11 ti = 1,tnod kgrid((ki-1)*tnod+ti) = kg(ki) tgrid((ki-1)*tnod+ti) = tg(ti) 11 CONTINUE 10 CONTINUE C 2.2 Create the Projection Matrix C -------------------------------- C Prepare the Polynomial CALL chebpol(pok,kgrid,dimgrid,xk) CALL chebpol(pot,tgrid,dimgrid,xt) DO 12 ki = 0,pok DO 13 ti = 0,pot DO 14 r = 1,dimgrid xmat(r,ki*(pot+1)+ti+1)=xk(r,ki+1)*xt(r,ti+1) 14 CONTINUE 13 CONTINUE 12 CONTINUE call mtrans(xmat,xmatt,dimgrid,psisize) call mmult(xmatt,xmat,xxmat,psisize,dimgrid,psisize) C In this case xxmat is diagonal, hence we use an inversion routine C for diagonal matrices. Otherwise use minv! call mdinv(xxmat,psisize,xxmat) call mmult(xxmat,xmatt,xr,psisize,psisize,dimgrid) C 2.3 Initialize the Coefficients of the Polynomial C ------------------------------------------------- C 2.3.1 Read the initial point from the specified file C and initialize psi C ---------------------------------------------------- OPEN (UNIT = 1,FILE ='psi0.in', STATUS = 'OLD') DO 20 r = 1,psisize READ (1, *, END = 21) psi0(r) psi(r) = psi0(r) 20 CONTINUE 21 psiread = r - 1 CLOSE(unit=1) C 2.3.2 Check the size C -------------------- IF ( psiread.NE.psisize ) THEN write(*,*) 'Wrong number of values for psi0, we need ', + psisize, ' not ',psiread stop ENDIF C ******************************************************************* C C 3 THE MAIN LOOP C =============== DO 100 psiter = 1,maxiter C 3.1 Compute the expected value (Step II) C ---------------------------------------- DO 110 r = 1,dimgrid lnk1 = scalup(kgrid(r),kmin,kmax) lnt1 = scalup(tgrid(r),tmin,tmax) lnk2 = savfun(lnk1,lnt1,psi,psisize,pok,pot) yvecr = numi(lnk2,lnt1,psi,psisize,pok,pot) yvec(r) = log(yvecr) 110 CONTINUE C 3.2 Update the polynomial coefficients (Steps III and IV) C --------------------------------------------------------- call mvmult(xr,yvec,psin,psisize,dimgrid) delpsi = norm(psi,psin,psisize) print *,'Iteration',psiter,' Change',delpsi IF (delpsi .LE. psitol) THEN goto 200 ELSE DO 120 r = 1,psisize psi(r) = lrate * psin(r) + + (1.0D0-lrate)*psi(r) 120 CONTINUE ENDIF C This ends the main loop. 100 CONTINUE C Save the Coefficients of the Polynomial C --------------------------------------- 200 OPEN (UNIT=3,FILE='c:\haefke\fortran\psi0.out',STATUS='unknown') DO 210 r=1,psisize write(3,*) psi(r) 210 CONTINUE CLOSE(UNIT = 3) C 3.3 Print Parameters C -------------------- print *,'' print *,'NEW VALUES FOR PSI' print *,psi print *,'OLD VALUES FOR PSI' print *,psi0 C ******************************************************************* C 4. Simulation C ============= IF (simyes .EQ. 1) THEN C 4.0 Start at the Steady state C ----------------------------- C ks(1) is steady state capital ks(1) = dlog(ks(1)) ts(1) = 0.0D0 C 4.1 Run the simulation C ---------------------- kmins = 1000000.0D0 tmins = 1000000.0D0 kmaxs = -1000000.0D0 tmaxs = -1000000.0D0 kmeans = 0.0D0 DO 300 r = 1,simtime epsi(r) = gasdev(idum) * sigma ts(r+1) = rho*ts(r) + epsi(r) lnk1 = ks(r) lnt1 = ts(r+1) lnk2 = savfun(lnk1,lnt1,psi,psisize,pok,pot) ks(r+1) = lnk2 C Compute some statistics immediately kmins = dmin1(ks(r+1),kmins) tmins = dmin1(ts(r+1),tmins) kmaxs = dmax1(ks(r+1),kmaxs) tmaxs = dmax1(ts(r+1),tmaxs) kmeans = kmeans + lnk2 300 CONTINUE kmeans = kmeans/dble(simtime) C 4.2 Evaluate the simulation C --------------------------- print *, ' ' print *,'Log of Steady state capital at ',ks(1) print *,'Mean Capital over the simulation',kmeans print *,'Std.Dev of shock is ',tsig 301 format(7X,4(A11,2X)) 302 format(1X,A4,2X,4(F11.4,2X)) 303 format(7X,2(A24,2x)) print *, ' ' print 303, 'Log Capital Stock','Log Product. Shock' print 301, 'Simulation','Grid','Simulation','Grid' print 302, 'Mins ',kmins, kmin, tmins, tmin print 302, 'Maxs ',kmaxs, kmax, tmaxs, tmax print *, ' ' C This ends the simulation if ENDIF C 5. Check the Accuracy C ==================== IF (accur.EQ.1) THEN C C Vary the capital stock C ---------------------- diffck = -1000000000.0D0 diffct = -1000000000.0D0 tobs = dint(dble(tnod)/2.0D0) obs = actime tgfix = tg(tobs) lnt1fix = scalup(tg(tobs),tmin,tmax) incr = (kmax-kmin)/(dble(obs)-1.0D0) lnk1v = kmin-incr DO 500 r = 1,obs lnk1v = lnk1v + incr lnk2v = savfun(lnk1v,lnt1fix,psi,psisize,pok,pot) actc = numi(lnk2v,lnt1fix,psi,psisize,pok,pot) c c calculate log consumption level using numerically c calculated conditional expectation and calculate c log consumption level using the approximation c actc = -log(actc)/nu appc = confun(lnk1v,lnt1fix,psi,psisize,pok,pot) diffcr = dabs(appc - actc) diffck = max(diffck,diffcr) 500 CONTINUE C Now we vary the technology shock C -------------------------------- kobs = dint(dble(knod)/2.0D0) kgfix = kg(kobs) lnk1fix = scalup(kgfix,kmin,kmax) incr = (tmax-tmin)/(dble(obs)-1.) lnt1v = tmin-incr DO 700 r = 1,obs lnt1v = lnt1v + incr lnk2v = savfun(lnk1fix,lnt1v,psi,psisize,pok,pot) actc = numi(lnk2v,lnt1v,psi,psisize,pok,pot) c c calculate log consumption level using numerically c calculated conditional expectation and calculate c log consumption level using the approximation c actc = -log(actc)/nu appc = confun(lnk1fix,lnt1v,psi,psisize,pok,pot) diffcr = dabs(appc - actc) diffct = max(diffct,diffcr) 700 CONTINUE print *,'Maximum percentage difference in consumption' print *,'varying Capital',diffck print *,'varying Shock ',diffct print *,'' C This ends the accuracy if ENDIF STOP END C ************************************************************* C Below are mathematical routines that are necessary to run pea2.f. C ************************************************************* C This is a subroutine which multiplies two matrices SUBROUTINE MMULT(X,Y,Z,N,M,K) C Argument Declarations IMPLICIT NONE INTEGER N,M,K DOUBLE PRECISION X(N,M),Y(M,K),Z(N,K) C Local Declarations INTEGER R,C,CC DOUBLE PRECISION ZRC C Loop over each matrix - element DO 10 R = 1,N DO 20 C = 1,K ZRC = 0 DO 30 CC = 1,M ZRC = ZRC + X(R,CC) * Y(CC,C) 30 CONTINUE Z(R,C) = ZRC 20 CONTINUE 10 CONTINUE RETURN END C ************************************************************* C C ************************************************************* C This is a subroutine which transposes a matrix SUBROUTINE MTRANS(X,XT,N,M) C Argument Declarations IMPLICIT NONE INTEGER N,M DOUBLE PRECISION X(N,M),XT(M,N) C Local Declarations INTEGER R,C C Loop over each matrix - element DO 10 R = 1,N DO 20 C = 1,M XT(C,R) = X(R,C) 20 CONTINUE 10 CONTINUE RETURN END C ************************************************************* C C ************************************************************* C This is a subroutine which inverts a DIAGONAL matrix SUBROUTINE MDINV(X,N,XI) C Argument Declarations IMPLICIT NONE INTEGER N DOUBLE PRECISION X(N,N),XI(N,N) C Local Declarations INTEGER R,C C Loop over each matrix - element DO 10 R = 1,N DO 20 C = 1,N IF (r.EQ.c) THEN XI(R,C) = 1/X(C,R) ELSE XI(R,C) = 0 ENDIF 20 CONTINUE 10 CONTINUE RETURN END C ************************************************************* C C ************************************************************* SUBROUTINE MINV(A,N,D) C C THIS SUBROUTINE COMPUTES THE INVERSE OF A MATRIX. 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 C ************************************************************* C C ************************************************************* C This is a subroutine which multiplies a matrix and a vector SUBROUTINE MVMULT(X,Y,Z,N,M) C Argument Declarations IMPLICIT NONE INTEGER N,M DOUBLE PRECISION X(N,M),Y(M),Z(N) C Local Declarations INTEGER R,C DOUBLE PRECISION ZR C Loop over each matrix - element DO 10 R = 1,N ZR = 0.0D0 DO 20 C = 1,M ZR = ZR + X(R,C) * Y(C) 20 CONTINUE Z(R) = ZR 10 CONTINUE RETURN END C ************************************************************* C C ************************************************************* C This is a subroutine which multiplies two vectors SUBROUTINE VMULT(X,Y,Z,N) C Argument Declarations IMPLICIT NONE INTEGER N DOUBLE PRECISION X(*),Y(*),Z C Local Declarations INTEGER R C Loop over each matrix - element Z = 0.D0 DO 10 R = 1,N Z = Z + X(R) * Y(R) 10 CONTINUE RETURN END C ************************************************************* C C ************************************************************* C This is a subroutine which computes the Euclidean norm (2norm) C of the distance between two vectors C REAL*8 FUNCTION NORM(X,Y,N) C Argument Declarations IMPLICIT NONE INTEGER N DOUBLE PRECISION X(N),Y(N) C Local Declarations INTEGER R DOUBLE PRECISION Z C Loop over each matrix - element NORM = 0 DO 10 R = 1,N Z = X(R) - Y(R) NORM = NORM + Z**2.0D0 10 CONTINUE NORM = NORM**(0.5D0) RETURN END C ************************************************************* c C THE FOLLOWING TWO FUNCTIONS GENERATE RANDOM NUMBERS C C ************************************************************* C REAL*8 FUNCTION GASDEV(IDUM) c c returns a normally distributed deviate with zero mean and unit c variance, using RAN1(IDUM) as the source of uniform deviates c 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 C ************************************************************* C C ************************************************************* REAL*8 FUNCTION RAN1(IDUM) c c returns a uniform random deviate 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 ************************************************************* C ************************************************************* C Below are the PEA related routines that are C necessary to run pea3.f. C ************************************************************* C ************************************************************* C C Create an ORD-th ORDer Chebyshev Polynomial C of the input vector p C SUBROUTINE CHEBPOL(ord,p,n,x) C Argument Declarations IMPLICIT none INTEGER ord,n DOUBLE PRECISION p(n), x(n,ord+1) C Local Declarations INTEGER r,c DO 05 r = 1,n x(r,1) = 1.0D0 x(r,2) = p(r) 05 CONTINUE IF (ord .ge. 2) then DO 10 c = 3,(ord+1) DO 20 r = 1,n x(r,c) = 2.0D0 * p(r) * x(r,c-1) + - x(r,c-2) 20 CONTINUE 10 CONTINUE ENDIF RETURN END C *********************************************************** C C *********************************************************** C Create n Chebyshev nodes SUBROUTINE CHEBNODE(N,X) C Argument Declarations IMPLICIT NONE INTEGER N DOUBLE PRECISION X(N) C Local Declarations INTEGER R DOUBLE PRECISION PI PARAMETER (PI=3.14159265D0) DO 10 R = 1,N X(R) = COS( ( PI * (R-0.5D0) )/dble(N) ) 10 CONTINUE RETURN END C *********************************************************** C C *********************************************************** C This function linearly transforms the (N,1) vector X C into [-1,1] C REAL*8 FUNCTION SCALDOWN(X,XMIN,XMAX) C Argument Declarations IMPLICIT NONE DOUBLE PRECISION X,XMIN,XMAX C Local Declarations DOUBLE PRECISION A,B,DENOM DENOM = XMAX - XMIN A = 2.0D0 / DENOM B = 1.0D0 - ( 2.0D0 * (XMAX / DENOM) ) SCALDOWN = B + A * X RETURN END C ************************************************************* C *********************************************************** C C This function linearly transforms the input C from [-1,1] into [XMIN,XMAX] C REAL*8 FUNCTION SCALUP(X,XMIN,XMAX) C Argument Declarations DOUBLE PRECISION X,XMIN,XMAX C Local Declarations DOUBLE PRECISION A,B A = XMIN + XMAX B = XMAX - XMIN SCALUP = ( A + B*X )/2.0D0 RETURN END C *********************************************************** C C *********************************************************** C returns the ln of consumption C REAL*8 FUNCTION confun(lnk1,lnt1,psi,psisize,pok,pot) IMPLICIT NONE C Common Declarations COMMON /consfcn/ nu,tmin,tmax,kmin,kmax INTEGER psisize,pok,pot DOUBLE PRECISION nu,tmin,tmax,kmin,kmax C Argument Declarations DOUBLE PRECISION lnk1,lnt1,psi(psisize) C Local Declarations INTEGER ki,ti,psiind DOUBLE PRECISION lnk1u,lnt1u,pol,scaldown DOUBLE PRECISION xk0,xk1,xk,xt0,xt1,xt,xmat lnk1u = scaldown(lnk1,kmin,kmax) lnt1u = scaldown(lnt1,tmin,tmax) C Prepare the Polynomial xk0 = 1.0D0 xk1 = lnk1u pol = 0.0D0 DO 12 ki= 0,pok if (ki .eq. 0) then xk = 1.0D0 elseif (ki .eq. 1) then xk = lnk1u else xk = 2.0D0 * lnk1u * xk1 - xk0 xk0 = xk1 xk1 = xk endif xt0 = 1.0D0 xt1 = lnt1u DO 13 ti = 0,pot if (ti .eq. 0) then xt = 1.0D0 elseif (ti .eq. 1) then xt = lnt1u else xt = 2.0D0 * lnt1u * xt1 - xt0 xt0 = xt1 xt1 = xt endif xmat = xk * xt psiind = ki * (pot+1) + ti + 1 pol = pol + xmat * psi(psiind) 13 CONTINUE 12 CONTINUE confun = -1.0D0/nu * pol RETURN END C *********************************************************** C C *********************************************************** C returns the ln of new capital C REAL*8 FUNCTION savfun(lnk1,lnt1,psi,psisize,pok,pot) IMPLICIT NONE C Common Declarations COMMON /consfcn/ nu,tmin,tmax,kmin,kmax COMMON /savfcn/ alpha,delta,dfactor INTEGER psisize,pok,pot DOUBLE PRECISION alpha,delta,dfactor,nu,tmin,tmax,kmin,kmax C Argument Declarations DOUBLE PRECISION lnk1,lnt1,lnk2,psi(psisize) C Local Declarations DOUBLE PRECISION lnc1,c1,k1,k2,t1,confun lnc1 = confun(lnk1,lnt1,psi,psisize,pok,pot) c1 = exp(lnc1) k1 = exp(lnk1) t1 = exp(lnt1) k2 = t1 * (k1**alpha) + (1-delta) * k1 - c1 IF ( k2.le.0 ) THEN write(*,*) 'consumption too large' stop ENDIF lnk2 = log(k2) savfun = lnk2 RETURN END C *********************************************************** C C *********************************************************** C Compute the actual marginal utility C SUBROUTINE aexp(epsi,lnk2,lnt1,psi,aex,psisize,pok,pot) C Common Declarations INTEGER psisize,pok,pot COMMON /consfcn/ nu,tmin,tmax,kmin,kmax COMMON /savfcn/ alpha,delta,dfactor COMMON /aexpfcn/ rho,sigma DOUBLE PRECISION alpha,delta,dfactor,nu,tmin,tmax,kmin,kmax, + rho,sigma C Argument Declarations DOUBLE PRECISION epsi,lnk2,lnt1,aex,psi(psisize) C Local Declarations DOUBLE PRECISION lnt2,lnc2,c2,k2,t2,x,confun lnt2 = rho * lnt1 + epsi lnc2 = confun(lnk2,lnt2,psi,psisize,pok,pot) c2 = exp(lnc2) k2 = exp(lnk2) t2 = exp(lnt2) x = alpha * t2 * (k2**(alpha-1)) + 1 - delta aex = dfactor * (c2**(-nu)) * x RETURN END C *********************************************************** C C *********************************************************** C Calculate the expected value of E(aexp(epsi)) using C Hermite Gaussian quadrature when epsi is a normally C distributed random variable with standard deviations sigma C REAL*8 FUNCTION numi(lnk2,lnt1,psi,psisize,pok,pot) C Common Declarations IMPLICIT NONE COMMON /consfcn/ nu,tmin,tmax,kmin,kmax COMMON /savfcn/ alpha,delta,dfactor COMMON /aexpfcn/ rho,sigma INTEGER psisize,pok,pot DOUBLE PRECISION alpha,delta,dfactor,nu,tmin,tmax,kmin,kmax, + rho,sigma C Argument Declarations DOUBLE PRECISION lnt1,lnk2 DOUBLE PRECISION psi(psisize) C Local Declarations INTEGER i, tnod PARAMETER (tnod = 11) DOUBLE PRECISION weight(tnod),nodes(tnod),nodet,valu(tnod) DOUBLE PRECISION s2,pi,spi,valui PARAMETER (pi=3.14159265D0) C Define the nodes and weights for the integration CALL gauher(nodes,weight,tnod) C C Note that there is a transformation of variables because C we have a normal density but use Hermite Gaussian Quadrature. C The Jacobian is 2**0.5*sig. The Jacobian cancels against C the 2**0.5*sig that is in the normal density C s2 = 2**.5 DO 10 i = 1,tnod nodet = sigma * s2 * nodes(i) CALL aexp(nodet,lnk2,lnt1,psi,valui,psisize,pok,pot) valu(i) = valui 10 CONTINUE C Compute the Integral spi = dble(pi**.5) numi = 0.0D0 CALL vmult(weight,valu,numi,tnod) numi = (1.0D0/spi) * numi RETURN END C *********************************************************** C C *********************************************************** SUBROUTINE gauher(x,w,n) C Compute the nodes and weights for Gauss-Hermite quadrature C Argument Declarations INTEGER n DOUBLE PRECISION W(n),X(n) C Local Declarations INTEGER i,its,j,m,MAXIT DOUBLE PRECISION EPS,PIM4,P1,P2,P3,PP,Z,Z1 PARAMETER (EPS=3.D-14,PIM4=.7511255444649425,MAXIT=10) m = dint(dble(n+1)/2.0D0) DO 10 i = 1,m IF (i.EQ.1) THEN z = sqrt(float(2*n+1))-1.85575 + * (2*n+1)**(-.16666667) ELSE IF (i.EQ.2) THEN z = z - 1.14 * n **.426/z ELSE IF (i.EQ.3) THEN z = 1.86 * z - .86 * x(1) ELSE IF (i.EQ.4) THEN z = 1.91 * z - .91 * x(2) ELSE z = 2 .* z - x(i-2) ENDIF DO 20 its = 1,MAXIT p1 = PIM4 p2 = 0.d0 DO 30 j = 1,n p3 = p2 p2 = p1 p1 = z * sqrt(2.d0/j) * p2 - + sqrt(dble(j-1)/dble(j))*p3 30 CONTINUE pp = sqrt(2.d0*n) * p2 z1 = z z = z1 - p1/pp IF (abs(z-z1) .le. EPS) THEN GOTO 1 ENDIF 20 CONTINUE pause 'too many iterations in gauher' 1 x(i) = z x(n+1-i)= -z w(i) = 2.d0/(pp*pp) w(n+1-i)= w(i) 10 CONTINUE RETURN END C *********************************************************** C C *********************************************************** C C Key Characteristics of the algorithm: C C 1. State variables are the logs of capital,ln(k), and C productivity,ln(t). C 2. The conditional expectation is approximated using C exp{ Pn[ln(k),ln(t)]}, where Pn[x] is a Chebyshev polynomial C of order n in x. C 3. Since Chebyshev polynomials are defined on the interval [-1,1] C the state variables are scaled to lie in this interval. C 4. The conditional expectation is calculated numerically using C Hermite Gaussian quadrature. C 5. A functional iteration scheme is used to find the parameters C of the approximating function. C C Steps of the algorithm: C C I. Choose initial values psi0 for the coefficients of the C approximating function. Note that for given values of psi C it is easy to calculate the values for the consumption and C capital choice as a function of the beginning-of-period C capital stock and the productivity shock by replacing the C true conditional expectation by the approximating function C in the Euler equation. C C II. At each grid point numerically calculate the conditional C expectation of C C dfactor*{confun[ln(k'),rho*ln(t)+sigma*epsi']} C *{alpha*(k')^(alpha-1)*exp(rho*ln(t)+sigma*innovation')+1-delta} C C where a ' indicates next period's value. C Note that this is a function of the random variable epsilon', C which is assumed to be a N(0,1) random variable so we can use C Hermite Gaussian quadrature. The conditional expectation is C also a function of ln(t) which is known at each grid point C and ln(k') which is a function of ln(t) and ln(k) and has C to be calculated before the conditional expectation is C calculated. psi0 is used to calculate next period's capital C stock, savfun[ln(k),ln(t)], next period's consumption level C for different values of the innovation epsi, C confun[ln(k'),ln(t')]. C C III. Project the log of the calculated conditional expectations C on Pn[ln(k),ln(t)]. Note that this is just a simple C linear regression and that taking logs is allowed C since there is no stochastic error term in the equation C that equates the conditional expectation to the approximating C function. C C IV. If the values for psi that come out of step 3 are close to C the ones used in step 2 then the algorithm has converged, C if not step 2 and 3 have to be repeated. C C C Functions used: C C 1.confun(lnk1,lnt1,psi): calculates the log of the consumption C level as a function of the unscaled levels of the state variables. C 2.savfun(lnk1,lnt1,psi): calculates the log of the capital stock C as a function of the unscaled levels of the state variables. C 3.aexp(epsi,lnk2,lnt1,psi): calculates the discounted marginal C utility of next period's marginal product of capital as a function C of the log of next period's capital stock (which is known today), C the log of today's productivity level, and next period's innovation C 4.scalup(x,xmin,xmax): a linear function to map a variable defined C on [-1,1] to [xmin,xmax]. C 5.scaldown(x,xmin,xmax): a linear function to map a variable defined C on [xmin,xmax] to [-1,1]. C 6.chebnode(n): function to define n Chebyshev nodes. C 7.chebpol(ord,p): function to define a (p+1) vector with C a constant and the p Chebyshev polynomial terms. C 8.makepoly(ord,x): For the (Tx2) matrix x and the (1x2) vector ord C this creates a [T x ord(1)*ord(2)] matrix where each row contains C the ord(1)*ord(2) bivariate Chebyshev polynomial terms. C 9.hernodes(n): function to calculate the (nx2) matrix with the C n Hermite nodes and weights. C 10numi(fun,lnk2,lnt1,psi): function to numerically integrate the C function fun(epsi) using Hermite Gaussian quadrature when C epsi is a normally distributed random variable with zero mean C and standard deviation sigma. C C C 1.2 Initialize algorithm parameters and algorithm options C --------------------------------------------------------- C C kmin: minimum capital stock as a fraction of steady state C capital stock C kmax: maximum capital stock as a fraction of steady state C capital stock C tmin: minimum log of productivity shock in number of C standard deviations away from 0 C tmax: maximum log of productivity shock in number of C standard deviations away from 0 C C Note that kmin,kmax,tmin, and tmax will be redefined in the C program as the actual boundaries of the grid for log capital C and log productivity shock C C maxiter: maximum number of iterations to reach fixed point C psitol: accuracy criterion C lrate: updating coefficient from one iteration to the next C C idum: seed for simulation (simulation not used for C solution procedure) C C simyes:if simyes = 1 then the program simulates a timeseries C accur: if accur = 1 then the program performs an accuracy test C C 2.1 Create the grid C ------------------- C All gridpoints are stored in two vectors kgrid,tgrid with dimensions C (knod*tnod) by 1. C Note that the gridpoints are just the standard Chebyshev nodes which C are in the interval [-1,1]. The subroutine scalup will tell you what C the value of the state variable is for any point in [-1,1]. C C C 2.2 Create the Projection Matrix C -------------------------------- C In step III of the algorithm we regress the calculated conditional C expectation terms on the polynomial terms. Since the regressors never C change we can compute and store the X matrix of the regression and C the inv(X'X)*X' matrix. C C 3.1 Compute the expected value (Step II) C ---------------------------------------- C Here we calculate the conditional expectation at each grid C point. First we calculate next period's capital stock as a C function of the ln of the current capital stock and the ln C of the technology shock. Second we numerically integrate the C function aexp over the random innovation to the technology C shock. psi is an input of the numerical integration since C it is used to evaluate the consumption at lnk2 and lnt2. C C 3.2 Update the polynomial coefficients (Steps III and IV) C --------------------------------------------------------- C Here we calculate new coefficients for the parametrized C expectation by regressing the ln of yvec(i) on the C polynomial terms. Then we measure the distance between the C new and the old psi and if the difference is not small C enough then we take a weighted average and repeat the loop. C C 4. Simulation C ============= C Here we show how the solution can be used to generate a C time-series. We also calculate the mean capital stock and some C extreme values to check whether the grid boundaries make any sense C C 5. Accuracy C =========== C Here we perform an accuracy test. In particular at a bunch of points we C calculate the approximation of the conditional expectation and we C numerically calculate the conditional expectation. Then we compare the C consumption levels implied by these two conditinal expectations. In C practice you would want to cover the whole state space where you want C your solution to be accurate. Here we just vary one state variable and C keep the other state variable at a grid point value. Note that when C knod = pok+1 and tnod = pot+1 then by construction you get an exact fit C at the grid points. C C============================================================================ C C EXAMPLE; Solution for psi when C C pok = 3, pot = 3, knod = 7, tnod = 7 C C alpha = 0.33D0 C dfactor = 0.95D0 C nu = 2.0D0 C delta = .95D0 C C sigma = .01D0 C rho = 0.95D0 C C C kmin = 0.8D0 C kmax = 1.2D0 C C tmin = -3.5D0 C tmax = 3.5D0 C C C 1.856033636012942 C -2.357064133658491E-001 C -2.439818435654846E-005 C 1.206701028617710E-006 C -1.136367830054171E-001 C 4.039760903712409E-005 C -8.929668254395991E-006 C 1.590293624380350E-007 C 3.307958486399994E-005 C 1.066761973693250E-005 C -2.725911314969676E-007 C 3.537088810880155E-009 C -1.998446325740655E-006 C 2.363115533506463E-007 C -5.022199010994641E-009 C 4.228266668466543E-011 C C