%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Probability Density of Standard Brownian Motion
%
%  This code simulates n trials of standard Brownian motion over
%  the interval [0,T]. The code then plots the density of the realizations
%  with N bins. 
%
%  Inputs:
%  1. n=number of realization to simulate
%  2. T=maximum time
%  3. N=number of bins used in histogram
%  4. M=number of temporal points used for Brownian motion
%
%  Ouputs:
%   1. Plot of emperical density and the exact density
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Inputs
n=100000;
T=1;
N=100;
M=100;

%% Construction of Discretization
t=linspace(0,T,M);
dt=t(2)-t(1);

B=zeros(1,M);
Z=zeros(1,n);

%% Initilizing Plot
figure('units','normalized','outerposition',[0 0 1 1]);
hold on;
xlabel('$$B_T$$','interpreter','latex')
ylabel('$$f(b)$$','interpreter','latex')
ax=gca;
ax.FontSize = 40;

%% Simulation all of the Brownian motions and recording there endpoint.
for j=1:n
    for i=1:M-1
        B(i+1)=B(i)+sqrt(dt)*randn;
    end
    Z(j)=B(end);
end

%% Plotting Everything
histogram(Z,M,'Normalization','pdf');

x=linspace(-3*T,3*T,n);
y=1/sqrt(2*pi*T)*exp(-x.^2/(2*T));
plot(x,y,'linewidth',4);


