c
c This program computes the price of a European option (call and put) 
c using the simple jump diffusion model of Kou (2002) discussed in Chapter 6 
c of the book Tsay (2002).
c
c The program is intended for use in my course ``Analysis of Financial Time 
c Series'' at GSB. I tested the program, but errors may still exist. Please 
c bring any error to my attention. 
c
c You may use the program at your own risk!  Ruey S. Tsay, December 2000.
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      implicit double precision(a-h,o-z)
c      
      write(6,1)
 1    format(1x,'Current and strike prices: ',$)
      read(5,*) pr, sk
      write(6,2)
 2    format(1x,'time to expiration (in years): ',$)
      read(5,*) dt
      write(6,3)
 3    format(1x,'Interest rate and volatility (std): ',$)
      read(5,*) free, sig
      if(sig .le. 0.0d0)then
       print*,'volatility must be positive!'
       stop
      endif
      write(6,4)
 4    format(1x,'lambda, kappa, eta: ',$)
      read(5,*) poi, rka, eta
      if((eta.le.0.0d0).or.(eta.gt.1.0d0))then
       print*,'eta is out of range!'
       stop
      endif
c
      nterm = 10
      write(6,5)
 5    format(1x,'terms of truncation? (default = 10): ',$)
      read(5,*) nterm
      if(nterm .le. 10) nterm = 10
      if(nterm .gt. 16) nterm = 16
c
      psi = dexp(rka)/(1.0d0-eta*eta)-1.0d0
      prsk = dlog(pr/sk)
      sigsrt = sig*dsqrt(dt)
      sig2h = sig*sig/2.0d0
      tlamb = poi*dt
      rtexp = dexp(-free*dt)
      ome1 = -prsk+tlamb*psi-(free-sig2h)*dt
      w1 = poi*psi*dt
      sigt = dexp(sig*sig*dt/(2.0d0*eta*eta))
      st2pi = dsqrt(6.28318530718d0)
      steta = sigsrt/eta
c
c----- compute the portion of geometric Brownian motion
      hp = (prsk+(free+sig2h-poi*psi)*dt)/sigsrt
      hm = (prsk+(free-sig2h-poi*psi)*dt)/sigsrt
      call gauss(hp,cdf1)
      call gauss(hm,cdf2)
      cgeo = dexp(-tlamb)*(pr*dexp(-w1)*cdf1-sk*rtexp*cdf2)
      print*,'Price due to the geometric part: ', cgeo
c---- compute the price for the jump process portion
c
      cj = 0.0d0
c
      do 100 n = 1, nterm
       tmp = dfloat(n)*rka/sigsrt
       bp = hp+tmp
       bm = hm+tmp
       ome = ome1-dfloat(n)*rka
       tmp1 = ome/sigsrt
       cp = steta+tmp1
       cm = steta-tmp1
       omexp = dexp(ome/eta)
       omexpi = dexp(-ome/eta)
       a11 = dexp(-w1+dfloat(n)*rka)*pr/2.0d0
       a22 = rtexp*omexpi*sigt*sk/2.0d0
       a33 = rtexp*omexp*sigt*sk/2.0d0
c
       call gauss(bp,cdf1)
       call gauss(bm,cdf2)
c
       fac = factor(n)
       npow = 2**(2*n-1)
       coef = (tlamb**n)*dexp(-tlamb)/fac
c
       a = 0.0d0
       do 50 j=1, n
        m = 2*n-j-1
        call comb(m,n-1,icnt)
        icnt = icnt*(2**j)
        coef1 = coef*dfloat(icnt)/dfloat(npow)
        a1 = (1.0d0/(1.0d0-eta)**j+1.0d0/(1.0d0+eta)**j)*a11*cdf1 
     &   -rtexp*sk*cdf2
        a2 = 0.0d0
        a3 = 0.0d0
        do 30 ii=1, j
         i = ii-1
          w2 = 1.0d0/(1.0d0-eta)**(j-i)-1.0d0
          w2 = w2*(steta**i)*hh(cm,i)/st2pi
          w3 = 1.0d0-1.0d0/(1.0d0+eta)**(j-i)
          w3 = w3*(steta**i)*hh(cp,i)/st2pi
          a2 = a2 + w2
          a3 = a3 + w3
 30      continue
         a2 = a2*a22
         a3 = a3*a33
         a = coef1*(a1+a2+a3)+a
 50       continue
         cj = cj + a
 100      continue
c
      print*,'price due to jump component: ', cj
      c = cj + cgeo
      print*,'Price of a call: ', c
c
c---- put option
c
      p = c+sk*dexp(-free*dt)-pr
c 
      print*,'Price of a put: ', p
c
      stop
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine comb(m,n,icnt)
c----- combites the combinations      
      implicit double precision(a-h,o-z)
      if((m .lt. n).or.(m .lt. 0))then
       icnt = 0
       return
      endif
      if(m.eq.n)then
       icnt = 1
       return
      endif
c
      n1 = 1
      do 5 j=n, 1, -1
 5       n1 = n1*j
c
      icnt = 1
      do 10 i=m, (m-n+1), -1
 10      icnt =icnt*i
      icnt = icnt/n1
      return
      end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision function factor(n)
      implicit double precision(a-h,o-z)
c
      if(n .lt. 0)then
        factor = 0.0d0
        return
       endif
       if(n.le. 1)then
        factor = 1.0d0
        return
       endif
       a = 1.0d0
       do 10 i=n, 1, -1
 10       a = a*dfloat(i)
       factor = a
c
       return
       end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision function hh(x,n)
c
      implicit double precision(a-h,o-z)
      dimension h(100)
      if(n .gt. 100)then
       print*,'dimension in Hn function is too low!'
       stop
      endif
c
      pisqt = dsqrt(6.28318530718d0)
      h(1) = dexp(-x*x/2.0d0)
      if(n .eq. -1)then
       hh = h(1)
       return
      endif
      tmp = -x
      call gauss(tmp,cdf)
      h(2) = pisqt*cdf
      if(n .eq. 0)then
        hh = h(2)
        return
      endif
      do 40 j=3, 2+n
 40        h(j) = (h(j-2)-x*h(j-1))/dfloat(j-2)
 50        continue
c
       hh = h(n+2)
c
       return
       end
ccccccccccccccccccccccccccccccccccccccc
       subroutine gauss(x,cdf)
c
       implicit double precision(a-h,o-z)
       dimension a(5)
c
       a(1) = 0.319381530d0
       a(2) = -0.356563782d0
       a(3) = 1.781477937d0
       a(4) = -1.821255978d0
       a(5) = 1.330274429d0
c
c---------------- n: is used to indicate negative value.
       n = 0
c
       if(x .lt. 0.0d0) then
        x = -x 
        n = 1
       endif
       xk = 1.0d0/(1.0d0+0.2316419d0*x)
       xf = dexp(-x*x/2.0d0)/dsqrt(6.28318530718d0)
       tmp = 0.0d0
       do 10 i=1, 5
 10       tmp = tmp + a(i)*xk**i
       cdf = 1.0d0-xf*tmp
       if(n.eq.1)cdf = 1.0d0-cdf
c
       return
       end
cccccccccccccccccccccccccccccccccccccc