% GBM_hist.m
% Use Euler-Maruyama scheme to draw price paths for equity obeying
%  Geometric Brownian Motion SDE: dS=mu*S*dt+sigma*S*dW, and show histogram
%  of distribution of x(T), where x=ln(S).
clear

Nrealiz=5000
S=100 % stock price at time 0
sig=0.2 % volatility
T=1 % end-time (in years)
r=0.06 % risk-free interest rate
div=0.03 % dividend yield
N=500 % number of time steps

dt=T/N;
nudt=(r-div-0.5*sig^2)*dt;
sigsdt=sig*sqrt(dt);

x0=log(S); % initial value for x=ln(S)

driftT=(r-div-0.5*sig^2)*T;
histstart=x0+driftT-2;
histend=x0+driftT+2;
histstep=0.05;
holdendx=zeros(Nrealiz,1); % vector to hold values of x(T) from each path
tic
for i=1:Nrealiz
    x=zeros(N,1);
    t=zeros(N,1);
    x(1)=x0; % initial condition
    t(1)=0;
    r=randn(N,1);
    for n=1:(N-1)
        x(n+1)=x(n)+nudt+sigsdt*r(n);
        t(n+1)=t(n)+dt;
    end
holdendx(i)=x(N);
end

binedges=histstart:histstep:histend;
hist_unnorm=hist(holdendx,binedges);
pdf=hist_unnorm/(length(holdendx)*histstep);
figure
plot(binedges,pdf,'b.');
hold on
gaussian=1/sqrt(2*pi*sig^2*T)*exp(-(binedges-x0-driftT).^2/(2*sig^2*T));
plot(binedges,gaussian,'r-')
toc