clear; clc;

s0 = 250;
X = 250;
r = 0.01;
v = 0.25;
T = 1;

M = 1e5;
c = 0;
for i = 1 : M
    ST = s0 * exp((r - 0.5 * v ^ 2) * T + v * sqrt(T) * randn(1));
    payoff = ST - X;
    if payoff > 0
        c = c + payoff;
    end
end
c = c / M * exp(-r * T)

%% closed-form formula in built-in function
[Call, Put] = blsprice(250, 250, 0.01, 1, 0.25);

%% vectorization
ST = s0 * exp((r - 0.5 * v ^ 2) * T + v * sqrt(T) * randn(1, M));
c = sum(ST(ST > X) - X) / M * exp(-r * T)



