clear all
close all
clc

T = 1;
dt = 1/250; % Discretization step size
N = round(T/dt); % Number of discretization steps
M = 1000000; % Number of simulated path
dW = randn(N,M)*sqrt(dt);
W = [zeros(1,M); cumsum(dW(1:end-1,:),1)];
integrand = W.*dW;
PL = sum(integrand);
E_PL = mean(PL)

h = figure()
histogram(PL, 100, 'normalization','pdf')
title (['pdf of RV \int_{0}^{T}(W_tdW_t) for T = ' num2str(T,'%.2f') '    Numerical Expected Value = ' num2str(E_PL,'%.3f') ])