%% Simulation

t0 = 0;  % starting time
T = linspace(0,1,50);  % maturity date in years
sigma = 0.3; % volatility
alpha = 0.3;
K = 100; % strike price
lambda = 1;
beta = 0.083287;
b = 0.4;
r = 0.05;
S0 = 100;

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

len_T = length(T);

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

for idx = 1:len_T
    [W, t] = wiener_process(n,T(idx),false,n_instances);
    S_bsm = S0*exp((r-0.5*sigma^2)*repmat(t',1,n_instances) + sigma*W);
    
    c_bs_e(idx) = exp(-r*T(idx))*mean(max(S_bsm(end,:)-K,0));
    c_bs_t(idx) = black_scholes(S0,r,T(idx),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, ~] = poisson_process(0, rn_lambda, T(idx), false); % times of jumps
        pow = 0;
        jj(k) = length(Tn)-1;
        Tn = [Tn(2:end); T(idx)];
        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) = S0*(sigma + 1)^pow*exp(r*t(i) - rn_lambda*sigma*t(i));
        end
    end
    c_pj_e(idx) = exp(-r*T(idx))*mean(max(S_pj(end,:)-K,0));
    c_pj_t(idx) = pure_jump_price(S0,r,T(idx),K,sigma,lambda,alpha);

    % Jump-diffusion
    S_jd = zeros(length(t),n_instances);
    jj = zeros(n_instances,1);
    alpha = 0.3;
    beta = 0.083287; % mu of lognormal = 1 + beta
    sigma = 0.3;
    lambda = 1;
    b = 0.4; % st. dev. of lognormal
    rn_lambda = lambda*(1+beta);
    % normal to lognormal tranformation
    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
        [Tn, Eps, N, ~] = poisson_process(0, rn_lambda, T(idx), false); % times of jumps
        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) = S0*exp(sigma*W(i,k) + ...
                (r - beta*lambda - 0.5*sigma^2)*t(i))*mult;
        end
    end
    c_jd_e(idx) = exp(-r*T(idx))*mean(max(S_jd(end,:)-K,0));
    c_jd_t(idx) = jump_diffusion_price(S0,r,T(idx),K,sigma,lambda,beta,b);

end

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

title('Simulated and theoretical option price')
xlabel('\it time to maturity')
ylabel('option price')
grid on;
lgd = legend("Pure Jump - simulation", ...
             "Pure Jump - theoretical", ...
             "Jump-diff - simulation", ...
             "Jump-diff - theoretical", ...
             "BSM - simulation", ...
             "BSM - theoretical", ...
             "Location", "northeast");
title(lgd, 'Simulation vs. theory')
set(gca, 'XDir','reverse')