c***********************************************************************
c       Sample program to generate equilibrium bond yields using
c       calls to ADISET and ADSTEP.  Two factor interest rate model
c       is Jacobs/Jones log model.
c***********************************************************************

      implicit double precision (A-H, K-L, O-Z)
      parameter ( IM =   20, IN =   20 )
      parameter ( IM1= IM+1, IN1= IN+1 )
      dimension   PARM(15), INFO(8)
      dimension   U(0:IM, 0:IN), V(0:IM, 0:IN), ARR(8, 0:IM, 0:IN),
     &            C(0:IM, 0:IN)
      save
           

C  Algorithm parmeters:
      data INFO / 0, 0, 0, 0, 0, 0, 0, 0 /
      data LMN, LMX, RMN, RMX / .00, .30, .00, .30 /
      data IFUT, KK  / 0, .10 /

C  Model params:  SIGR  SIGL  RHO  KAPR  KAPL  LAMR  LAML  LBAR
      data PARM / .72,  .19, -.28, 2.60, .064, -.23,-.008, .085, 7*0. /

      HL = ( LMX - LMN ) / IM
      HR = ( RMX - RMN ) / IN


C  Enter desired time to maturity from console:

 100  print*, ' Enter maturity and coupon interval in yrs (0 to quit):'
      read *, TMAX, DTCOUP
      if ( TMAX .le. 0.0 ) stop

      call COPYC ( 1d0, U, IM1, IN1 )
      call COPYC ( 0d0, C, IM1, IN1 )

C     Adjust KK to be integer fraction of DTCOUP
      KK = DTCOUP / max ( 1d0, dnint ( DTCOUP/KK ) )

      call ADISET(IM,IN,IMU,LMN,LMX,RMN,RMX,KK,IFN,IFUT,INFO,PARM,ARR)

      do 200 T = 0.0D0, TMAX-KK/2, KK
C        Pay coupon if time has come
         if ( mod ( T+KK/2, DTCOUP ) .lt. KK ) then
            call ADDC( 1d0, C, IM1, IN1 )
         endif
C        Step back KK
         call ADSTEP ( U, V, ARR )
         call ADSTEP ( C, V, ARR )
 200  continue

C     Calculate equilibrium coupon rates and store in U()
      do 300 I = 0, IM
         do 300 J = 0, IN
         COUPON = (1.0 - U(I,J)) / ( C(I,J) * DTCOUP )
         U(I,J) = COUPON * 100
 300  continue

C     Print out table of results (10 by 10)
      print*
      print 9010, TMAX
9010  format(' Equilibrium par coupon rates on bonds of maturity ',F6.3)

      write(*,9000) ( (RMN + I*HR), I = 0, IN, IN/10 )
9000  format('   R:  ', 12(1x,f5.3) )
      print*,' L:'

      do 400 I = 0, IM, IM/10
 400     write(*,9020) (LMN + I*HL), (U(I,J), J=0, IM, IM/10)
9020  format(1x,f5.3,2x,12(1x,f5.2) )
      print*
      go to 100
      end

C***********************************************************************
C**** UTILITY SUBROUTINES CALLED ***************************************
 
      subroutine COPYC( X, U, IM, IN )
      implicit double precision (A-H, K-L, O-Z)
      dimension U(IM,*)
      do 100 I = 1, IM
         do 100 J = 1, IN
 100     U(I,J) = X
      return
      end

      subroutine ADDC( X, U, IM, IN )
      implicit double precision (A-H, K-L, O-Z)
      dimension U(IM,*)
      do 100 I = 1, IM
         do 100 J = 1, IN
 100     U(I,J) = U(I,J) + X
      return
      end

C **********************************************************************
C     DEFINE FUNCTIONS WHICH WILL BE USED TO GENERATE THE
C     COEFFICIENTS IN PDE .
C
C     A*ULL + B*URR + C*ULR + D*UL + E*UR + F*U + COUP(T) + UT = 0
C **********************************************************************
 
      double precision function COEFFX()
      implicit double precision (A-H,J-Z)
      implicit integer (I)
      dimension PARM(*)

      entry FNAKA(IFN,L,R,PARM)
          SIGL = PARM(2)
          FNAKA= SIGL*SIGL*L*L/2.0D0
          return

      entry FNAKB(IFN,L,R,PARM)
          SIGR = PARM(1)
          FNAKB= SIGR*SIGR*R*R/2.0D0
          return

      entry FNAKC(IFN,L,R,PARM)
          SIGR = PARM(1)
          SIGL = PARM(2)
          RHO  = PARM(3)
          FNAKC= RHO*SIGR*SIGL*R*L
          return

      entry FNAKD(IFN,L,R,PARM)
          KAPL = PARM(5)
          LAML = PARM(7)
          LBAR = PARM(8)
          if ( L .gt. 0.0 ) then
             FNAKD = -KAPL * L * log(L/LBAR) - LAML * sqrt(R) * L
          else
             FNAKD = 0.0
          endif
          return

      entry FNAKE(IFN,L,R,PARM)
          SIGR = PARM(1)
          KAPR = PARM(4)
          KAPL = PARM(5)
          LAMR = PARM(6)
          LBAR = PARM(8)
          if ( R .LE. 0.0 ) then
             FNAKE = 0.0
             return
          endif
          if ( L .GT. 0.0 ) then
             G = log(L)
             else
             G = -1.0D60
          endif
          FNAKE= KAPR * R * (G - log(R)) - LAMR * sqrt(R) * R
          return

      entry FNAKF(IFN,L,R,PARM)
          FNAKF= -R
          return
      end

C**** FUNCTIONS CALLED BY ADSTEP GIVING BOUNDARY VALUES IF KNOWN *******
 
      double precision function CXXXF()
      implicit double precision (A-H,K-L,O-Z)
      entry FNULMN(R)
         FNULMN = 0.0
      return
      entry FNULMX(R)
         FNULMX = 0.0
      return
      entry FNURMN(L)
         FNURMN = 0.0
      return
      entry FNURMX(L)
         FNURMX = 0.0
      return
      entry FNLNRN
         FNLNRN = 0.0
      return
      entry FNLXRN
         FNLXRN = 0.0
      return
      entry FNLNRX
         FNLNRX = 0.0
      return
      entry FNLXRX
         FNLXRX = 0.0
      return
      end