
r     = 0.05; 
sig   = 0.5;
T     = 2.5;
S0    = 4;
K     = 3;

lam_j = 1;
mu_j  = 0.1;
del_j = 0.4;

eps = 1e-20;

v = merton_call(lam_j,mu_j,del_j,r,sig,T,S0,K,'value');

disp(sprintf('value = %f ',v))

disp(sprintf('v1 (analytic), v2 (CVT)'))

v1 = merton_call(lam_j,mu_j,del_j,r,sig,T,S0,K,'delta');
v2 = merton_call(lam_j,mu_j,del_j,r,sig,T,S0+i*eps,K,'value');
v2 = imag(v2/eps);
disp(sprintf('delta: v1, v2 = %f %f',v1,v2))

v1 = merton_call(lam_j,mu_j,del_j,r,sig,T,S0,K,'vega');
v2 = merton_call(lam_j,mu_j,del_j,r,sig+i*eps,T,S0,K,'value');
v2 = imag(v2/eps);
disp(sprintf('vega:  v1, v2 = %f %f',v1,v2))

v1 = merton_call(lam_j,mu_j,del_j,r,sig,T,S0,K,'gamma');
v2 = merton_call(lam_j,mu_j,del_j,r,sig,T,S0+i*eps,K,'delta');
v2 = imag(v2/eps);
disp(sprintf('gamma: v1, v2 = %f %f',v1,v2))

v1 = merton_call(lam_j,mu_j,del_j,r,sig,T,S0,K,'mu_j');
v2 = merton_call(lam_j,mu_j+i*eps,del_j,r,sig,T,S0,K,'value');
v2 = imag(v2/eps);
disp(sprintf('mu_j:  v1, v2 = %f %f',v1,v2))

v1 = merton_call(lam_j,mu_j,del_j,r,sig,T,S0,K,'del_j');
v2 = merton_call(lam_j,mu_j,del_j+i*eps,r,sig,T,S0,K,'value');
v2 = imag(v2/eps);
disp(sprintf('del_j: v1, v2 = %f %f',v1,v2))

v1 = merton_call(lam_j,mu_j,del_j,r,sig,T,S0,K,'lam_j');
v2 = merton_call(lam_j+i*eps,mu_j,del_j,r,sig,T,S0,K,'value');
v2 = imag(v2/eps);
disp(sprintf('lam_j: v1, v2 = %f %f',v1,v2))
