%% Simulation

t0 = 0;  % starting time
T = 1;  % maturity date in years
sigma = 0.3; % volatility
alpha = 0.3;
K = 100; % strike price
lambda = 1;
beta = 0.083287; % 1+beta = mean of log-normal jump
b = 0.4; % variance of log-normal jump
r = 0.05;

n = 1000;% steps
n_instances = 1000;% instances

St = (50:150);
len_St = length(St);

% BSM
c_bs_e = zeros(len_St,1); % empirical prices
c_bs_t = zeros(len_St,1); % theoretical prices
% PJ
c_pj_e = zeros(len_St,1);
c_pj_t = zeros(len_St,1);
% JD
c_jd_e = zeros(len_St,1);
c_jd_t = zeros(len_St,1);

for idx = 1:len_St
    [W, t] = wiener_process(n,T,false,n_instances);
    S_bsm = St(idx)*exp((r-0.5*sigma^2)*repmat(t',1,n_instances) + sigma*W);
    
    c_bs_e(idx) = exp(-r*T)*mean(max(S_bsm(end,:)-K,0));
    c_bs_t(idx) = black_scholes(St(idx),r,T,K,sigma,'call');

    % Pure Jump
    S_pj = zeros(length(t),n_instances);
    rn_lambda = lambda - (alpha - r)/sigma;
    sigma = 0.3;
    jj = zeros(n_instances,1);
    for k = 1:n_instances
        [Tn, Eps, N, T] = poisson_process(0, rn_lambda, T, false); % times of jumps
        pow = 0;
        jj(k) = length(Tn)-1;
        Tn = [Tn(2:end); T];
        for i = 1:length(t)
            for j = 2:length(Tn)
                if Tn(j-1) <= t(i) && t(i) < Tn(j) 
                    pow = j-1;
                end
            end
            S_pj(i,k) = St(idx)*(sigma + 1)^pow*exp(r*t(i) ...
                - rn_lambda*sigma*t(i));
        end
    end
    c_pj_e(idx) = exp(-r*T)*mean(max(S_pj(end,:)-K,0));
    c_pj_t(idx) = pure_jump_price(St(idx),r,T,K,sigma,lambda,alpha);

    % Jump-diffusion
    S_jd = zeros(length(t),n_instances);
    jj = zeros(n_instances,1);
    rn_lambda = lambda*(1+beta);
    % normal to log-normal tranformation with given params of the LN dist.
    ln_mu = log((1+beta)^2/(sqrt((1+beta)^2 + b^2)));
    ln_sigma = sqrt(log(1 + b^2/(1+beta)^2));
    tes = [];
    mm = ones(n_instances, 1);
    for k = 1:n_instances
        % times of jumps
        [Tn, Eps, N, T] = poisson_process(0, rn_lambda, T, false); 
        mult = 1;  % multiplicator
        jj(k) = length(Tn)-1;
        for kk = 1:jj(k)
            pom = exp(ln_mu + ln_sigma*randn());
            tes = cat(1,tes,mult);
            mult = pom*mult;
            mm(k) = mult;
        end
        for i = 1:length(t)
        S_jd(i,k) = St(idx)*exp(sigma*W(i,k) + ...
                (r - beta*lambda - 0.5*sigma^2)*t(i))*mult;
        end
    end
    c_jd_e(idx) = exp(-r*T)*mean(max(S_jd(end,:)-K,0));
    c_jd_t(idx) = jump_diffusion_price(St(idx),r,T,K,sigma,lambda,beta,b);

end

plot(St,c_pj_e); hold on; plot(St,c_pj_t);
plot(St,c_jd_e); hold on; plot(St,c_jd_t);
plot(St,c_bs_e); hold on; plot(St,c_bs_t);

title('Simulated and theoretical option price')
xlabel('\it spot price')
ylabel('option price')
grid on;
lgd = legend("Pure Jump - simulation", ...
             "Pure Jump - theoretical", ...
             "Jump-diff - simulation", ...
             "Jump-diff - theoretical", ...
             "BSM - simulation", ...
             "BSM - theoretical", ...
             "Location", "northwest");
title(lgd, 'Simulation vs. theory')




