function [c] = pure_jump_price(S0, r, T, K, sigma, lambda, alpha)
% Computes the theoretical price of the european call option
% via the pure jump model
%
%[C]
%   Returns the price of the option
%
%PURE_JUMP_PRICE(S0,R,T,K,SIGMA,LAMBDA,ALPHA)
%   S0  Prize of the stock now
%   R   Interest rate
%   T   Maturity date
%   K   Strike price
%   SIGMA   Volatility
%   LAMBDA  Frequency of jumps per unit of time
%   ALPHA   Mean rate of return

%% Input parameters handling
if nargin < 7; alpha = 0.2; 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 = max([S0 - K, 0]);
    return
end

%% Computation task
rn_lambda = lambda - (alpha - r)/sigma;
n = 100;  % summation index (max is 170 due to the n! limit)
c = option_price(S0, rn_lambda, T, r, sigma, K, n);

end

function [y] = positive(x)
    % represents (x)^+ operator
    len_x = length(x);
    x = reshape(x, [len_x, 1]);
    y = max([x zeros(length(x),1)], [], 2);
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, sigma, K, n)
    i = (0:n);
    a = x*exp(-lambda*sigma*T)*(sigma + 1).^i';
    b = K*exp(-r*T);
    p = poisson_prob(lambda, T, n);
    elements = positive(a - b).*p;
    c = sum(elements);
end

