%% EXERCISE 1 — Convergence (Black–Scholes European call)
%   For a European call in Black–Scholes (S0=100, K=100, r=5%, sigma=20%, T=1):
%   (i) Show convergence of the Monte Carlo estimator with respect to the number of paths N.
%   (ii) Show convergence with respect to the number of time steps M.
%   (iii) Plot the function K -> C0(K) for K in [20, 150].
%

clear; clc; close all; rng(2727);

% Parameters
S0    = 100;
r     = 0.05;
sigma = 0.20;
T     = 1.00;

%% Convergence vs number of paths N (fixed M)
M_fix   = 100;
N_list  = [10:100:50000];
C_N     = zeros(size(N_list));


K = 100;  
for i = 1:numel(N_list)
    N  = N_list(i);
    ST = bs_terminal(S0, r, sigma, T, M_fix, N); %calculate ST
    Y  = exp(-r*T)*max(ST - K, 0); %compute the payoff
    C_N(i) = mean(Y); %calculate the price
end
C_BS = bs_call(S0,K,r,sigma,T);

figure('Name','(i) Convergence vs N');
set(gca,'XScale','log'); grid on;
plot(N_list, C_N, 'o-','LineWidth',1.2)
xlabel('N (log scale)'); ylabel('Price estimate');
title('European Call: Convergence vs Number of Paths (Fixed M)');
yline(C_BS,'--','BS closed','LineWidth',1.2,'Color','red');

%% Convergence vs number of time steps M (fixed N)
N_fix   = 5e4;
M_list  = [1:10:300];
C_M     = zeros(size(M_list));

for j = 1:numel(M_list)
    ST   = bs_terminal(S0, r, sigma, T, M_list(j), N_fix);
    Y    = exp(-r*T)*max(ST - K, 0);
    C_M(j) = mean(Y);
end

figure('Name','(ii) Convergence vs M');
plot(M_list, C_M, 'o-','LineWidth',1.2); hold on; grid on;
yline(C_BS,'--','BS closed','LineWidth',1.2);
xlabel('M (time steps)'); ylabel('Price estimate');
title('European Call: Convergence vs Time Steps (N=5e4 fixed)');

%% Function K -> C0(K)
N_price = 1e5; 
M_price = 100;
K_vec   = [20:5:150];
C_mc    = zeros(size(K_vec)); 
C_bs    = C_mc;

for k = 1:numel(K_vec)
    ST      = bs_terminal(S0, r, sigma, T, M_price, N_price);
    C_mc(k) = exp(-r*T)*mean(max(ST - K_vec(k), 0));
    C_bs(k) = bs_call(S0, K_vec(k), r, sigma, T);
end

figure('Name','(iii) K -> C0(K)');
plot(K_vec, C_mc, 'o-','LineWidth',1.2); hold on;
plot(K_vec, C_bs, '--','LineWidth',1.2);
grid on; xlabel('K'); ylabel('C_0(K)');
legend('MC (M=100, N=1e5)','Black–Scholes closed','Location','northeast');
title('Strike Sweep: K --> C_0(K)');

%% Functions 
function ST = bs_terminal(S0, r, sigma, T, M, N)
% Terminal price via exact lognormal steps on a grid of M steps
dt = T/M; drift = (r - 0.5*sigma^2)*dt; vol = sigma*sqrt(dt);
logS = log(S0) + sum(drift + vol.*randn(N,M), 2); 
%sum across columns (over time steps (2)) to get total log-increment for each path.
ST   = exp(logS);
end

function C = bs_call(S0, K, r, sigma, T)
% Black–Scholes closed form
d1 = (log(S0/K) + (r + 0.5*sigma^2)*T) / (sigma*sqrt(T));
d2 = d1 - sigma*sqrt(T);
C  = S0*normcdf(d1) - K*exp(-r*T)*normcdf(d2);
end