C     Model solves version of Aiyagari (1994) with search and
C     unemployment insurance.
C     To run, program requires input file xkpts150.in and links to
C     files roots.f and matrix4.f.

      implicit real*8 (a-h,o-z)
      parameter (nk2pts=5000,nkpts=150,nread=0,nread2=0,nread3=0,
     &           toler=1.0D-03,mxiter=20,update=0.25D+00)
      common/cfs/f0(nk2pts,3),f1(nk2pts,3)
      common/cv2/y2(nkpts,3),optk12(nkpts),optk02(nkpts,3),
     &           opta2(nkpts,3)
      common/cgam/gam0,gam1,gam2
      common/cprice/r,tau
      dimension x(4),x2(4),y(4),fv(4),fv1(4),fv2(4),fv3(4),
     &          fv4(4),xjacob(4,4),xjcinv(4,4)

      one = 1.0D+00
      if (nread2 .eq. 1) then
         open(2,file='ui1d.f1',form='formatted')
         do 66 i = 1,nk2pts
            read(2,547) f1(i,1),f1(i,2),f1(i,3)
 547        format(3f18.10)
 66      continue
         close(2)
      else
         do 49 i = 1,nk2pts
            do 48 j = 1,3
               f1(i,j) = 0.0D+00
 48         continue
 49      continue
         f1(275,1) = one/3.0D+00
         f1(275,2) = one/3.0D+00
         f1(275,3) = one/3.0D+00
      endif   

      if (nread .eq. 1) then
         open(2,file='ui1c.in',form='formatted')
         do 98 i = 1,nkpts
            read(2,700) dum,y2(i,1),y2(i,2),y2(i,3),
     &                  optk12(i),optk02(i,1),optk02(i,2),optk02(i,3),
     &                  opta2(i,1),opta2(i,2),opta2(i,3)
 700        format(11f18.10)
 98      continue
         close(2)
      else
         do 97 i = 1,nkpts
            do 96 j = 1,3
               y2(i,j) = 0.0D+00
               opta2(i,j) = 0.0D+00
               optk02(i,j) = 0.0D+00
 96         continue
            optk12(i) = 0.0D+00
 97      continue
      endif
      
c     Newton-Raphson with numerical derivatives to solve for tax rate
c     and gammas that match calibration equations

	if (nread3 .eq. 1) then
         open(2,file='ui1c.out3',form='formatted')
         read(2,777) dum,(x2(i),i=1,4)
         close(2)
	else
         x2(1) = 0.0360812391D+00
         x2(2) = 10.7421647551D+00
         x2(3) = 2.9911138434D+00
         x2(4) = 1.2683391377D+00
      endif

      do 55 ijk = 1,mxiter

         do 45 i = 1,4
            x(i) = x2(i)
            y(i) = x(i)
 45      continue

         call fvec(x,fv)

         y(1) = y(1)+1.0D-07*x(1)
         call fvec(y,fv1)
         y(1) = x(1)

         y(2) = y(2)+1.0D-07*x(2)
         call fvec(y,fv2)
         y(2) = x(2)

         y(3) = y(3)+1.0D-07*x(3)
         call fvec(y,fv3)
         y(3) = x(3)

         y(4) = y(4)+1.0D-07*x(4)
         call fvec(y,fv4)
         y(4) = x(4)

         do 15 i = 1,4
            xjacob(i,1) = (fv1(i)-fv(i))/(1.0D-07*x(1))
            xjacob(i,2) = (fv2(i)-fv(i))/(1.0D-07*x(2))
            xjacob(i,3) = (fv3(i)-fv(i))/(1.0D-07*x(3))
            xjacob(i,4) = (fv4(i)-fv(i))/(1.0D-07*x(4))
 15      continue
         call inv4(xjacob,xjcinv,nsing)

         do 17 i = 1,4
            xjf = 0.0D+00
            do 16 j = 1,4
               xjf = xjf+xjcinv(i,j)*fv(j)
 16         continue
            x2(i) = x(i)-xjf*update
 17      continue

         xmxdif = 0.0D+00
         do 59 i = 1,4
            xdiff = dabs(x2(i)-x(i))
            if (xdiff .gt. xmxdif) xmxdif = xdiff
 59      continue

         write(6,594) (x(i),i=1,4)
         write(6,594) (x2(i),i=1,4)
         write(6,594) (fv(i),i=1,4)
 594     format(4f18.10)

         open(2,file='fvec.out',form='formatted')
         write(2,594) (x(i),i=1,4)
         write(2,594) (x2(i),i=1,4)
         write(2,594) (fv(i),i=1,4)
         close(2)

         if (xmxdif .lt. toler) goto 4

 55   continue

 4    open(2,file='ui1c.out3',form='formatted')
      write(2,777) r,(x2(i),i=1,4)
 777  format(5f18.10)
      close(2)

      end

      subroutine fvec(x,fv)
      implicit real*8 (a-h,o-z)
      parameter (beta=0.99D+00,delta=0.025D+00,one=1.0D+00)
      common/csurpl/surpls,u0,u1,unemp
      dimension x(4),fv(4)
      common/cprice/r,tau
      common/cgam/gam0,gam1,gam2
      external zeta2

      tau = x(1)
      gam0 = x(2)
      gam1 = x(3)
      gam2 = x(4)

c     Brent's method to compute market clearing r

      x1 = 0.03504D+00
      x2 = 0.035085D+00

      r = zbrent(zeta2,x1,x2,1.0D-10)

      fv(1) = surpls
      fv(2) = 0.698D+00-u0
      fv(3) = 0.155D+00-u1
      fv(4) = 0.074D+00-unemp

      return
      end

      double precision function zeta2(x)
      implicit real*8 (a-h,o-z)
      common/csurpl/surpls,u0,u1,unemp
      common/cprice/r,tau
      common/cgam/gam0,gam1,gam2

      r = x

      call dozeta(r,tau,gam0,gam1,gam2,rnew,surpls,u0,u1,u2,
     &            unemp)

      zeta2 = rnew-r

      return
      end

      subroutine dozeta(r,tau,gam0,gam1,gam2,rnew,surpls,u0,u1,u2,
     &                  unemp)
      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,
     &           theta=0.5D+00,ncalop=50,mxiter=2000,
     &           nk2pts=5000,b=0.17D+00)
      common/cv2/y2(nkpts,3),optk12(nkpts),optk02(nkpts,3),
     &           opta2(nkpts,3)
      common/cvs/v0(nkpts),v0dp(nkpts),v1(nkpts),
     &           v1dp(nkpts),v2(nkpts),v2dp(nkpts)
      common/cfs/f0(nk2pts,3),f1(nk2pts,3)      
      dimension y(nkpts,3),ydp(nkpts,3),optag(nk2pts,3),
     &          optk1(nkpts),xkpts(nkpts),optkg1(nk2pts),
     &          opta(nkpts,3),vgrid(nk2pts,3),
     &          optk0(nkpts,3),xk2pts(nk2pts),optkg0(nk2pts,3)

      one = 1.0D+00
      xkcut1 = 0.5D+00
      xkcut2 = 1.5D+00

      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)
      xklow = xkbor
      xkhigh = 1255.0D+00
      xk2pts(1) = xkpts(1)
	xkinc = (xkpts(nkpts)-xkpts(1))/dble(real(nk2pts-1))
	do 93 i = 2,nk2pts
	   xk2pts(i) = xk2pts(i-1)+xkinc
 93   continue

      w = (one-alpha)*((r/alpha)**(alpha/(alpha-one)))

c     value iteration stage

      do 30 kkk = 1,mxiter
                 
         do 12 i = 1,nkpts
            do 11 j = 1,3
               y(i,j) = y2(i,j)
               optk0(i,j) = optk02(i,j)
               opta(i,j) = opta2(i,j)
 11         continue
            optk1(i) = optk12(i)
 12      continue

         call interp(xkpts,y,ydp,v0,v1,v2,v0dp,v1dp,v2dp)

         do 13 i = 1,nkpts
            xkcur = xkpts(i)

            neps = 1
            nhist = 0
            if (xkcur .le. xkcut1 .or. kkk .le. ncalop) then
               call calopt(xkcur,neps,xkp1,fval1,niter,
     &                     1.0D-10,r,w,xkhigh,xklow,grad1,xkpts,
     &                     tau,nhist)
            else
               call altopt(xkcur,neps,xkp1,fval1,niter,
     &                     1.0D-10,optk1(i),grad1,xkpts,r,w,
     &                     tau,nhist)
            endif
            optk12(i) = xkp1

            neps = 0
            nhist = 0
            if (xkcur .le. xkcut2 .or. kkk .le. ncalop) then
               call calopt(xkcur,neps,xkp0,fval00,niter,
     &                     1.0D-10,r,w,xkhigh,xklow,grad0,xkpts,
     &                     tau,nhist)
            else
               call altopt(xkcur,neps,xkp0,fval00,niter,
     &                     1.0D-10,optk0(i,2),grad0,xkpts,r,w,
     &                     tau,nhist)
            endif
            optk02(i,1) = xkp0

            neps = 0
            nhist = 1
            if (xkcur .le. xkcut2 .or. kkk .le. ncalop) then
               call calopt(xkcur,neps,xkp0,fval01,niter,
     &                     1.0D-10,r,w,xkhigh,xklow,grad0,xkpts,
     &                     tau,nhist)
            else
               call altopt(xkcur,neps,xkp0,fval01,niter,
     &                     1.0D-10,optk0(i,1),grad0,xkpts,r,w,
     &                     tau,nhist)
            endif
            optk02(i,2) = xkp0

            neps = 0
            nhist = 2
            if (xkcur .le. xkcut2 .or. kkk .le. ncalop) then
               call calopt(xkcur,neps,xkp0,fval02,niter,
     &                     1.0D-10,r,w,xkhigh,xklow,grad0,xkpts,
     &                     tau,nhist)
            else
               call altopt(xkcur,neps,xkp0,fval02,niter,
     &                     1.0D-10,optk0(i,3),grad0,xkpts,r,w,
     &                     tau,nhist)
            endif
            optk02(i,3) = xkp0

            call altop2(gam0,fval1,fval00,niter,1.0D-10,
     &                  aval1,opta(i,1),grada1,val1)
            opta2(i,1) = aval1
            y2(i,1) = val1
            
            call altop2(gam1,fval1,fval01,niter,1.0D-10,
     &                  aval0,opta(i,2),grada0,val0)
            opta2(i,2) = aval0
            y2(i,2) = val0

            call altop2(gam2,fval1,fval02,niter,1.0D-10,
     &                  aval0,opta(i,3),grada0,val0)
            opta2(i,3) = aval0
            y2(i,3) = val0

 13      continue

         vmxdif = 0.0D+00
         dmxdif = 0.0D+00
         emxdif = 0.0D+00
         amxdif = 0.0D+00
         do 15 i = 1,nkpts
            do 14 j = 1,3
               vdiff = dabs(y2(i,j)-y(i,j))
               if (vdiff .gt. vmxdif) vmxdif = vdiff
               ediff = dabs(optk02(i,j)-optk0(i,j))
               if (ediff .gt. emxdif) emxdif = ediff
               adiff = dabs(opta2(i,j)-opta(i,j))
               if (adiff .gt. amxdif) amxdif = adiff
 14         continue
            ddiff = dabs(optk12(i)-optk1(i))
            if (ddiff .gt. dmxdif) dmxdif = ddiff
 15      continue

         if (vmxdif .lt. 1.0D-06 .and. dmxdif .lt. 1.0D-08
     &       .and. emxdif .lt. 1.0D-08 .and. amxdif .lt. 1.0D-08
     &       .and. kkk .gt. 50) goto 8

         if (multpl(kkk,1000) .eq. 1) then
            write(6,888) kkk,mxiter,vmxdif,dmxdif,emxdif,amxdif
 888        format(2i5,4f15.10)
         endif

 30   continue

 8    write(6,888) kkk,mxiter,vmxdif,dmxdif,emxdif,amxdif

      open(3,file='ui1c.out',form='formatted')
      do 45 i = 1,nkpts
         write(3,456) xkpts(i),y2(i,1),y2(i,2),y2(i,3),optk12(i),
     &                optk02(i,1),optk02(i,2),optk02(i,3),
     &                opta2(i,1),opta2(i,2),opta2(i,3)
 456     format(11f18.10)
 45   continue
      close(3)
      
c     resolve over fine grid for invariant distribution

      do 77 i = 1,nk2pts
         
         xkcur = xk2pts(i)
         
         neps = 1
         nhist = 0
         if (xkcur .le. xkcut1) then
            call calopt(xkcur,neps,xkp1,fval1,niter,1.0D-10,
     &                  r,w,xkhigh,xklow,grad1,xkpts,tau,nhist)
         else
            call altopt(xkcur,neps,xkp1,fval1,niter,1.0D-10,
     &                  xkcur,grad1,xkpts,r,w,tau,nhist)
         endif
         optkg1(i) = xkp1

         neps = 0
         nhist = 0
         if (xkcur .le. xkcut2) then
            call calopt(xkcur,neps,xkp0,fval00,niter,1.0D-10,
     &                  r,w,xkhigh,xklow,grad0,xkpts,tau,nhist)
         else
            call altopt(xkcur,neps,xkp0,fval00,niter,1.0D-10,
     &                  xkcur,grad0,xkpts,r,w,tau,nhist)
         endif
         optkg0(i,1) = xkp0

         neps = 0
         nhist = 1
         if (xkcur .le. xkcut2) then
            call calopt(xkcur,neps,xkp0,fval01,niter,1.0D-10,
     &                  r,w,xkhigh,xklow,grad0,xkpts,tau,nhist)
         else
            call altopt(xkcur,neps,xkp0,fval01,niter,1.0D-10,
     &                  xkcur,grad0,xkpts,r,w,tau,nhist)
         endif
         optkg0(i,2) = xkp0

         neps = 0
         nhist = 2
         if (xkcur .le. xkcut2) then
            call calopt(xkcur,neps,xkp0,fval02,niter,1.0D-10,
     &                  r,w,xkhigh,xklow,grad0,xkpts,tau,nhist)
         else
            call altopt(xkcur,neps,xkp0,fval02,niter,1.0D-10,
     &                  xkcur,grad0,xkpts,r,w,tau,nhist)
         endif
         optkg0(i,3) = xkp0

         call altop2(gam0,fval1,fval00,niter,1.0D-10,aval1,
     &               0.25D+00,grada1,val0)
         optag(i,1) = aval1
         vgrid(i,1) = val0
         call altop2(gam1,fval1,fval01,niter,1.0D-10,aval0,
     &               0.4D+00,grada0,val1)
         optag(i,2) = aval0
         vgrid(i,2) = val1
         call altop2(gam2,fval1,fval02,niter,1.0D-10,aval0,
     &               0.4D+00,grada0,val2)
         optag(i,3) = aval0
         vgrid(i,3) = val2

 77   continue

      open(5,file='vgrid.ui1c',form='formatted')
      do 71 i = 1,nk2pts
         write(5,603) xk2pts(i),(vgrid(i,j),j=1,3)
 603     format(4f18.10)
 71   continue
      close(5)

      open(9,file='ui1cgrd.out',form='formatted')
      do 99 i = 1,nk2pts
         write(9,400) xk2pts(i),optkg1(i),optkg0(i,1),optkg0(i,2),
     &                optkg0(i,3),optag(i,1),
     &                optag(i,2),optag(i,3)
 400     format(8f18.10)
 99   continue
      close(9)

      write(6,7984)
 7984 format('computing invariant distribution')

      call calcf(xkpts,gam0,gam1,gam2,optkg1,optkg0,optag,
     &           xkmean,unemp,u0,u1,u2,xk2pts)
      unemp0 = u0*unemp
      unemp1 = u1*unemp
      unemp2 = u2*unemp

      surpls = tau*w*hfix*(one-unemp)-
     &         (theta*(unemp0+unemp1)+b*unemp2)*w*hfix

      write(6,509) xkmean,unemp
 509  format(2f18.10)

      open(21,file='ui1c.out2',form='formatted')
      write(21,338) xkmean,gam0,gam1,gam2
      write(21,338) rnew,r,tau,surpls
      write(21,338) unemp,u0,u1,u2
 338  format(5f18.10)
      close(21)

      open(22,file='ui1c.f1',form='formatted')
      do 81 i = 1,nk2pts
         write(22,499) f1(i,1),f1(i,2),f1(i,3)
 499     format(3f18.10)
 81      continue
      close(22)

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C   subroutines and functions
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

c     evaluates spline interpolation

      subroutine splint(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

c     computes second derivatives for spline interpolation

      subroutine interp(x1,y,ydp,y1,y2,y3,y1dp,y2dp,y3dp)
      implicit real*8 (a-h,o-z)
      parameter (nkpts=150,n=nkpts)

      dimension x1(n),y(n,3),a(n),b(n),c(n),r(n),u(n),xx(n),
     &          yy(n),ydp(n,3),y1(n),y2(n),y1dp(n),y2dp(n),
     &          y3(n),y3dp(n)

      do 72 i = 1,n
         xx(i) = x1(i)
 72   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,3
         do 30 i = 1,n
            yy(i) = y(i,k)
            if (k .eq. 1) y1(i) = yy(i)
            if (k .eq. 2) y2(i) = yy(i)
            if (k .eq. 3) y3(i) = 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) y1dp(i) = u(i)
            if (k .eq. 2) y2dp(i) = u(i)
            if (k .eq. 3) y3dp(i) = u(i)
 22      continue
 81   continue

      do 56 i = 1,n
         do 55 j = 1,3
            if (j .eq. 1) ydp(i,j) = y1dp(i)
            if (j .eq. 2) ydp(i,j) = y2dp(i)
            if (j .eq. 3) ydp(i,j) = y3dp(i)
 55      continue
 56   continue

      return
      end

      subroutine tridag(a,b,c,r,u,n)
      implicit real*8 (a-h,o-z)
      parameter (toler=1.0D-12,nmax=150)
      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

      integer function multpl(i,m)
      implicit real*8 (a-h,o-z)
      parameter (toler=1.0D-12)

      remain = dble(real(i))/dble(real(m))-i/m

      if (remain .le. toler) then
         multpl = 1
      else
         multpl = 0
      endif

      return
      end

c     solves first-order condition by bisection

      subroutine calopt(xk,neps,ynew,value,niter,toler,
     &                  r,w,ykhigh,yklow,grad,xkpts,tau,nhist)
      implicit real*8 (a-h,o-z)
      parameter (mxiter=200,nkpts=150,hfix=0.3271D+00,b=0.17D+00,
     &           delta=0.025D+00,beta=0.99D+00,theta=0.5D+00)

      common/cvs/v0(nkpts),v0dp(nkpts),
     &           v1(nkpts),v1dp(nkpts),
     &           v2(nkpts),v2dp(nkpts)

      dimension xkpts(nkpts)

      one = 1.0D+00
      niter = 0

      if (neps .eq. 1) then
         y = r*xk+w*hfix*(one-tau)
      elseif (neps .eq. 0 .and. nhist .lt. 2) then
         y = r*xk+theta*w*hfix
      elseif (neps .eq. 0 .and. nhist .eq. 2) then
         y = r*xk+b*w*hfix
      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 = -one/c

         if (neps .eq. 1) then
            call splint(xkpts,v0,v0dp,xkcur,y,yp,ydp,khi,klo,
     &                  0,nkpts,0,1,0)
         elseif (neps .eq. 0 .and. nhist .eq. 0) then
            call splint(xkpts,v1,v1dp,xkcur,y,yp,ydp,khi,klo,
     &                  0,nkpts,0,1,0)
         elseif (neps .eq. 0 .and. nhist .gt. 0) then
            call splint(xkpts,v2,v2dp,xkcur,y,yp,ydp,khi,klo,
     &                  0,nkpts,0,1,0)
         endif
         
         grad = upc+beta*yp

         if (grad .ge. 0.0D+00) then
            xklow = xkcur
         else
            xkhigh = xkcur
         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

      utilc = dlog(c)
      if (neps .eq. 1) then
         call splint(xkpts,v0,v0dp,xkcur,yval,yp,ydp,khi,klo,
     &               0,nkpts,1,0,0)
      elseif (neps .eq. 0 .and. nhist .eq. 0) then
         call splint(xkpts,v1,v1dp,xkcur,yval,yp,ydp,khi,klo,
     &               0,nkpts,1,0,0)
      elseif (neps .eq. 0 .and. nhist .gt. 0) then
         call splint(xkpts,v2,v2dp,xkcur,yval,yp,ydp,khi,klo,
     &               0,nkpts,1,0,0)
      endif

      value = utilc+beta*yval
      ynew = xkcur

      return
      end

c     solves first-order condition by Newton-Raphson

      subroutine altopt(xk,neps,ynew,value,niter,toler,
     &                  startk,grad,xkpts,r,w,tau,nhist)
      implicit real*8 (a-h,o-z)
      parameter (mxiter=200,hfix=0.3271D+00,theta=0.5D+00,
     &           delta=0.025D+00,beta=0.99D+00,nkpts=150,
     &           b=0.17D+00)

      common/cvs/v0(nkpts),v0dp(nkpts),
     &           v1(nkpts),v1dp(nkpts),
     &           v2(nkpts),v2dp(nkpts)

      dimension xkpts(nkpts)

      one = 1.0D+00
      niter = 0

      if (neps .eq. 1) then
         y = r*xk+w*hfix*(one-tau)
      elseif (neps .eq. 0 .and. nhist .lt. 2) then
         y = r*xk+theta*w*hfix
      elseif (neps .eq. 0 .and. nhist .eq. 2) then
         y = r*xk+b*w*hfix
      endif

      xk2p = startk

      do 10 i = 1,mxiter
         xk2 = xk2p
         x = xk2-(one-delta)*xk
         c = y-x

         if (neps .eq. 1) then
            call splint(xkpts,v0,v0dp,xk2,y,yp,ydp,khi,klo,0,
     &                  nkpts,0,1,1)
         elseif (neps .eq. 0 .and. nhist .eq. 0) then
            call splint(xkpts,v1,v1dp,xk2,y,yp,ydp,khi,klo,0,
     &                  nkpts,0,1,1)
         elseif (neps .eq. 0 .and. nhist .gt. 0) then
            call splint(xkpts,v2,v2dp,xk2,y,yp,ydp,khi,klo,0,
     &                  nkpts,0,1,1)
         endif

         upc = -one/c
         udpc = -one/(c*c)

         grad = upc+beta*yp
         hess = udpc+beta*ydp

         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) then
         utilc = dlog(c)

         if (neps .eq. 1) then
            call splint(xkpts,v0,v0dp,xk2,yval,yp,ydp,khi,klo,
     &                  0,nkpts,1,0,0)
         elseif (neps .eq. 0 .and. nhist .eq. 0) then
            call splint(xkpts,v1,v1dp,xk2,yval,yp,ydp,khi,klo,
     &                  0,nkpts,1,0,0)
         elseif (neps .eq. 0 .and. nhist .gt. 0) then
            call splint(xkpts,v2,v2dp,xk2,yval,yp,ydp,khi,klo,
     &                  0,nkpts,1,0,0)
         endif

         ynew = xk2
         value = utilc+beta*yval
      else
         write(6,123) xk,neps,nhist,c
 123     format('failure in altopt: negative c',f8.3,2i3,f12.8)
         pause
      endif

      return
      end

c     finds optimal effort level using Newton-Raphson

      subroutine altop2(gamma,val1,val0,niter,toler,aval,
     &                  starta,grad,value)
      implicit real*8 (a-h,o-z)
      parameter (mxiter=200)

      one = 1.0D+00
      niter = 0

      if (val1 .le. val0) then
         aval = 0.0D+00
         value = val0
         return
      endif

      aprime = starta

      do 10 i = 1,mxiter
         acur = aprime 
         grad = -2.0D+00*acur+gamma*dexp(-gamma*acur)*
     &          (val1-val0)
         hess = -2.0D+00-gamma*gamma*dexp(-gamma*acur)*
     &          (val1-val0)

         aprime = acur-grad/hess

         if (dabs(aprime-acur) .le. toler) then
            niter = i
            goto 9
         endif

 10   continue

 9    aval = aprime
      prob = one-dexp(-gamma*aval)
      value = -aval*aval+prob*val1+(one-prob)*val0

      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     computes invariant distribution using histogram
c     updates current distribution using decision rules over fine grid

      subroutine calcf(xkpts,gam0,gam1,gam2,optkg1,optkg0,optag,
     &           xkmean,unemp,u0,u1,u2,xk2pts)
      implicit real*8 (a-h,o-z)
      parameter (nkpts=150,nk2pts=5000,mxiter=10000000)
      common/cfs/f0(nk2pts,3),f1(nk2pts,3)
      dimension xkpts(nkpts),xk2pts(nk2pts),optkg1(nk2pts),
     &          optkg0(nk2pts,3),optag(nk2pts,3)

      one = 1.0D+00

      do 70 k = 1,mxiter

         do 12 i = 1,nk2pts
            do 11 j = 1,3
               f0(i,j) = f1(i,j)
               f1(i,j) = 0.0D+00
 11         continue
 12      continue

         do 68 ii = 1,nk2pts
            
            prob0 = one-dexp(-gam0*optag(ii,1))
            prob1 = one-dexp(-gam1*optag(ii,2))
            prob2 = one-dexp(-gam2*optag(ii,3))

            call hunt(xk2pts,nk2pts,optkg1(ii),jlo1)
            call hunt(xk2pts,nk2pts,optkg0(ii,1),jlo00)
            call hunt(xk2pts,nk2pts,optkg0(ii,2),jlo01)
            call hunt(xk2pts,nk2pts,optkg0(ii,3),jlo02)

            if (jlo1 .lt. 1) then
               jlo1 = 1
               wght1 = 1.0D+00
            elseif (jlo1 .ge. nk2pts) then
               jlo1 = nk2pts-1
               wght1 = 0.0D+00
            else   
               wght1 = one-(optkg1(ii)-xk2pts(jlo1))/
     &                 (xk2pts(jlo1+1)-xk2pts(jlo1))
            endif

            if (jlo00 .lt. 1) then
               jlo00 = 1
               wght00 = 1.0D+00
            elseif (jlo00 .ge. nk2pts) then
               jlo00 = nk2pts-1
               wght00 = 0.0D+00
            else
               wght00 = one-(optkg0(ii,1)-xk2pts(jlo00))/
     &                  (xk2pts(jlo00+1)-xk2pts(jlo00))
            endif

            if (jlo01 .lt. 1) then
               jlo01 = 1
               wght01 = 1.0D+00
            elseif (jlo01 .ge. nk2pts) then
               jlo01 = nk2pts-1
               wght01 = 0.0D+00
            else
               wght01 = one-(optkg0(ii,2)-xk2pts(jlo01))/
     &                  (xk2pts(jlo01+1)-xk2pts(jlo01))
            endif

            if (jlo02 .lt. 1) then
               jlo02 = 1
               wght02 = 1.0D+00
            elseif (jlo02 .ge. nk2pts) then
               jlo02 = nk2pts-1
               wght02 = 0.0D+00
            else
               wght02 = one-(optkg0(ii,3)-xk2pts(jlo02))/
     &                  (xk2pts(jlo02+1)-xk2pts(jlo02))
            endif

            f1(jlo1,1) = prob0*wght1*f0(ii,1)+f1(jlo1,1)
            f1(jlo1+1,1) = prob0*(one-wght1)*f0(ii,1)+f1(jlo1+1,1)
            f1(jlo1,1) = prob1*wght1*f0(ii,2)+f1(jlo1,1)
            f1(jlo1+1,1) = prob1*(one-wght1)*f0(ii,2)+f1(jlo1+1,1)
            f1(jlo1,1) = prob2*wght1*f0(ii,3)+f1(jlo1,1)
            f1(jlo1+1,1) = prob2*(one-wght1)*f0(ii,3)+f1(jlo1+1,1)

            f1(jlo00,2) = (one-prob0)*wght00*f0(ii,1)+f1(jlo00,2)
            f1(jlo00+1,2) = (one-prob0)*(one-wght00)*f0(ii,1)+
     &                      f1(jlo00+1,2)
            f1(jlo01,3) = (one-prob1)*wght01*f0(ii,2)+f1(jlo01,3)
            f1(jlo01+1,3) = (one-prob1)*(one-wght01)*f0(ii,2)+
     &                      f1(jlo01+1,3)
            f1(jlo02,3) = (one-prob2)*wght02*f0(ii,3)+f1(jlo02,3)
            f1(jlo02+1,3) = (one-prob2)*(one-wght02)*f0(ii,3)+
     &                      f1(jlo02+1,3)

 68      continue

         fmxdif = 0.0D+00
         fsum = 0.0D+00
         do 98 i = 1,nk2pts
            do 97 j = 1,3
               fsum = fsum+f1(i,j)
               fdiff = dabs(f1(i,j)-f0(i,j))
               if (fdiff .gt. fmxdif) fmxdif = fdiff
 97         continue
 98      continue
         if (dabs(fsum-one) .gt. 1.0D-06) then
            write(6,543) k,fsum
 543        format('does not integrate to 1: ',i10,f18.10)
            read *
         endif

         if (multpl(k,10000) .eq. 1) then
            write(6,703) k,mxiter,fmxdif
 703        format(2i10,f18.10)
         endif

         if (fmxdif .le. 1.0D-08) goto 1
 70   continue

 1    sumk = 0.0D+00
      do 75 i = 1,nk2pts
         do 74 j = 1,3
            sumk = sumk+xk2pts(i)*f1(i,j)
 74      continue
 75   continue
      xkmean = sumk
      unemp = 0.0D+00
      u0 = 0.0D+00
      u1 = 0.0D+00
      u2 = 0.0D+00
      do 88 i = 1,nk2pts
         unemp = unemp+dexp(-gam0*optag(i,1))*f1(i,1)+
     &           dexp(-gam1*optag(i,2))*f1(i,2)+
     &           dexp(-gam2*optag(i,3))*f1(i,3)
         u0 = u0+dexp(-gam0*optag(i,1))*f1(i,1)
         u1 = u1+dexp(-gam1*optag(i,2))*f1(i,2)
         u2 = u2+dexp(-gam2*optag(i,3))*f1(i,3)
 88   continue
      u0 = u0/unemp
      u1 = u1/unemp
      u2 = u2/unemp

      return
      end