function [c] = jump_diffusion_price(S0, r, T, K, sigma, lambda, beta, b)
% Computes the theoretical price of the european call option
% via the Jump-diffusion model
%
%[C]
%   Returns the price of the option
%
%JUMP_DIFFUSION_PRICE(S0,R,T,K,SIGMA,LAMBDA,BETA,B)
%   S0  Prize of the stock now
%   R   Interest rate
%   T   Maturity date
%   K   Strike price
%   SIGMA   Volatility
%   LAMBDA  Freqeuncy of jumps per unit of time
%   BETA    Mean size of jumps
%   B   Variance of jump's size

%% Input parameters handling
if nargin < 8; b = 0.4; end
if nargin < 7; beta = 0.083287; end
if nargin < 6; lambda = 1; end
if nargin < 5; sigma = 0.3; end
if nargin < 4; K = 100; end
if nargin < 3; T = 1; end
if nargin < 2; r = 0.05; end
if nargin < 1; S0 = 100; end

if T == 0
    c = black_scholes(S0,r,T,K,sigma,'call');
    return
end

%% Computation task
rn_lambda = lambda*(1+beta);  % risk-neutral lambda
n = 50;  % summation index (max is 170 due to the n! limit)
r_n = r_element(r, lambda, beta, T, n);
sigma_n = sigma_element(sigma, b, T, n);
c = option_price(S0, rn_lambda, T, r_n, sigma_n, K, n);

end

function s_n = sigma_element(sigma, b, T, n)
    i = (0:n)';
    s_n = sqrt(sigma^2 + b^2*i./T);
end

function r_n = r_element(r, lambda, beta, T, n)
    i = (0:n)';
    r_n = r - lambda.*(beta) + (log(1+beta)./T).*i;
end

function [p] = poisson_prob(lambda, T, n)
    i = (0:n)';
    p = (lambda.^i.*T.^i)./factorial(i).*exp(-lambda*T);
end

function [c] = option_price(x, lambda, T, r_n, sigma_n, K, n)
    p = poisson_prob(lambda, T, n);
    bsm = black_scholes(x, r_n, T, K, sigma_n, 'call');
    elements = p.*bsm;
    c = sum(elements);
end

