C Solves Krusell and Smith (1998) economy as in Young (2005) C Uses piecewise-histogram approximation to parameterize cross-sectional distribution. C Uses hybrid cubic spline-polynomial interpolation scheme to parameterize value function. C To run: compile with links to helph36.f, matrix4.f, minim.f, random2.f, and statfunc.f C Required input files: xkpts150.in, saving.f1 C Hot start input files: jedc2.out3, jedc2.out C Written: ERY, June 2002 C Last edited: ERY, April 2005 implicit real*8 (a-h,o-z) parameter (nk2pts=5000,nkpts=150,nread=1,nsample=1000, & nmupts=4,ncapt=10000,nm2pts=30,durug=1.5D+00, & unempg=0.04D+00,durgd=8.0D+00,unempb=0.1D+00, & durbd=8.0D+00,zgood=1.01D+00,zbad=0.99D+00, & durub=2.5D+00,alpha=0.36D+00,nread3=1,ndod=0, & hfix=0.3271D+00,toler=1.0D-05, & theta=0.25D+00,whome=0.0D+00,mxiter=10) common/cfs/f1(nk2pts,2) common/cv2/y2(nkpts,nmupts,4),optk12(nkpts,nmupts,2), & optk02(nkpts,nmupts,2) common/cprice/r(nmupts,2),w(nmupts,2) common/cgrid/xkpts(nkpts),xmupts(nmupts) common/ccoef/a(2),b(2),anew(2),bnew(2) common/cprob/pi(4,4) common/nshock/nz(ncapt+1) common/cmgrid/xm2pts(nm2pts) common/ckgrid/xk2pts(nk2pts) common/cdecg/optkg1(nk2pts,nm2pts,2),optkg0(nk2pts,nm2pts,2) allocatable :: uniz(:),hansens(:),dseeds(:) one = 1.0D+00 c Set Markov chain for (z,epsilon) pgg00 = (durug-one)/durug pbb00 = (durub-one)/durub pbg00 = 1.25D+00*pbb00 pgb00 = 0.75D+00*pgg00 pgg01 = (unempg-unempg*pgg00)/(one-unempg) pbb01 = (unempb-unempb*pbb00)/(one-unempb) pbg01 = (unempb-unempg*pbg00)/(one-unempg) pgb01 = (unempg-unempb*pgb00)/(one-unempb) pgg = (durgd-one)/durgd pgb = one-(durbd-one)/durbd pgg10 = one-(durug-one)/durug pbb10 = one-(durub-one)/durub pbg10 = one-1.25D+00*pbb00 pgb10 = one-0.75D+00*pgg00 pgg11 = one-(unempg-unempg*pgg00)/(one-unempg) pbb11 = one-(unempb-unempb*pbb00)/(one-unempb) pbg11 = one-(unempb-unempg*pbg00)/(one-unempg) pgb11 = one-(unempg-unempb*pgb00)/(one-unempb) pbg = one-(durgd-one)/durgd pbb = (durbd-one)/durbd pi(1,1) = pgg*pgg11 pi(1,2) = pgg*pgg01 pi(1,3) = pbg*pbg11 pi(1,4) = pbg*pbg01 pi(2,1) = pgg*pgg10 pi(2,2) = pgg*pgg00 pi(2,3) = pbg*pbg10 pi(2,4) = pbg*pbg00 pi(3,1) = pgb*pgb11 pi(3,2) = pgb*pgb01 pi(3,3) = pbb*pbb11 pi(3,4) = pbb*pbb01 pi(4,1) = pgb*pgb10 pi(4,2) = pgb*pgb00 pi(4,3) = pbb*pbb10 pi(4,4) = pbb*pbb00 write(6,592) ((pi(i,j),j=1,4),i=1,4) 592 format(4f18.10) c Draw random z sequence allocate (uniz(ncapt+1)) dseed = 1894.0D+00 call ggubs(dseed,ncapt+1,uniz) nz(1) = 1 do 52 i = 1,ncapt if (nz(i) .eq. 1) then if (uniz(i) .le. pgg) then nz(i+1) = 1 else nz(i+1) = 2 endif else if (uniz(i) .le. pgb) then nz(i+1) = 1 else nz(i+1) = 2 endif endif 52 continue deallocate (uniz) c Set grids and determine prices open(2,file='xkpts150.in',form='formatted') do 999 i = 1,nkpts read(2,570) xkpts(i) 570 format(f18.10) 999 continue close(2) xk2pts(1) = xkpts(1) xkinc = (400.0D+00-xkpts(1))/dble(real(nk2pts-1)) do 393 i = 2,nk2pts xk2pts(i) = xk2pts(i-1)+xkinc 393 continue xmupts(1) = 10.0D+00 xmupts(2) = 10.75D+00 xmupts(3) = 11.5D+00 xmupts(4) = 12.25D+00 xminc = (xmupts(4)-xmupts(1))/dble(real(nm2pts-1)) xm2pts(1) = xmupts(1) do 58 i = 2,nm2pts xm2pts(i) = xm2pts(i-1)+xminc 58 continue do 38 i = 1,nmupts r(i,1) = alpha*zgood*xmupts(i)**(alpha-one)* & ((one-unempg)*hfix)**(one-alpha) w(i,1) = (one-alpha)*zgood*xmupts(i)**alpha* & ((one-unempg)*hfix)**(-alpha) r(i,2) = alpha*zbad*xmupts(i)**(alpha-one)* & ((one-unempb)*hfix)**(one-alpha) w(i,2) = (one-alpha)*zbad*xmupts(i)**alpha* & ((one-unempb)*hfix)**(-alpha) write(6,787) r(i,1),w(i,1),r(i,2),w(i,2) 787 format(4f18.10) 38 continue c Input initial conditions if (nread .eq. 1) then open(2,file='jedc2.out',form='formatted') do 98 i = 1,nkpts do 97 j = 1,nmupts read(2,700) dum,dum,y2(i,j,1),y2(i,j,2), & y2(i,j,3),y2(i,j,4), & optk12(i,j,1),optk12(i,j,2), & optk02(i,j,1),optk02(i,j,2) 700 format(10f18.10) 97 continue 98 continue close(2) else open(2,file='xkpts150.in',form='formatted') do 57 i = 1,nkpts read(2,390) xkpts(i) 390 format(f18.10) 57 continue close(2) do 56 i = 1,nkpts do 96 j = 1,nmupts y2(i,j,1) = dlog(xkpts(i)+ & (one-alpha)/alpha* & xmupts(j)/(one-unempg)) y2(i,j,2) = dlog(xkpts(i)+ & whome/r(j,1)) y2(i,j,3) = dlog(xkpts(i)+ & (one-alpha)/alpha* & xmupts(j)/(one-unempb)) y2(i,j,4) = dlog(xkpts(i)+ & whome/r(j,2)) optk12(i,j,1) = xkpts(i) optk12(i,j,2) = xkpts(i) optk02(i,j,1) = xkpts(i) optk02(i,j,2) = xkpts(i) 96 continue 56 continue endif c Iterate over coefficients for fixed law of motion if (nread3 .eq. 1) then open(2,file='jedc2.out3',form='formatted') read(2,303) (a(i),i=1,2) read(2,303) (b(i),i=1,2) 303 format(2f18.10) close(2) else a(1) = 0.093D+00 a(2) = 0.965D+00 b(1) = 0.084D+00 b(2) = 0.966D+00 endif iter = 1 1 call fvec(a,b,anew,bnew) write(6,303) (a(i),i=1,2) write(6,303) (anew(i),i=1,2) write(6,303) (b(i),i=1,2) write(6,303) (bnew(i),i=1,2) iter = iter+1 diffmx = 0.0D+00 do 74 i = 1,2 adiff = dabs(anew(i)-a(i)) if (adiff .gt. diffmx) diffmx = adiff bdiff = dabs(bnew(i)-b(i)) if (bdiff .gt. diffmx) diffmx = bdiff 74 continue if (diffmx .lt. toler .or. iter .gt. mxiter) goto 2 do 63 i = 1,2 a(i) = (one-theta)*a(i)+theta*anew(i) b(i) = (one-theta)*b(i)+theta*bnew(i) 63 continue open(9,file='jedc2.tmp',form='formatted') write(9,303) (a(i),i=1,2) write(9,303) (b(i),i=1,2) close(9) goto 1 2 open(9,file='jedc2.out3',form='formatted') write(9,303) (a(i),i=1,2) write(9,303) (b(i),i=1,2) close(9) c Optional loop that compute distribution of den Haan-Marcet statistic c Usually skip because simulation length is so long that test is always rejected if (ndod .eq. 1) then allocate (dseeds(nsample),hansens(nsample)) eseed = 34958.0D+00 call ggubs(eseed,nsample,dseeds) hansens = 0.0D+00 do 39 j = 1,nsample dseed = dseeds(j)*1000.0D+00 call ggubs(dseed,ncapt+1,uniz) nz(1) = 1 do 53 i = 1,ncapt if (nz(i) .eq. 1) then if (uniz(i) .le. pgg) then nz(i+1) = 1 else nz(i+1) = 2 endif else if (uniz(i) .le. pgb) then nz(i+1) = 1 else nz(i+1) = 2 endif endif 53 continue call calcf(optkg1,optkg0,xk2pts,xm2pts,a,b,anew,bnew, & xkmean,sdmean,skmean,s4mean,rsqa,rsqb, & sderra,sderrb, & gmean,p1mean,p5mean,p10mean,p20mean, & p30mean,hansen,erfmx1,erfmx2) hansens(j) = hansen if (multpl(j,10) .eq. 1) print *,j 39 continue open(2,file='jedc2.hans',form='formatted') do 88 i = 1,nsample write(2,366) hansens(i) 366 format(f18.10) 88 continue close(2) deallocate (hansens,dseeds) endif end c Solves value function, simulates, computes implied coefficients subroutine fvec(a,b,anew,bnew) implicit real*8 (a-h,o-z) parameter (nkpts=150,alpha=0.36D+00,delta=0.025D+00, & hfix=0.3271D+00,xkbor=0.0D+00,beta=0.99D+00, & mxiter=2000,nk2pts=5000,nm2pts=30,gamma=1.0D+00, & mxpol=200,nmupts=4,unempb=0.1D+00, & unempg=0.04D+00,zgood=1.01D+00,zbad=0.99D+00, & whome=0.0D+00) common/cfs/f1(nk2pts,2) common/cv2/y2(nkpts,nmupts,4),optk12(nkpts,nmupts,2), & optk02(nkpts,nmupts,2) common/cprice/r(nmupts,2),w(nmupts,2) common/cgrid/xkpts(nkpts),xmupts(nmupts) common/cvs/v11(nkpts,nmupts,2),v11dp(nkpts,nmupts,2), & v10(nkpts,nmupts,2),v10dp(nkpts,nmupts,2), & v01(nkpts,nmupts,2),v01dp(nkpts,nmupts,2), & v00(nkpts,nmupts,2),v00dp(nkpts,nmupts,2) common/cws/w11(nkpts,2),w11dp(nkpts,2), & w10(nkpts,2),w10dp(nkpts,2), & w01(nkpts,2),w01dp(nkpts,2), & w00(nkpts,2),w00dp(nkpts,2) common/cprob/pi(4,4) common/cmgrid/xm2pts(nm2pts) common/ckgrid/xk2pts(nk2pts) common/cvgrid/vgrid(nk2pts,nm2pts,4) common/cdecg/optkg1(nk2pts,nm2pts,2),optkg0(nk2pts,nm2pts,2) allocatable :: y(:,:,:),optk1(:,:,:),optk0(:,:,:) dimension a(2),b(2),anew(2),bnew(2),r0(2),w0(2) one = 1.0D+00 xkcut1 = 0.5D+00 xkcut2 = 0.5D+00 xklow = xkbor xkhigh = 1255.0D+00 vmxdif = 100.0D+00 c Value function iteration allocate (y(nkpts,nmupts,4),optk1(nkpts,nmupts,2), & optk0(nkpts,nmupts,2)) do 30 kkk = 1,mxiter if (vmxdif .lt. 1.0D-00) then ncalop = 0 else ncalop = 1 endif y = y2 optk1 = optk12 optk0 = optk02 call interp(xkpts,xmupts,y,v11,v10,v01,v00, & v11dp,v10dp,v01dp,v00dp,a,b) do 13 i = 1,nkpts do 567 j = 1,nmupts xkcur = xkpts(i) nmu = j neps = 1 nagg = 1 if (xkcur .le. xkcut1 .or. ncalop .eq. 1) then call calopt(xkcur,neps,xkp1,fval1,niter, & 1.0D-10,r,w,xkhigh,xklow, & grad1,xkpts,xmupts,nmu,nagg) else call altopt(xkcur,neps,xkp1,fval1,niter, & 1.0D-10,optk1(i,j,1),grad1, & xkpts,r,w,xmupts,nmu,nagg,ndone) if (ndone .eq. 0) then call calopt(xkcur,neps,xkp1,fval1,niter, & 1.0D-10,r,w,xkhigh,xklow, & grad1,xkpts,xmupts,nmu,nagg) endif endif optk12(i,j,1) = xkp1 y2(i,j,1) = fval1 neps = 0 nagg = 1 if (xkcur .le. xkcut2 .or. ncalop .eq. 1) then call calopt(xkcur,neps,xkp2,fval2,niter, & 1.0D-10,r,w,xkhigh,xklow, & grad0,xkpts,xmupts,nmu,nagg) else call altopt(xkcur,neps,xkp2,fval2,niter, & 1.0D-10,optk0(i,j,1),grad0, & xkpts,r,w,xmupts,nmu,nagg,ndone) if (ndone .eq. 1) then call calopt(xkcur,neps,xkp2,fval2,niter, & 1.0D-10,r,w,xkhigh,xklow, & grad0,xkpts,xmupts,nmu,nagg) endif endif optk02(i,j,1) = xkp2 y2(i,j,2) = fval2 neps = 1 nagg = 2 if (xkcur .le. xkcut1 .or. ncalop .eq. 1) then call calopt(xkcur,neps,xkp3,fval3,niter, & 1.0D-10,r,w,xkhigh,xklow, & grad1,xkpts,xmupts,nmu,nagg) else call altopt(xkcur,neps,xkp3,fval3,niter, & 1.0D-10,optk1(i,j,2),grad1, & xkpts,r,w,xmupts,nmu,nagg,ndone) if (ndone .eq. 1) then call calopt(xkcur,neps,xkp3,fval3,niter, & 1.0D-10,r,w,xkhigh,xklow, & grad1,xkpts,xmupts,nmu,nagg) endif endif optk12(i,j,2) = xkp3 y2(i,j,3) = fval3 neps = 0 nagg = 2 if (xkcur .le. xkcut2 .or. ncalop .eq. 1) then call calopt(xkcur,neps,xkp4,fval4,niter, & 1.0D-10,r,w,xkhigh,xklow, & grad0,xkpts,xmupts,nmu,nagg) else call altopt(xkcur,neps,xkp4,fval4,niter, & 1.0D-10,optk0(i,j,2),grad0, & xkpts,r,w,xmupts,nmu,nagg,ndone) if (ndone .eq. 0) then call calopt(xkcur,neps,xkp4,fval4,niter, & 1.0D-10,r,w,xkhigh,xklow, & grad0,xkpts,xmupts,nmu,nagg) endif endif optk02(i,j,2) = xkp4 y2(i,j,4) = fval4 567 continue 13 continue c Howard's improvement if (kkk .gt. 5) then npol = 1 else npol = 0 endif if (npol .eq. 0) goto 2 do 43 ii = 1,mxpol call interp(xkpts,xmupts,y2,v11,v10,v01,v00, & v11dp,v10dp,v01dp,v00dp,a,b) do 42 i = 1,nkpts do 41 j = 1,nmupts xkcur = xkpts(i) con1 = (r(j,1)+one-delta)*xkcur+w(j,1)*hfix- & optk12(i,j,1) call splint(xkpts,v11,v11dp,optk12(i,j,1), & yval1,yp, & ydp,khi,klo,0,nkpts,1,0,0,j,1) call splint(xkpts,v10,v10dp,optk12(i,j,1), & yval2,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,1) call splint(xkpts,v01,v01dp,optk12(i,j,1), & yval3,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,1) call splint(xkpts,v00,v00dp,optk12(i,j,1), & yval4,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,1) if (dabs(gamma-one) .lt. 1.0D-04) then utilc = dlog(con1) else utilc = one/(one-gamma)*(con1**(one-gamma)-one) endif y2(i,j,1) = utilc+beta*(pi(1,1)*yval1+ & pi(1,2)*yval2+pi(1,3)*yval3+ & pi(1,4)*yval4) con2 = (r(j,1)+one-delta)*xkcur-optk02(i,j,1)+ & whome call splint(xkpts,v11,v11dp,optk02(i,j,1), & yval1,yp, & ydp,khi,klo,0,nkpts,1,0,0,j,1) call splint(xkpts,v10,v10dp,optk02(i,j,1), & yval2,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,1) call splint(xkpts,v01,v01dp,optk02(i,j,1), & yval3,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,1) call splint(xkpts,v00,v00dp,optk02(i,j,1), & yval4,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,1) if (dabs(gamma-one) .lt. 1.0D-04) then utilc = dlog(con2) else utilc = one/(one-gamma)*(con2**(one-gamma)-one) endif y2(i,j,2) = utilc+beta*(pi(2,1)*yval1+ & pi(2,2)*yval2+pi(2,3)*yval3+ & pi(2,4)*yval4) con3 = (r(j,2)+one-delta)*xkcur+w(j,2)*hfix- & optk12(i,j,2) call splint(xkpts,v11,v11dp,optk12(i,j,2), & yval1,yp, & ydp,khi,klo,0,nkpts,1,0,0,j,2) call splint(xkpts,v10,v10dp,optk12(i,j,2), & yval2,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,2) call splint(xkpts,v01,v01dp,optk12(i,j,2), & yval3,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,2) call splint(xkpts,v00,v00dp,optk12(i,j,2), & yval4,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,2) if (dabs(gamma-one) .lt. 1.0D-04) then utilc = dlog(con3) else utilc = one/(one-gamma)*(con3**(one-gamma)-one) endif y2(i,j,3) = utilc+beta*(pi(3,1)*yval1+ & pi(3,2)*yval2+pi(3,3)*yval3+ & pi(3,4)*yval4) con4 = (r(j,2)+one-delta)*xkcur-optk02(i,j,2)+ & whome call splint(xkpts,v11,v11dp,optk02(i,j,2), & yval1,yp, & ydp,khi,klo,0,nkpts,1,0,0,j,2) call splint(xkpts,v10,v10dp,optk02(i,j,2), & yval2,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,2) call splint(xkpts,v01,v01dp,optk02(i,j,2), & yval3,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,2) call splint(xkpts,v00,v00dp,optk02(i,j,2), & yval4,yp, & ydp,khi,klo,1,nkpts,1,0,0,j,2) if (dabs(gamma-one) .lt. 1.0D-04) then utilc = dlog(con4) else utilc = one/(one-gamma)*(con4**(one-gamma)-one) endif y2(i,j,4) = utilc+beta*(pi(4,1)*yval1+ & pi(4,2)*yval2+pi(4,3)*yval3+ & pi(4,4)*yval4) 41 continue 42 continue 43 continue c Check for convergence using sup-norm 2 vmxdif = 0.0D+00 dmxdif = 0.0D+00 emxdif = 0.0D+00 do 15 i = 1,nkpts do 14 j = 1,nmupts do 80 k = 1,4 vdiff = dabs(y2(i,j,k)-y(i,j,k)) if (vdiff .gt. vmxdif) vmxdif = vdiff 80 continue do 79 k = 1,2 ddiff = dabs(optk12(i,j,k)-optk1(i,j,k)) if (ddiff .gt. dmxdif) dmxdif = ddiff ediff = dabs(optk02(i,j,k)-optk0(i,j,k)) if (ediff .gt. emxdif) emxdif = ediff 79 continue 14 continue 15 continue if (vmxdif .lt. 1.0D-06 .and. dmxdif .lt. 1.0D-08 & .and. emxdif .lt. 1.0D-08 .and. kkk .gt. 2) goto 8 if (multpl(kkk,100) .eq. 1) then write(6,888) kkk,mxiter,vmxdif,dmxdif,emxdif 888 format(2i5,3f15.10) endif 30 continue 8 write(6,888) kkk,mxiter,vmxdif,dmxdif,emxdif deallocate (y,optk1,optk0) open(3,file='jedc2.out',form='formatted') do 45 i = 1,nkpts do 44 j = 1,nmupts write(3,456) xkpts(i),xmupts(j),y2(i,j,1), & y2(i,j,2),y2(i,j,3),y2(i,j,4), & optk12(i,j,1),optk12(i,j,2), & optk02(i,j,1),optk02(i,j,2) 456 format(10f18.10) 44 continue 45 continue close(3) c Solve for demand functions (more important for elastic labor supply case) do 78 j = 1,nm2pts xmu = xm2pts(j) r0(1) = zgood*alpha*xm2pts(j)**(alpha-one)* & ((one-unempg)*hfix)**(one-alpha) w0(1) = zgood*(one-alpha)*xm2pts(j)**alpha* & ((one-unempg)*hfix)**(-alpha) r0(2) = zbad*alpha*xm2pts(j)**(alpha-one)* & ((one-unempb)*hfix)**(one-alpha) w0(2) = zbad*(one-alpha)*xm2pts(j)**alpha* & ((one-unempb)*hfix)**(-alpha) call inter2(xkpts,xmupts,y2,w11,w10,w01,w00, & w11dp,w10dp,w01dp,w00dp,a,b,xmu) do 77 i = 1,nk2pts xkcur = xk2pts(i) neps = 1 nagg = 1 if (xkcur .le. xkcut1) then call calopt2(xkcur,neps,xkp1,fval1,niter, & 1.0D-10,r0,w0,xkhigh,xklow, & grad1,xkpts,xmupts,xmu,nagg) else if (i .eq. 1) then sk = xkcur else sk = optkg1(i-1,j,1) endif call altopt2(xkcur,neps,xkp1,fval1,niter, & 1.0D-10,sk,grad1,xkpts,r0,w0,xmupts, & xmu,nagg,ndone) if (ndone .eq. 0) then call calopt2(xkcur,neps,xkp1,fval1,niter, & 1.0D-10,r0,w0,xkhigh,xklow, & grad1,xkpts,xmupts,xmu,nagg) endif endif optkg1(i,j,1) = xkp1 vgrid(i,j,1) = fval1 neps = 0 nagg = 1 if (xkcur .le. xkcut2) then call calopt2(xkcur,neps,xkp2,fval2,niter, & 1.0D-10,r0,w0,xkhigh,xklow, & grad0,xkpts,xmupts,xmu,nagg) else if (i .eq. 1) then sk = xkcur else sk = optkg0(i-1,j,1) endif call altopt2(xkcur,neps,xkp2,fval2,niter, & 1.0D-10,sk,grad0,xkpts,r0,w0,xmupts, & xmu,nagg,ndone) if (ndone .eq. 0) then call calopt2(xkcur,neps,xkp2,fval2,niter, & 1.0D-10,r0,w0,xkhigh,xklow, & grad0,xkpts,xmupts,xmu,nagg) endif endif optkg0(i,j,1) = xkp2 vgrid(i,j,2) = fval2 neps = 1 nagg = 2 if (xkcur .le. xkcut1) then call calopt2(xkcur,neps,xkp3,fval3,niter, & 1.0D-10,r0,w0,xkhigh,xklow, & grad1,xkpts,xmupts,xmu,nagg) else if (i .eq. 1) then sk = xkcur else sk = optkg1(i-1,j,2) endif call altopt2(xkcur,neps,xkp3,fval3,niter, & 1.0D-10,sk,grad1,xkpts,r0,w0,xmupts, & xmu,nagg,ndone) if (ndone .eq. 0) then call calopt2(xkcur,neps,xkp3,fval3,niter, & 1.0D-10,r0,w0,xkhigh,xklow, & grad1,xkpts,xmupts,xmu,nagg) endif endif optkg1(i,j,2) = xkp3 vgrid(i,j,3) = fval3 neps = 0 nagg = 2 if (xkcur .le. xkcut2) then call calopt2(xkcur,neps,xkp4,fval4,niter, & 1.0D-10,r0,w0,xkhigh,xklow, & grad0,xkpts,xmupts,xmu,nagg) else if (i .eq. 1) then sk = xkcur else sk = optkg0(i-1,j,2) endif call altopt2(xkcur,neps,xkp4,fval4,niter, & 1.0D-10,sk,grad0,xkpts,r0,w0,xmupts, & xmu,nagg,ndone) if (ndone .eq. 0) then call calopt2(xkcur,neps,xkp4,fval4,niter, & 1.0D-10,r0,w0,xkhigh,xklow, & grad0,xkpts,xmupts,xmu,nagg) endif endif optkg0(i,j,2) = xkp4 vgrid(i,j,4) = fval4 77 continue 78 continue open(5,file='jedc2.vgr',form='formatted') do 71 i = 1,nk2pts do 70 j = 1,nm2pts write(5,603) xk2pts(i),xm2pts(j), & (vgrid(i,j,k),k=1,4) 603 format(6f18.10) 70 continue 71 continue close(5) open(9,file='jedc2.grd',form='formatted') do 99 i = 1,nk2pts do 98 j = 1,nm2pts write(9,400) xk2pts(i),xm2pts(j),optkg1(i,j,1), & optkg1(i,j,2),optkg0(i,j,1), & optkg0(i,j,2) 400 format(6f18.10) 98 continue 99 continue close(9) write(6,7984) 7984 format('computing simulation') call calcf(optkg1,optkg0,xk2pts,xm2pts,a,b,anew,bnew, & xkmean,sdmean,skmean,s4mean,rsqa,rsqb,sderra,sderrb, & gmean,p1mean,p5mean,p10mean,p20mean,p30mean,hansen, & erfmx1,erfmx2) write(6,799) 799 format('moments of ergodic distribution') write(6,949) xkmean,sdmean,skmean,s4mean,gmean write(6,949) p1mean,p5mean,p10mean,p20mean,p30mean write(6,949) hansen,erfmx1/xkmean,erfmx2/xkmean 744 format(f18.10) 949 format(5f15.7) write(6,544) 544 format('measures of fit') write(6,949) rsqa,rsqb,sderra,sderrb open(2,file='jedc2.tmp2',form='formatted') write(2,949) rsqa,rsqb,sderra,sderrb close(2) return end C------------------------------------------------------------------- subroutine splint(xa,ya,y2a,x,y,yp,ydp,khi,klo,kvals,n, & ndoy,ndoyp,ndoydp,nmu,nagg) implicit real*8 (a-h,o-z) parameter (nmupts=4,m=nmupts) dimension xa(n),ya(n,m,2),y2a(n,m,2) one = 1.0D+00 if (kvals .eq. 0) then klo = 1 khi = n 1 if ((khi-klo) .gt. 1) then k = (khi+klo)/2 if (xa(k) .gt. x) then khi = k else klo = k endif goto 1 endif endif h = xa(khi)-xa(klo) a = (xa(khi)-x)/h b = (x-xa(klo))/h asq = a*a bsq = b*b if (ndoy .eq. 1) & y = a*ya(klo,nmu,nagg)+b*ya(khi,nmu,nagg)+ & ((asq*a-a)*y2a(klo,nmu,nagg)+ & (bsq*b-b)*y2a(khi,nmu,nagg))*(h*h)/6.0D+00 if (ndoyp .eq. 1) & yp = (ya(khi,nmu,nagg)-ya(klo,nmu,nagg))/h- & (3.0D+00*asq-one)/6.0D+00*h*y2a(klo,nmu,nagg)+ & (3.0D+00*bsq-one)/6.0D+00*h*y2a(khi,nmu,nagg) if (ndoydp .eq. 1) & ydp = a*y2a(klo,nmu,nagg)+b*y2a(khi,nmu,nagg) return end subroutine splin2(xa,ya,y2a,x,y,yp,ydp,khi,klo,kvals,n, & ndoy,ndoyp,ndoydp,nagg) implicit real*8 (a-h,o-z) parameter (nmupts=4,m=nmupts) dimension xa(n),ya(n,2),y2a(n,2) one = 1.0D+00 if (kvals .eq. 0) then klo = 1 khi = n 1 if ((khi-klo) .gt. 1) then k = (khi+klo)/2 if (xa(k) .gt. x) then khi = k else klo = k endif goto 1 endif endif h = xa(khi)-xa(klo) a = (xa(khi)-x)/h b = (x-xa(klo))/h asq = a*a bsq = b*b if (ndoy .eq. 1) & y = a*ya(klo,nagg)+b*ya(khi,nagg)+ & ((asq*a-a)*y2a(klo,nagg)+ & (bsq*b-b)*y2a(khi,nagg))*(h*h)/6.0D+00 if (ndoyp .eq. 1) & yp = (ya(khi,nagg)-ya(klo,nagg))/h- & (3.0D+00*asq-one)/6.0D+00*h*y2a(klo,nagg)+ & (3.0D+00*bsq-one)/6.0D+00*h*y2a(khi,nagg) if (ndoydp .eq. 1) & ydp = a*y2a(klo,nagg)+b*y2a(khi,nagg) return end c Constructs second derivatives of value functions, evaluated c at K' determined from laws of motion subroutine interp(x1,x2,y,y11,y10,y01,y00,y11dp,y10dp, & y01dp,y00dp,a0,b0) implicit real*8 (a-h,o-z) parameter (nkpts=150,n=nkpts,nmupts=4,m=nmupts) dimension x1(n),y(n,m,4),a(n),b(n),c(n),r(n),u(n),xx(n), & yy(n),y11(n,m,2),y10(n,m,2),y11dp(n,m,2), & y10dp(n,m,2),x2(m),y01(n,m,2),y00(n,m,2), & y01dp(n,m,2),y00dp(n,m,2),a0(2),b0(2), & xmup(m,2),xa(m),ya(m) do 72 i = 1,n xx(i) = x1(i) 72 continue do 70 i = 1,m xmup(i,1) = dexp(a0(1)+a0(2)*dlog(x2(i))) xmup(i,2) = dexp(b0(1)+b0(2)*dlog(x2(i))) 70 continue do 71 i = 1,m xa(i) = x2(i) 71 continue a(1) = 0.0D+00 a(n) = -1.0D+00 b(1) = 1.0D+00 b(n) = 1.05D+00 c(1) = -1.05D+00 c(n) = 0.0D+00 r(1) = 0.0D+00 r(n) = 0.0D+00 do 81 k = 1,4 do 80 ii = 1,m do 79 jj = 1,2 do 30 i = 1,n do 29 j = 1,m ya(j) = y(i,j,k) 29 continue call polint(xa,ya,m,xmup(ii,jj),yval,dy) yy(i) = yval if (k .eq. 1) y11(i,ii,jj) = yy(i) if (k .eq. 2) y10(i,ii,jj) = yy(i) if (k .eq. 3) y01(i,ii,jj) = yy(i) if (k .eq. 4) y00(i,ii,jj) = yy(i) 30 continue do 21 i = 2,n-1 a(i) = (xx(i)-xx(i-1))/6.0D+00 b(i) = (xx(i+1)-xx(i-1))/3.0D+00 c(i) = (xx(i+1)-xx(i))/6.0D+00 r(i) = (yy(i+1)-yy(i))/(xx(i+1)-xx(i))- & (yy(i)-yy(i-1))/(xx(i)-xx(i-1)) 21 continue call tridag(a,b,c,r,u,n) do 22 i = 1,n if (k .eq. 1) y11dp(i,ii,jj) = u(i) if (k .eq. 2) y10dp(i,ii,jj) = u(i) if (k .eq. 3) y01dp(i,ii,jj) = u(i) if (k .eq. 4) y00dp(i,ii,jj) = u(i) 22 continue 79 continue 80 continue 81 continue return end subroutine inter2(x1,x2,y,y11,y10,y01,y00,y11dp,y10dp, & y01dp,y00dp,a0,b0,xmu) implicit real*8 (a-h,o-z) parameter (nkpts=150,n=nkpts,nmupts=4,m=nmupts) dimension x1(n),y(n,m,4),a(n),b(n),c(n),r(n),u(n),xx(n), & yy(n),y11(n,2),y10(n,2),y11dp(n,2),ya(m), & y10dp(n,2),x2(m),y01(n,2),y00(n,2),xa(m), & y01dp(n,2),y00dp(n,2),a0(2),b0(2),xmup(2) do 72 i = 1,n xx(i) = x1(i) 72 continue one = 1.0D+00 xmup(1) = dexp(a0(1)+a0(2)*dlog(xmu)) xmup(2) = dexp(b0(1)+b0(2)*dlog(xmu)) do 71 i = 1,m xa(i) = x2(i) 71 continue a(1) = 0.0D+00 a(n) = -1.0D+00 b(1) = 1.0D+00 b(n) = 1.05D+00 c(1) = -1.05D+00 c(n) = 0.0D+00 r(1) = 0.0D+00 r(n) = 0.0D+00 do 81 k = 1,4 do 79 jj = 1,2 do 30 i = 1,n do 29 j = 1,m ya(j) = y(i,j,k) 29 continue call polint(xa,ya,m,xmup(jj),yval,dy) yy(i) = yval if (k .eq. 1) y11(i,jj) = yy(i) if (k .eq. 2) y10(i,jj) = yy(i) if (k .eq. 3) y01(i,jj) = yy(i) if (k .eq. 4) y00(i,jj) = yy(i) 30 continue do 21 i = 2,n-1 a(i) = (xx(i)-xx(i-1))/6.0D+00 b(i) = (xx(i+1)-xx(i-1))/3.0D+00 c(i) = (xx(i+1)-xx(i))/6.0D+00 r(i) = (yy(i+1)-yy(i))/(xx(i+1)-xx(i))- & (yy(i)-yy(i-1))/(xx(i)-xx(i-1)) 21 continue call tridag(a,b,c,r,u,n) do 22 i = 1,n if (k .eq. 1) y11dp(i,jj) = u(i) if (k .eq. 2) y10dp(i,jj) = u(i) if (k .eq. 3) y01dp(i,jj) = u(i) if (k .eq. 4) y00dp(i,jj) = u(i) 22 continue 79 continue 81 continue return end subroutine polint(xa,ya,n,x,y,dy) implicit real*8 (a-h,o-z) parameter (nmax=5,toler=1.0D-12) dimension xa(n),ya(n),c(nmax),d(nmax) ns = 1 dif = dabs(x-xa(1)) do 11 i = 1,n dift = dabs(x-xa(i)) if (dift .lt. dif) then ns = i dif = dift endif c(i) = ya(i) d(i) = ya(i) 11 continue y = ya(ns) ns = ns-1 do 13 m = 1,n-1 do 12 i = 1,n-m ho = xa(i)-x hp = xa(i+m)-x w = c(i+1)-d(i) den = ho-hp if (dabs(den) .lt. toler) then write(6,100) 100 format('error in polint') pause endif den = w/den d(i) = hp*den c(i) = ho*den 12 continue if ((2*ns) .lt. (n-m)) then dy = c(ns+1) else dy = d(ns) ns = ns-1 endif y = y+dy 13 continue return end subroutine tridag(a,b,c,r,u,n) implicit real*8 (a-h,o-z) parameter (toler=1.0D-12,nmax=5000) dimension gam(nmax),a(n),b(n),c(n),r(n),u(n) bet = b(1) u(1) = r(1)/bet do 11 j = 2,n gam(j) = c(j-1)/bet bet = b(j)-a(j)*gam(j) if (dabs(bet) .lt. toler) then write(6,100) 100 format('failure in tridag') pause endif u(j) = (r(j)-a(j)*u(j-1))/bet 11 continue do 12 j = n-1,1,-1 u(j) = u(j)-gam(j+1)*u(j+1) 12 continue return end c Solves foc using bisection subroutine calopt(xk,neps,ynew,value,niter,toler, & r,w,ykhigh,yklow,grad,xkpts, & xmupts,nmu,nagg) implicit real*8 (a-h,o-z) parameter (mxiter=200,nkpts=150,hfix=0.3271D+00, & delta=0.025D+00,beta=0.99D+00,nmupts=4, & nprint=0,whome=0.0D+00,gamma=1.0D+00) common/cprob/pi(4,4) common/cvs/v11(nkpts,nmupts,2),v11dp(nkpts,nmupts,2), & v10(nkpts,nmupts,2),v10dp(nkpts,nmupts,2), & v01(nkpts,nmupts,2),v01dp(nkpts,nmupts,2), & v00(nkpts,nmupts,2),v00dp(nkpts,nmupts,2) dimension xkpts(nkpts),xmupts(nmupts),r(nmupts,2), & w(nmupts,2) one = 1.0D+00 niter = 0 if (neps .eq. 1) then y = r(nmu,nagg)*xk+w(nmu,nagg)*hfix else y = r(nmu,nagg)*xk+whome endif if (neps .eq. 1 .and. nagg .eq. 1) then npi = 1 elseif (neps .eq. 0 .and. nagg .eq. 1) then npi = 2 elseif (neps .eq. 1 .and. nagg .eq. 2) then npi = 3 elseif (neps .eq. 0 .and. nagg .eq. 2) then npi = 4 endif xklow = yklow xkhigh = ykhigh do 10 i = 1,mxiter xkcur = (xkhigh+xklow)/2.0D+00 x = xkcur-(one-delta)*xk if (x .ge. y) then xkhigh = xkcur goto 10 endif c = y-x upc = -c**(-gamma) call splint(xkpts,v11,v11dp,xkcur,yval,yp1,ydp,khi,klo, & 0,nkpts,0,1,0,nmu,nagg) call splint(xkpts,v10,v10dp,xkcur,yval,yp2,ydp,khi,klo, & 1,nkpts,0,1,0,nmu,nagg) call splint(xkpts,v01,v01dp,xkcur,yval,yp3,ydp,khi,klo, & 1,nkpts,0,1,0,nmu,nagg) call splint(xkpts,v00,v00dp,xkcur,yval,yp4,ydp,khi,klo, & 1,nkpts,0,1,0,nmu,nagg) grad = upc+beta*(pi(npi,1)*yp1+pi(npi,2)*yp2+ & pi(npi,3)*yp3+pi(npi,4)*yp4) if (grad .ge. 0.0D+00) then xklow = xkcur else xkhigh = xkcur endif if (nprint .eq. 1) then write(6,955) xk,xmupts(nmu),neps,nagg,xkcur,grad 955 format(2f12.7,2i3,2f18.14) write(6,757) yp1,yp2,yp3,yp4 757 format(4f18.10) write(6,444) c 444 format(f18.10) pause endif if ((xkhigh-xklow) .le. toler) then niter = i goto 9 endif 10 continue 9 xkcur = (xkhigh+xklow)/2.0D+00 x = xkcur-(one-delta)*xk c = y-x if (c .le. 0.0D+00) then write(6,787) c 787 format('failure in calopt: negative c',f15.10) pause endif if (dabs(gamma-one) .lt. 1.0D-04) then utilc = dlog(c) else utilc = one/(one-gamma)*(c**(one-gamma)-one) endif call splint(xkpts,v11,v11dp,xkcur,yval1,yp,ydp,khi,klo, & 0,nkpts,1,0,0,nmu,nagg) call splint(xkpts,v10,v10dp,xkcur,yval2,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nmu,nagg) call splint(xkpts,v01,v01dp,xkcur,yval3,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nmu,nagg) call splint(xkpts,v00,v00dp,xkcur,yval4,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nmu,nagg) value = utilc+beta*(pi(npi,1)*yval1+pi(npi,2)*yval2+ & pi(npi,3)*yval3+pi(npi,4)*yval4) ynew = xkcur return end subroutine calopt2(xk,neps,ynew,value,niter,toler, & r,w,ykhigh,yklow,grad,xkpts, & xmupts,xmu,nagg) implicit real*8 (a-h,o-z) parameter (mxiter=200,nkpts=150,hfix=0.3271D+00, & delta=0.025D+00,beta=0.99D+00,nmupts=4, & nprint=0,whome=0.0D+00,gamma=1.0D+00) common/cprob/pi(4,4) common/cws/w11(nkpts,2),w11dp(nkpts,2), & w10(nkpts,2),w10dp(nkpts,2), & w01(nkpts,2),w01dp(nkpts,2), & w00(nkpts,2),w00dp(nkpts,2) dimension xkpts(nkpts),xmupts(nmupts),r(2),w(2) one = 1.0D+00 niter = 0 if (neps .eq. 1) then y = r(nagg)*xk+w(nagg)*hfix else y = r(nagg)*xk+whome endif if (neps .eq. 1 .and. nagg .eq. 1) then npi = 1 elseif (neps .eq. 0 .and. nagg .eq. 1) then npi = 2 elseif (neps .eq. 1 .and. nagg .eq. 2) then npi = 3 elseif (neps .eq. 0 .and. nagg .eq. 2) then npi = 4 endif xklow = yklow xkhigh = ykhigh do 10 i = 1,mxiter xkcur = (xkhigh+xklow)/2.0D+00 x = xkcur-(one-delta)*xk if (x .ge. y) then xkhigh = xkcur goto 10 endif c = y-x upc = -c**(-gamma) call splin2(xkpts,w11,w11dp,xkcur,yval,yp1,ydp,khi,klo, & 0,nkpts,0,1,0,nagg) call splin2(xkpts,w10,w10dp,xkcur,yval,yp2,ydp,khi,klo, & 1,nkpts,0,1,0,nagg) call splin2(xkpts,w01,w01dp,xkcur,yval,yp3,ydp,khi,klo, & 1,nkpts,0,1,0,nagg) call splin2(xkpts,w00,w00dp,xkcur,yval,yp4,ydp,khi,klo, & 1,nkpts,0,1,0,nagg) grad = upc+beta*(pi(npi,1)*yp1+pi(npi,2)*yp2+ & pi(npi,3)*yp3+pi(npi,4)*yp4) if (grad .ge. 0.0D+00) then xklow = xkcur else xkhigh = xkcur endif if (nprint .eq. 1) then write(6,955) xk,xmu,neps,nagg,grad,xkcur 955 format(2f12.7,2i3,2f18.14) pause endif if ((xkhigh-xklow) .le. toler) then niter = i goto 9 endif 10 continue 9 xkcur = (xkhigh+xklow)/2.0D+00 x = xkcur-(one-delta)*xk c = y-x if (c .le. 0.0D+00) then write(6,787) c 787 format('failure in calopt: negative c',f15.10) pause endif if (dabs(gamma-one) .lt. 1.0D-04) then utilc = dlog(c) else utilc = one/(one-gamma)*(c**(one-gamma)-one) endif call splin2(xkpts,w11,w11dp,xkcur,yval1,yp,ydp,khi,klo, & 0,nkpts,1,0,0,nagg) call splin2(xkpts,w10,w10dp,xkcur,yval2,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nagg) call splin2(xkpts,w01,w01dp,xkcur,yval3,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nagg) call splin2(xkpts,w00,w00dp,xkcur,yval4,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nagg) value = utilc+beta*(pi(npi,1)*yval1+pi(npi,2)*yval2+ & pi(npi,3)*yval3+pi(npi,4)*yval4) ynew = xkcur return end c Solves foc using Newton-Raphson subroutine altopt(xk,neps,ynew,value,niter,toler, & startk,grad,xkpts,r,w,xmupts,nmu, & nagg,ndone) implicit real*8 (a-h,o-z) parameter (mxiter=200,hfix=0.3271D+00,nmupts=4, & delta=0.025D+00,beta=0.99D+00,nkpts=150, & whome=0.0D+00,nprint=0,gamma=1.0D+00) common/cprob/pi(4,4) common/cvs/v11(nkpts,nmupts,2),v11dp(nkpts,nmupts,2), & v10(nkpts,nmupts,2),v10dp(nkpts,nmupts,2), & v01(nkpts,nmupts,2),v01dp(nkpts,nmupts,2), & v00(nkpts,nmupts,2),v00dp(nkpts,nmupts,2) dimension xkpts(nkpts),xmupts(nmupts),r(nmupts,2), & w(nmupts,2) one = 1.0D+00 niter = 0 ndone = 1 if (neps .eq. 1) then y = r(nmu,nagg)*xk+w(nmu,nagg)*hfix else y = r(nmu,nagg)*xk+whome endif if (neps .eq. 1 .and. nagg .eq. 1) then npi = 1 elseif (neps .eq. 0 .and. nagg .eq. 1) then npi = 2 elseif (neps .eq. 1 .and. nagg .eq. 2) then npi = 3 elseif (neps .eq. 0 .and. nagg .eq. 2) then npi = 4 endif xk2p = startk do 10 i = 1,mxiter xk2 = xk2p x = xk2-(one-delta)*xk c = y-x upc = -c**(-gamma) udpc = -gamma*c**(-gamma-one) call splint(xkpts,v11,v11dp,xk2,yval,yp1,ydp1,khi,klo, & 0,nkpts,0,1,1,nmu,nagg) call splint(xkpts,v10,v10dp,xk2,yval,yp2,ydp2,khi,klo, & 1,nkpts,0,1,1,nmu,nagg) call splint(xkpts,v01,v01dp,xk2,yval,yp3,ydp3,khi,klo, & 1,nkpts,0,1,1,nmu,nagg) call splint(xkpts,v00,v00dp,xk2,yval,yp4,ydp4,khi,klo, & 1,nkpts,0,1,1,nmu,nagg) grad = upc+beta*(pi(npi,1)*yp1+pi(npi,2)*yp2+ & pi(npi,3)*yp3+pi(npi,4)*yp4) hess = udpc+beta*(pi(npi,1)*ydp1+pi(npi,2)*ydp2+ & pi(npi,3)*ydp3+pi(npi,4)*ydp4) xk2p = xk2-grad/hess if (nprint .eq. 1) then write(6,797) xk,neps,nagg write(6,788) xk2,xk2p,grad,hess 797 format(f12.8,2i3) 788 format(4f18.10) pause endif if (dabs(xk2p-xk2) .le. toler) then niter = i goto 9 endif 10 continue 9 if (c .gt. 0.0D+00 .and. xk2 .ge. 0.0D+00 .and. & dabs(grad) .lt. 1.0D-06) then call splint(xkpts,v11,v11dp,xk2,yval1,yp,ydp,khi,klo, & 0,nkpts,1,0,0,nmu,nagg) call splint(xkpts,v10,v10dp,xk2,yval2,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nmu,nagg) call splint(xkpts,v01,v01dp,xk2,yval3,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nmu,nagg) call splint(xkpts,v00,v00dp,xk2,yval4,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nmu,nagg) ynew = xk2 c = y+(one-delta)*xk-ynew if (dabs(gamma-one) .lt. 1.0D-04) then utilc = dlog(c) else utilc = one/(one-gamma)*(c**(one-gamma)-one) endif value = utilc+beta*(pi(npi,1)*yval1+pi(npi,2)*yval2+ & pi(npi,3)*yval3+pi(npi,4)*yval4) else ndone = 0 endif return end subroutine altopt2(xk,neps,ynew,value,niter,toler, & startk,grad,xkpts,r,w,xmupts,xmu, & nagg,ndone) implicit real*8 (a-h,o-z) parameter (mxiter=200,hfix=0.3271D+00,nmupts=4, & delta=0.025D+00,beta=0.99D+00,nkpts=150, & whome=0.0D+00,gamma=1.0D+00) common/cprob/pi(4,4) common/cws/w11(nkpts,2),w11dp(nkpts,2), & w10(nkpts,2),w10dp(nkpts,2), & w01(nkpts,2),w01dp(nkpts,2), & w00(nkpts,2),w00dp(nkpts,2) dimension xkpts(nkpts),xmupts(nmupts),r(2),w(2) one = 1.0D+00 niter = 0 ndone = 1 if (neps .eq. 1) then y = r(nagg)*xk+w(nagg)*hfix else y = r(nagg)*xk+whome endif if (neps .eq. 1 .and. nagg .eq. 1) then npi = 1 elseif (neps .eq. 0 .and. nagg .eq. 1) then npi = 2 elseif (neps .eq. 1 .and. nagg .eq. 2) then npi = 3 elseif (neps .eq. 0 .and. nagg .eq. 2) then npi = 4 endif xk2p = startk do 10 i = 1,mxiter xk2 = xk2p x = xk2-(one-delta)*xk c = y-x upc = -c**(-gamma) udpc = -gamma*c**(-gamma-one) call splin2(xkpts,w11,w11dp,xk2,yval,yp1,ydp1,khi,klo, & 0,nkpts,0,1,1,nagg) call splin2(xkpts,w10,w10dp,xk2,yval,yp2,ydp2,khi,klo, & 1,nkpts,0,1,1,nagg) call splin2(xkpts,w01,w01dp,xk2,yval,yp3,ydp3,khi,klo, & 1,nkpts,0,1,1,nagg) call splin2(xkpts,w00,w00dp,xk2,yval,yp4,ydp4,khi,klo, & 1,nkpts,0,1,1,nagg) grad = upc+beta*(pi(npi,1)*yp1+pi(npi,2)*yp2+ & pi(npi,3)*yp3+pi(npi,4)*yp4) hess = udpc+beta*(pi(npi,1)*ydp1+pi(npi,2)*ydp2+ & pi(npi,3)*ydp3+pi(npi,4)*ydp4) xk2p = xk2-grad/hess if (dabs(xk2p-xk2) .le. toler) then niter = i goto 9 endif 10 continue 9 if (c .gt. 0.0D+00 .and. xk2 .ge. 0.0D+00 .and. & dabs(grad) .lt. 1.0D-06) then call splin2(xkpts,w11,w11dp,xk2,yval1,yp,ydp,khi,klo, & 0,nkpts,1,0,0,nagg) call splin2(xkpts,w10,w10dp,xk2,yval2,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nagg) call splin2(xkpts,w01,w01dp,xk2,yval3,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nagg) call splin2(xkpts,w00,w00dp,xk2,yval4,yp,ydp,khi,klo, & 1,nkpts,1,0,0,nagg) ynew = xk2 c = y+(one-delta)*xk-ynew if (dabs(gamma-one) .lt. 1.0D-04) then utilc = dlog(c) else utilc = one/(one-gamma)*(c**(one-gamma)-one) endif value = utilc+beta*(pi(npi,1)*yval1+pi(npi,2)*yval2+ & pi(npi,3)*yval3+pi(npi,4)*yval4) else ndone = 0 endif return end subroutine hunt(xx,n,x,jlo) implicit real*8 (a-h,o-z) dimension xx(n) logical ascnd ascnd = xx(n) .ge. xx(1) if (jlo .le. 0 .or. jlo .gt. n) then jlo = 0 jhi = n+1 goto 3 endif inc = 1 if (x .ge. xx(jlo) .eqv. ascnd) then 1 jhi = jlo+inc if (jhi .gt. n) then jhi = n+1 elseif (x .ge. xx(jhi) .eqv. ascnd) then jlo = jhi inc = inc+inc goto 1 endif else jhi = jlo 2 jlo = jhi-inc if (jlo .lt. 1) then jlo = 0 elseif (x .lt. xx(jlo) .eqv. ascnd) then jhi = jlo inc = inc+inc goto 2 endif endif 3 if (jhi-jlo .eq. 1) then if (x .eq. xx(n)) jlo = n-1 if (x .eq. xx(1)) jlo = 1 return endif jm = (jhi+jlo)/2 if (x .ge. xx(jm) .eqv. ascnd) then jlo = jm else jhi = jm endif goto 3 end c Simulation routine subroutine calcf(optkg1,optkg0,xk2pts,xm2pts, & a0,b0,a0new,b0new,xkmean,sdmean, & skmean,s4mean,rsqa,rsqb,sderra,sderrb,gmean, & p1mean,p5mean,p10mean,p20mean,p30mean,hansen, & erfmx1,erfmx2) implicit real*8 (a-h,o-z) parameter (nk2pts=5000,ncapt=10000,ndrop=500,nm2pts=30, & unempb=0.1D+00,unempg=0.04D+00,durub=2.5D+00, & durug=1.5D+00,nread2=1,alpha=0.36D+00, & hfix=0.3271D+00,delta=0.025D+00,zgood=1.01D+00, & zbad=0.99D+00,gamma=1.0D+00,beta=0.99D+00) common/cfs/f1(nk2pts,2) common/cprob/pi(4,4) common/nshock/nz(ncapt+1) common/cvgrid/vgrid(nk2pts,nm2pts,4) common/cgmm/rent(ncapt),cons(ncapt) common/vgmm/vhat(2,2),wght(2,2) external gmm dimension a0(2),b0(2),a0new(2),b0new(2),atinv(28,28), & optkg1(nk2pts,nm2pts,2),abdev(ncapt),temp(1,28), & optkg0(nk2pts,nm2pts,2),fracs(ncapt),btt(1,28), & xk2pts(nk2pts),xm2pts(nm2pts),xdata3(ncapt,28), & xkvals(ncapt+1),f0(nk2pts,2),xdata2(ncapt,28), & xdata(ncapt,2),ydata(ncapt),stdxk(ncapt),ttt(1,1), & skewxk(ncapt),xkurtk(ncapt),ginis(ncapt), & pct1(ncapt),pct5(ncapt),pct10(ncapt),pct20(ncapt), & pct30(ncapt),toppct(5),tops(5),bt(28,1),at(28,28), & vmus(ncapt),xdatan(ncapt,6),ydatan(ncapt),bhat(6), & xdatas(ncapt,8),p0new(8),q0new(8),p(3,2),y(3) one = 1.0D+00 xkinc = xk2pts(2)-xk2pts(1) xminc = xm2pts(2)-xm2pts(1) toppct(1) = 0.01D+00 toppct(2) = 0.05D+00 toppct(3) = 0.1D+00 toppct(4) = 0.2D+00 toppct(5) = 0.3D+00 c Individual state transition probabilities pgg00 = (durug-one)/durug pbb00 = (durub-one)/durub pbg00 = 1.25D+00*pbb00 pgb00 = 0.75D+00*pgg00 pgg01 = (unempg-unempg*pgg00)/(one-unempg) pbb01 = (unempb-unempb*pbb00)/(one-unempb) pbg01 = (unempb-unempg*pbg00)/(one-unempg) pgb01 = (unempg-unempb*pgb00)/(one-unempb) pgg10 = one-(durug-one)/durug pbb10 = one-(durub-one)/durub pbg10 = one-1.25D+00*pbb00 pgb10 = one-0.75D+00*pgg00 pgg11 = one-(unempg-unempg*pgg00)/(one-unempg) pbb11 = one-(unempb-unempb*pbb00)/(one-unempb) pbg11 = one-(unempb-unempg*pbg00)/(one-unempg) pgb11 = one-(unempg-unempb*pgb00)/(one-unempb) xkmean = 0.0D+00 sdmean = 0.0D+00 skmean = 0.0D+00 s4mean = 0.0D+00 gmean = 0.0D+00 p1mean = 0.0D+00 p5mean = 0.0D+00 p10mean = 0.0D+00 p20mean = 0.0D+00 p30mean = 0.0D+00 c Initial distribution open(2,file='saving.f1',form='formatted') do 66 i = 1,nk2pts read(2,547) dum,f1(i,1),f1(i,2) 547 format(3f18.10) 66 continue close(2) fsum = 0.0D+00 do 38 i = 1,nk2pts fsum = fsum+f1(i,1)+f1(i,2) 38 continue if (dabs(fsum-one) .gt. 1.0D-04) then write(6,575) fsum 575 format('initial distribution integrates to :',f18.10) pause endif xkbcur = 0.0D+00 do 45 i = 1,nk2pts xkbcur = xkbcur+xk2pts(i)*(f1(i,1)+f1(i,2)) 45 continue c Simulation loop open(7,file='jedc2.sim',form='formatted') do 12 ii = 1,ncapt nzcur = nz(ii) if (nzcur .eq. 1) then nz2 = 1 nz3 = 2 else nz2 = 3 nz3 = 4 endif if (nzcur .eq. 1) then if (nz(ii+1) .eq. 1) then prob0 = pgg11 prob1 = pgg10 else prob0 = pbg11 prob1 = pbg10 endif else if (nz(ii+1) .eq. 1) then prob0 = pgb11 prob1 = pgb10 else prob0 = pbb11 prob1 = pbb10 endif endif xkvals(ii) = xkbcur varcur = 0.0D+00 skwcur = 0.0D+00 xkurt = 0.0D+00 unemp = 0.0D+00 frpoor = 0.0D+00 absval = 0.0D+00 do 82 i = 1,nk2pts varcur = varcur+(f1(i,1)+f1(i,2))* & ((xk2pts(i)-xkbcur)*(xk2pts(i)-xkbcur)) skwcur = skwcur+(f1(i,1)+f1(i,2))* & ((xk2pts(i)-xkbcur)*(xk2pts(i)-xkbcur)* & (xk2pts(i)-xkbcur)) xkurt = xkurt+(f1(i,1)+f1(i,2))* & ((xk2pts(i)-xkbcur)*(xk2pts(i)-xkbcur)* & (xk2pts(i)-xkbcur)*(xk2pts(i)-xkbcur)) unemp = unemp+f1(i,2) absval = absval+dabs(xk2pts(i)-xkbcur)*(f1(i,1)+f1(i,2)) if (xk2pts(i) .lt. 0.1D+00*xkbcur) frpoor = & frpoor+f1(i,1)+f1(i,2) 82 continue stdxk(ii) = dsqrt(varcur) skewxk(ii) = skwcur/(stdxk(ii)*varcur) xkurtk(ii) = xkurt/(varcur*varcur) fracs(ii) = frpoor abdev(ii) = absval call getgini(f1,xk2pts,xkbcur,ginis(ii),toppct,tops) pct1(ii) = tops(1) pct5(ii) = tops(2) pct10(ii) = tops(3) pct20(ii) = tops(4) pct30(ii) = tops(5) if (nz(ii) .eq. 1) then rent(ii) = alpha*zgood*xkvals(ii)**(alpha-one)* & (hfix*(one-unempg))**(one-alpha)-delta wage = (one-alpha)*zgood*xkvals(ii)**alpha* & (hfix*(one-unempg))**(-alpha) gdp = zgood*xkvals(ii)**alpha*(hfix*(one-unempg))** & (one-alpha) else rent(ii) = alpha*zbad*xkvals(ii)**(alpha-one)* & (hfix*(one-unempb))**(one-alpha)-delta wage = (one-alpha)*zbad*xkvals(ii)**alpha* & (hfix*(one-unempb))**(-alpha) gdp = zbad*xkvals(ii)**alpha*(hfix*(one-unempb))** & (one-alpha) endif f0 = f1 f1 = 0.0D+00 c Determine location of current K in K grid call hunt(xm2pts,nm2pts,xkbcur,mcor) if (mcor .lt. 1) then capa = 0.0D+00 mcor = 1 elseif (mcor .ge. nm2pts) then capa = 1.0D+00 mcor = nm2pts-1 else capa = (xkbcur-xm2pts(mcor))/xminc endif c Use decision rules to update distribution vmu = 0.0D+00 do 88 i = 1,nk2pts xk1 = optkg1(i,mcor,nzcur)+capa* & (optkg1(i,mcor+1,nzcur)-optkg1(i,mcor,nzcur)) xk0 = optkg0(i,mcor,1)+capa* & (optkg0(i,mcor+1,nzcur)-optkg0(i,mcor,nzcur)) v1 = vgrid(i,mcor,nz2)+capa* & (vgrid(i,mcor+1,nz2)-vgrid(i,mcor,nz2)) v0 = vgrid(i,mcor,nz3)+capa* & (vgrid(i,mcor+1,nz3)-vgrid(i,mcor,nz3)) call hunt(xk2pts,nk2pts,xk1,jlo1) if (jlo1 .lt. 1) then wght1 = 1.0D+00 jlo1 = 1 elseif (jlo1 .ge. nk2pts) then wght1 = 0.0D+00 jlo1 = nk2pts-1 else wght1 = one-(xk1-xk2pts(jlo1))/xkinc endif call hunt(xk2pts,nk2pts,xk0,jlo0) if (jlo0 .lt. 1) then wght0 = 1.0D+00 jlo0 = 1 elseif (jlo0 .ge. nk2pts) then wght0 = 0.0D+00 jlo0 = nk2pts-1 else wght0 = one-(xk0-xk2pts(jlo0))/xkinc endif f1(jlo1,1) = prob0*wght1*f0(i,1)+f1(jlo1,1) f1(jlo1+1,1) = prob0*(one-wght1)*f0(i,1)+ & f1(jlo1+1,1) f1(jlo0,1) = prob1*wght0*f0(i,2)+f1(jlo0,1) f1(jlo0+1,1) = prob1*(one-wght0)*f0(i,2)+ & f1(jlo0+1,1) f1(jlo1,2) = (one-prob0)*wght1*f0(i,1)+f1(jlo1,2) f1(jlo1+1,2) = (one-prob0)*(one-wght1)* & f0(i,1)+f1(jlo1+1,2) f1(jlo0,2) = (one-prob1)*wght0*f0(i,2)+f1(jlo0,2) f1(jlo0+1,2) = (one-prob1)*(one-wght0)* & f0(i,2)+f1(jlo0+1,2) vmu = vmu+v1*f0(i,1)+v0*f0(i,2) 88 continue vmus(ii) = vmu c Compute K' from new distribution f1 xkbcur = 0.0D+00 fsum = 0.0D+00 do 83 i = 1,nk2pts xkbcur = xkbcur+xk2pts(i)*(f1(i,1)+f1(i,2)) fsum = fsum+f1(i,1)+f1(i,2) 83 continue if (dabs(fsum-one) .gt. 1.0D-04) then write(6,783) ii,fsum 783 format('does not integrate to one in period ',i6,f18.10) pause endif xinv = xkbcur-(one-delta)*xkvals(ii) cons(ii) = gdp-xinv write(7,504) ii,nz(ii),xkvals(ii), & stdxk(ii),skewxk(ii),xkurtk(ii),unemp,fracs(ii), & abdev(ii),rent(ii),wage,ginis(ii),gdp,cons(ii), & xinv,pct1(ii),pct5(ii),pct10(ii),pct20(ii), & pct30(ii),vmus(ii) 504 format(i6,i3,19f15.9) if (ii .gt. ndrop) then xkmean = xkmean+xkvals(ii)/dble(real(ncapt-ndrop-1)) sdmean = sdmean+stdxk(ii)/dble(real(ncapt-ndrop-1)) skmean = skmean+skewxk(ii)/dble(real(ncapt-ndrop-1)) s4mean = s4mean+xkurtk(ii)/dble(real(ncapt-ndrop-1)) gmean = gmean+ginis(ii)/dble(real(ncapt-ndrop-1)) p1mean = p1mean+pct1(ii)/dble(real(ncapt-ndrop-1)) p5mean = p5mean+pct5(ii)/dble(real(ncapt-ndrop-1)) p10mean = p10mean+pct10(ii)/dble(real(ncapt-ndrop-1)) p20mean = p20mean+pct20(ii)/dble(real(ncapt-ndrop-1)) p30mean = p30mean+pct30(ii)/dble(real(ncapt-ndrop-1)) endif 12 continue close(7) c OLS regressions to determine implied law of motion npts = 0 do 34 i = ndrop+1,ncapt-1 if (nz(i) .eq. 1) then npts = npts+1 xdata(npts,1) = one xdata(npts,2) = dlog(xkvals(i)) ydata(npts) = dlog(xkvals(i+1)) endif 34 continue call ols2(ydata,xdata,ncapt,a0new,rsqa,sderra,npts) npts = 0 do 33 i = ndrop+1,ncapt-1 if (nz(i) .eq. 2) then npts = npts+1 xdata(npts,1) = one xdata(npts,2) = dlog(xkvals(i)) ydata(npts) = dlog(xkvals(i+1)) endif 33 continue call ols2(ydata,xdata,ncapt,b0new,rsqb,sderrb,npts) c den Haan-Marcet statistic with chi-square(28) distribution npts = 0 do 39 i = ndrop+1,ncapt-1 if (nz(i) .eq. 1) then npts = npts+1 xdata2(npts,1) = one xdata2(npts,2) = dlog(xkvals(i)) xdata2(npts,3) = xdata2(npts,2)*xdata2(npts,2) xdata2(npts,4) = dlog(stdxk(i)) xdata2(npts,5) = xdata2(npts,4)*xdata2(npts,4) xdata2(npts,6) = xdata2(npts,2)*xdata2(npts,4) xdata2(npts,7) = skewxk(i) xdata2(npts,8) = xdata2(npts,7)*xdata2(npts,7) xdata2(npts,9) = dlog(xkurtk(i)) xdata2(npts,10) = xdata2(npts,9)*xdata2(npts,9) xdata2(npts,11) = xdata2(npts,2)*xdata2(npts,7) xdata2(npts,12) = xdata2(npts,2)*xdata2(npts,9) xdata2(npts,13) = xdata2(npts,4)*xdata2(npts,7) xdata2(npts,14) = xdata2(npts,4)*xdata2(npts,9) xdata2(npts,15) = xdata2(npts,7)*xdata2(npts,9) xdata2(npts,16) = dlog(fracs(i)) xdata2(npts,17) = xdata2(npts,16)*xdata2(npts,16) xdata2(npts,18) = xdata2(npts,16)*xdata2(npts,2) xdata2(npts,19) = xdata2(npts,16)*xdata2(npts,4) xdata2(npts,20) = xdata2(npts,16)*xdata2(npts,7) xdata2(npts,21) = xdata2(npts,16)*xdata2(npts,9) xdata2(npts,22) = dlog(abdev(i)) xdata2(npts,23) = xdata2(npts,22)*xdata2(npts,22) xdata2(npts,24) = xdata2(npts,22)*xdata2(npts,2) xdata2(npts,25) = xdata2(npts,22)*xdata2(npts,4) xdata2(npts,26) = xdata2(npts,22)*xdata2(npts,7) xdata2(npts,27) = xdata2(npts,22)*xdata2(npts,9) xdata2(npts,28) = xdata2(npts,22)*xdata2(npts,16) endif 39 continue npts = 0 do 98 i = ndrop+1,ncapt-1 if (nz(i) .eq. 2) then npts = npts+1 xdata3(npts,1) = one xdata3(npts,2) = dlog(xkvals(i)) xdata3(npts,3) = xdata3(npts,2)*xdata3(npts,2) xdata3(npts,4) = dlog(stdxk(i)) xdata3(npts,5) = xdata3(npts,4)*xdata3(npts,4) xdata3(npts,6) = xdata3(npts,2)*xdata3(npts,4) xdata3(npts,7) = skewxk(i) xdata3(npts,8) = xdata3(npts,7)*xdata3(npts,7) xdata3(npts,9) = dlog(xkurtk(i)) xdata3(npts,10) = xdata3(npts,9)*xdata3(npts,9) xdata3(npts,11) = xdata3(npts,2)*xdata3(npts,7) xdata3(npts,12) = xdata3(npts,2)*xdata3(npts,9) xdata3(npts,13) = xdata3(npts,4)*xdata3(npts,7) xdata3(npts,14) = xdata3(npts,4)*xdata3(npts,9) xdata3(npts,15) = xdata3(npts,7)*xdata3(npts,9) xdata3(npts,16) = dlog(fracs(i)) xdata3(npts,17) = xdata3(npts,16)*xdata3(npts,16) xdata3(npts,18) = xdata3(npts,16)*xdata3(npts,2) xdata3(npts,19) = xdata3(npts,16)*xdata3(npts,4) xdata3(npts,20) = xdata3(npts,16)*xdata3(npts,7) xdata3(npts,21) = xdata3(npts,16)*xdata3(npts,9) xdata3(npts,22) = dlog(abdev(i)) xdata3(npts,23) = xdata3(npts,22)*xdata3(npts,22) xdata3(npts,24) = xdata3(npts,22)*xdata3(npts,2) xdata3(npts,25) = xdata3(npts,22)*xdata3(npts,4) xdata3(npts,26) = xdata3(npts,22)*xdata3(npts,7) xdata3(npts,27) = xdata3(npts,22)*xdata3(npts,9) xdata3(npts,28) = xdata3(npts,22)*xdata3(npts,16) endif 98 continue nz1 = 0 nz0 = 0 bt = 0.0D+00 at = 0.0D+00 capt = dble(real(ncapt-ndrop-1)) erfmx1 = 0.0D+00 erfmx2 = 0.0D+00 do 78 i = ndrop,ncapt-1 if (nz(i) .eq. 1) then nz1 = nz1+1 error = dlog(xkvals(i+1))-a0new(1)- & a0new(2)*dlog(xkvals(i)) do 77 j = 1,28 bt(j,1) = bt(j,1)+error*xdata2(nz1,j)/capt btt(1,j) = bt(j,1) do 76 k = 1,28 at(j,k) = at(j,k)+error*error*xdata2(nz1,j)* & xdata2(nz1,k)/capt 76 continue 77 continue if (dabs(error) .ge. erfmx1) erfmx1 = dabs(error) else nz0 = nz0+1 error = dlog(xkvals(i+1))-b0new(1)- & b0new(2)*dlog(xkvals(i)) do 79 j = 1,28 bt(j,1) = bt(j,1)+error*xdata3(nz0,j)/capt btt(1,j) = bt(j,1) do 80 k = 1,28 at(j,k) = at(j,k)+error*error*xdata3(nz0,j)* & xdata3(nz0,k)/capt 80 continue 79 continue if (dabs(error) .ge. erfmx2) erfmx2 = dabs(error) endif 78 continue call inv28(at,atinv,nsing) temp = matmul(btt,atinv) ttt = matmul(temp,bt) hansen = ttt(1,1)*capt c Naive robustness test npts = 0 do 70 i = ndrop,ncapt-1 npts = npts+1 if (nz(i) .eq. 1) then error = dlog(xkvals(i+1))-a0new(1)-a0new(2)*dlog(xkvals(i)) else error = dlog(xkvals(i+1))-b0new(1)-b0new(2)*dlog(xkvals(i)) endif xdatan(npts,1) = dlog(abdev(i)) xdatan(npts,2) = dlog(stdxk(i)) xdatan(npts,3) = dlog(skewxk(i)) xdatan(npts,4) = dlog(xkurtk(i)) xdatan(npts,5) = dlog(ginis(i)) xdatan(npts,6) = dlog(fracs(i)) ydatan(npts) = error 70 continue call ols6(ydatan,xdatan,ncapt,bhat,rsqe,sderre,npts) write(6,289) 289 format('naive test statistics') write(6,729) (bhat(i),i=1,6) 729 format(3f18.10) write(6,729) rsqe,sderre npts = 0 do 354 i = ndrop+1,ncapt-1 if (nz(i) .eq. 1) then npts = npts+1 xdatas(npts,1) = one xdatas(npts,2) = dlog(xkvals(i)) xdatas(npts,3) = dlog(stdxk(i)) xdatas(npts,4) = dlog(skewxk(i)) xdatas(npts,5) = dlog(xkurtk(i)) xdatas(npts,6) = dlog(ginis(i)) xdatas(npts,7) = dlog(fracs(i)) xdatas(npts,8) = dlog(abdev(i)) ydata(npts) = dlog(stdxk(i+1)) endif 354 continue call ols8(ydata,xdatas,ncapt,p0new,rsqp,sderrp,npts) npts = 0 do 353 i = ndrop+1,ncapt-1 if (nz(i) .eq. 2) then npts = npts+1 xdatas(npts,1) = one xdatas(npts,2) = dlog(xkvals(i)) xdatas(npts,3) = dlog(stdxk(i)) xdatas(npts,4) = dlog(skewxk(i)) xdatas(npts,5) = dlog(xkurtk(i)) xdatas(npts,6) = dlog(ginis(i)) xdatas(npts,7) = dlog(fracs(i)) xdatas(npts,8) = dlog(abdev(i)) ydata(npts) = dlog(stdxk(i+1)) endif 353 continue call ols8(ydata,xdatas,ncapt,q0new,rsqq,sderrq,npts) write(6,299) 299 format('naive test statistics for Sigma') write(6,739) (p0new(i),i=1,4) 739 format(4f18.10) write(6,739) (p0new(i),i=5,8) write(6,729) rsqp,sderrp write(6,739) (q0new(i),i=1,4) write(6,739) (q0new(i),i=5,8) write(6,729) rsqq,sderrq c GMM estimation of preference parameters from aggregate consumption and return data vhat(1,1) = one vhat(1,2) = 0.0D+00 vhat(2,1) = 0.0D+00 vhat(2,2) = one call inv2(vhat,wght,nsing) p(1,1) = beta p(1,2) = gamma p(2,1) = 0.8D+00 p(2,2) = gamma p(3,1) = beta p(3,2) = 4.0D+00 do 67 i = 1,3 y(i) = gmm(p(i,:)) 67 continue call amoeba(p,y,3,2,2,1.0D-08,gmm,niter) bhat2 = p(1,1) ghat2 = p(1,2) vhat = 0.0D+00 do 56 i = ndrop+1,ncapt-1 err1 = one-bhat2*(cons(i+1)/cons(i))**(-ghat2)*(one+rent(i+1)) err2 = err1*rent(i) vhat(1,1) = vhat(1,1)+err1*err1 vhat(1,2) = vhat(1,2)+err1*err2 vhat(2,1) = vhat(2,1)+err2*err1 vhat(2,2) = vhat(2,2)+err2*err2 56 continue call inv2(vhat,wght,nsing) p(1,1) = bhat2 p(1,2) = ghat2 p(2,1) = 0.8D+00 p(2,2) = ghat2 p(3,1) = bhat2 p(3,2) = 4.0D+00 do 60 i = 1,3 y(i) = gmm(p(i,:)) 60 continue call amoeba(p,y,3,2,2,1.0D-08,gmm,niter) bhat2 = p(1,1) ghat2 = p(1,2) write(6,974) 974 format('GMM estimates') write(6,737) beta,bhat2,gamma,ghat2 write(6,737) y(1) 737 format(4f18.10) return end c Computes gini coefficient for current distribution f1 subroutine getgini(f1,xk2pts,xkmean,gini,toppct,tops) implicit real*8 (a-h,o-z) parameter (nk2pts=5000) dimension xpt(nk2pts),ypt(nk2pts),aa(nk2pts),bb(nk2pts), & cc(nk2pts),rr(nk2pts),ydp(nk2pts),xk2pts(nk2pts), & f1(nk2pts,2),tops(5),toppct(5) one = 1.0D+00 sumx = 0.0D+00 sumy = 0.0D+00 xpt(1) = 0.0D+00 ypt(1) = 0.0D+00 nhpts = 1 do 89 j = 2,nk2pts sumx = sumx+f1(j-1,1)+f1(j-1,2) sumy = sumy+(f1(j-1,1)+f1(j-1,2))*xk2pts(j-1)/ & xkmean if (f1(j-1,1)+f1(j-1,2) .gt. 1.0D-10) then nhpts = nhpts+1 xpt(nhpts) = sumx ypt(nhpts) = sumy endif 89 continue call calcsd(xpt,ypt,ydp,aa,bb,cc,rr,nhpts,-one,-one) call integr(xpt,ypt,ydp,nhpts,xint) gini = one-2.0D+00*xint do 10 i = 1,5 call spling(xpt,ypt,ydp,one-toppct(i),yval,yp,ydp,khi,klo, & 0,nhpts,1,0,0) tops(i) = one-yval 10 continue return end subroutine integr(xa,ya,y2a,n,xint) implicit real*8 (a-h,o-z) dimension xa(n),ya(n),y2a(n) one = 1.0D+00 one24 = one/24.0D+00 sum = 0.0D+00 do 10 i = 1,n-1 diff = xa(i+1)-xa(i) sum = sum+diff*( (ya(i)+ya(i+1))/2.0D+00 - & one24*diff*diff*(y2a(i)+y2a(i+1))) 10 continue xint = sum return end subroutine calcsd(x,y,ydp,a,b,c,r,n,an,c1) implicit real*8 (a-h,o-z) dimension x(n),y(n),ydp(n),a(n),b(n),c(n),r(n) a(1) = 0.0D+00 a(n) = an b(1) = 1.0D+00 b(n) = 1.0D+00 c(1) = c1 c(n) = 0.0D+00 r(1) = 0.0D+00 r(n) = 0.0D+00 do 21 i = 2,n-1 a(i) = (x(i)-x(i-1))/6.0D+00 b(i) = (x(i+1)-x(i-1))/3.0D+00 c(i) = (x(i+1)-x(i))/6.0D+00 r(i) = (y(i+1)-y(i))/(x(i+1)-x(i))- & (y(i)-y(i-1))/(x(i)-x(i-1)) 21 continue call tridag(a,b,c,r,ydp,n) return end subroutine spling(xa,ya,y2a,x,y,yp,ydp,khi,klo,kvals,n, & ndoy,ndoyp,ndoydp) implicit real*8 (a-h,o-z) dimension xa(n),ya(n),y2a(n) one = 1.0D+00 if (kvals .eq. 0) then klo = 1 khi = n 1 if ((khi-klo) .gt. 1) then k = (khi+klo)/2 if (xa(k) .gt. x) then khi = k else klo = k endif goto 1 endif endif h = xa(khi)-xa(klo) a = (xa(khi)-x)/h b = (x-xa(klo))/h asq = a*a bsq = b*b if (ndoy .eq. 1) & y = a*ya(klo)+b*ya(khi)+ & ((asq*a-a)*y2a(klo)+ & (bsq*b-b)*y2a(khi))*(h*h)/6.0D+00 if (ndoyp .eq. 1) & yp = (ya(khi)-ya(klo))/h- & (3.0D+00*asq-one)/6.0D+00*h*y2a(klo)+ & (3.0D+00*bsq-one)/6.0D+00*h*y2a(khi) if (ndoydp .eq. 1) & ydp = a*y2a(klo)+b*y2a(khi) return end double precision function gmm(theta) implicit real*8 (a-h,o-z) parameter (ncapt=10000,ndrop=500,one=1.0D+00) common/nshock/nz(ncapt+1) common/cgmm/rent(ncapt),cons(ncapt) common/vgmm/vhat(2,2),wght(2,2) dimension theta(2),fval(ncapt,2) beta = theta(1) gamma = theta(2) npts = 0 do 10 i = ndrop+1,ncapt-1 npts = npts+1 fval(npts,1) = one-beta*(cons(i+1)/cons(i))**(-gamma)* & (one+rent(i+1)) fval(npts,2) = fval(npts,1)*rent(i) 10 continue gmm = 0.0D+00 do 11 i = 1,npts gmm = gmm+fval(i,1)*wght(1,1)*fval(i,1)+ & fval(i,1)*wght(1,2)*fval(i,2)+ & fval(i,2)*wght(2,1)*fval(i,1)+ & fval(i,2)*wght(2,2)*fval(i,2) 11 continue gmm = gmm*dble(real(npts)) return end