%% Comparison between the binomial formula weights and lognormal

%%%%%%%%%%%%% Problem parameters %%%%%%%%%%%%%%
S0 = 100; tau = 5; r = 0.05; sigma = 0.3; 
N = 100;  dt = tau / N;
mu = 0.1;
u = exp(mu*dt + sigma*sqrt(dt));
d = exp(mu*dt - sigma*sqrt(dt));
q = (exp(r*dt)-d) / (u-d); % "risk-neutral probability"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Final level of the tree
ST = S0 * d.^(N - (0:N)) .* u.^(0:N);

% Calculate weights (the function nchoosek doesn't take vectors)
weights = zeros(1,N+1);
weights(1) = (1-q)^N*1;
for j=1:N
    weights(j+1) = weights(j) * (q/(1-q)) * (N-j+1)/j;
end
spc = ST * (sqrt(u/d) - sqrt(d/u));

% Plot the weights we want to approximate
figure;
plot(ST, weights./spc, '*', 'Linewidth', 3);
xlim([0, 500]);  hold on;

%% Plot the histogram of a simulated (!) lognormal variable
z = randn(1,1e6);
ss = S0*exp((mu-sigma^2/2)*tau + sigma*sqrt(tau)*z);
edges = [ST*sqrt(d/u), ST(end)*sqrt(u/d)];
histogram(ss, edges, 'Normalization', 'pdf');

%% Try to fit lognormal pdf
lognormal = @(y, m, s) 1./ (s*y*sqrt(2*pi)) ...
    .* exp(- (log(y) - m).^2 / (2*s^2));
plot(ST, lognormal(ST/S0, (mu-sigma^2/2)*tau,  sigma*sqrt(tau))/S0, 'k');


%% Look for a better fit

plot(ST, lognormal(ST/S0, (0.05-sigma^2/2)*tau,  sigma*sqrt(tau))/S0, 'k');
