%
% function V = merton_call(lam_j,mu_j,del_j,r,sigma,T,S,K,opt)
%
% European call based on Merton model of geometric Brownian motion
% plus jumps with a Poisson rate and log-Normal size
%
% lam_j - jump frequency
% mu_j  - log of jump mean
% del_j - variability of jump
% r     - interest rate
% sigma - volatility
% T     - time interval
% S     - asset value(s)
% K     - strike price(s)
% opt   - 'value', 'delta', 'gamma', 'vega', 'lam_j', 'mu_j', 'del_j'
% V     - option value(s)
%

function V = merton_european_call(lam_j,mu_j,del_j,r,sigma,T,S,K,opt)

if nargin ~= 9
  error('wrong number of arguments');
end

V = 0;

for n = 0:50
  sig_n = sqrt(sigma^2 + n*del_j^2);
  S0_n  = S*exp( (1-exp(mu_j))*lam_j*T + n*mu_j);

  if n==0
    prob = exp(-lam_j*T);
  elseif lam_j==0
    prob = 0;
  else
    prob = exp(-lam_j*T + n*log(lam_j*T) - sum(log(1:n)));
  end

  switch opt
    case 'value'
      val = european_call(r,sig_n,T,S0_n,K,opt);
      V = V + prob*val;
    case 'delta'
      val = european_call(r,sig_n,T,S0_n,K,opt);
      V = V + prob*val*S0_n./S;
    case 'gamma'
      val = european_call(r,sig_n,T,S0_n,K,opt);
      V = V + prob*val*(S0_n./S).^2;
    case 'vega'
      val = european_call(r,sig_n,T,S0_n,K,opt);
      V = V + prob*val*sigma/sig_n;

    case 'lam_j'
      if n==1
        prob2 = exp(-lam_j*T);
      elseif lam_j==0 | n==0
        prob2 = 0;
      else
        prob2 = exp(-lam_j*T + (n-1)*log(lam_j*T) - sum(log(1:n-1)));
      end
      val = european_call(r,sig_n,T,S0_n,K,'value');
      V = V + T*(prob2-prob)*val;
      val = european_call(r,sig_n,T,S0_n,K,'delta');
      V = V + prob*val*S0_n*(1-exp(mu_j))*T;
    case 'mu_j'
      val = european_call(r,sig_n,T,S0_n,K,'delta');
      V = V + prob*val*S0_n*(-exp(mu_j)*lam_j*T + n);
    case 'del_j'
      val = european_call(r,sig_n,T,S0_n,K,'vega');
      V = V + prob*val*n*del_j/sig_n;
  end
end
