Hermite Polynomials and Gaussian distribution Links

Gaussian functions from Wiki

Gaussian integrals from Wiki

Normal distribution from Wiki

Construct Hermite Polynomials from Gaussian distribution functions (see copy below)

Hermite polynomials from Wiki

Addition on Hermite polynomials


Construct Hermite Polynomials from Gaussian distribution functions, copied from link above:

Let us look at the Hermite polynomials in somewhat more detail. Consider a Gaussian function, .

.

We write

.

This defines the Hermite polynomials.

.

The parity of the Hermite polynomials is (-1)ⁿ.

H_n(z) is a nth degree polynomial in z.

Proof:

This statement holds for n = 1 and n = 2.

Let . Then

.

If H_n-1(z) is a polynomial of degree n-1, the H_n(z) is a polynomial of degree n. The statement therefore holds for all n.

The generating function

Consider F(z+l) = exp(-(z+l)²).

A Taylor series expansion, , yields

.

Let l’ = -l. Then

.

. (Taylor series expansion)

is called the generating function for the Hermite polynomials H_n(z).

Recurrence relations

We have already found

.

Other recurrence relations exist.

and combining the first two relations

.

[ is obtained by differentiating with respect to z,

, and equating terms of equal power of l.

is obtained by differentiating with respect to l.]

Summary

The Hermite polynomials H_n(z) are defined by

or .

The generating function of the Hermite polynomials is

.

The recurrence relations are

, ,

and

.

The differential equation satisfied by the Hermite polynomials is

.

The normalization and orthogonality condition are

Test strong convergence of Milstein Mthd for SDE

Below is the file

******************************************

%MILSTRONG  Test strong convergence of Milstein: vectorized
%
% Solves   dX = r*X*(K-X) dt + beta*X dW,  X(0) = Xzero,
%        where r = 2, K= 1, beta = 1 and Xzero = 0.5.
%
% Discretized Brownian path over [0,1] has dt = 2^(-11).
% Milstein uses timesteps 128*dt, 64*dt, 32*dt, 16*dt (also dt for reference).
%
% Examines strong convergence at T=1:  E | X_L - X(T) |.
% Code is vectorized: all paths computed simultaneously.

randn('state',100)
r = 2; K = 1; beta = 0.25; Xzero = 0.5;   % problem parameters
T = 1; N = 2^(11); dt = T/N;                       % from 0 to T=1. N = 2^11 steps, each of length dt            
M = 500;                                                            % number of paths sampled
R = [1; 16; 32; 64; 128];                               % Milstein stepsizes are R*dt, 5 R's for 5 paths

dW = sqrt(dt)*randn(M,N);                 % Brownian increments for each path each timestep
Xmil = zeros(M,5);                                  % preallocate array: M paths w. corresp. 5 R
for p = 1:5                                                   % Consider each time step separately
     Dt = R(p)*dt; L = N/R(p);                % Dt = new time step
                                                                       % L timesteps of size Dt = R dt
     Xtemp = Xzero*ones(M,1);            % All M paths begin at Xzero
     for j = 1:L                          % For each step
          Winc = sum(dW(:,R(p)*(j-1)+1:R(p)*j),2); 
          Xtemp = Xtemp + Dt*r*Xtemp.*(K-Xtemp) + beta*Xtemp.*Winc ...
              + 0.5*beta^2*Xtemp.*(Winc.^2 - Dt); % See below
end
     Xmil(:,p) = Xtemp;       % store Milstein solution at t = 1
end

Xref = Xmil(:,1);                                       % Reference solution: "exact" soln when timestep = dt 
Xerr = abs(Xmil(:,2:5) - repmat(Xref,1,4));   % Error in each path
mean(Xerr);                                              % Mean pathwise errors
Dtvals = dt*R(2:5);                                % Milstein timesteps used

subplot(224)                                            % lower RH picture
loglog(Dtvals,mean(Xerr),'b*-'), hold on
loglog(Dtvals,Dtvals,'r--'), hold off  % reference slope of 1
axis([1e-3 1e-1 1e-4 1])
xlabel('\Delta t')
ylabel('Sample average of | X(T) - X_L |')
title('milstrong.m','FontSize',10)

%%%% Least squares fit of error = C * Dt^q %%%%
A = [ones(4,1), log(Dtvals)]; rhs = log(mean(Xerr)');
sol = A\rhs; q = sol(2)
resid = norm(A*sol - rhs)

******************************************

Notes:

-

-

Brownian Path Simulation using Matlab

%BPATH1  Brownian path simulation

randn('state',100)           % set the state of randn to be 100
T = 1; N = 500; dt = T/N;    % Length = 1, do 500 steps, each step measures dt=1/500
dW = zeros(1,N);             % preallocate arrays ...
W = zeros(1,N);              % for efficiency

dW(1) = sqrt(dt)*randn;      % first approximation outside the loop ...
W(1) = dW(1);                % since W(0) = 0 is not allowed
for j = 2:N
 dW(j) = sqrt(dt)*randn;   % general increment, dW ~ sqrt(dt)
 W(j) = W(j-1) + dW(j);    % then go a step of dW(j) from W(j-1) to W(j)
end

plot([0:dt:T],[0,W],'r-')    % plot W against t: plot from time 0 to time T, increment dt
xlabel('t','FontSize',16)
ylabel('W(t)','FontSize',16,'Rotation',0)

Note: Matlab starts array from 1, not 0. So in order to include the starting point W(0) = 0 into the graph, we need to include it in the plotting command by doing plot([0:dt:T],[0,W])

Following is a more elegant and more efficient version of bpath1.m :
The "for" loop is replaced by a "vectorized" command
(KEY to writing efficient MATLAB)code (See DJHigham & NJHigha, MATLAB guide, SIAM, Philly 2000)


%BPATH2  Brownian path simulation: vectorized

randn('state',100)          % set the state of randn
T = 1; N = 500; dt = T/N;

dW = sqrt(dt)*randn(1,N);   % generate the array of increments
W = cumsum(dW);             % cumulative sum, W(j)= dW(1) + ... + dW(j)

plot([0:dt:T],[0,W],'r-')   % plot W against t, red line
xlabel('t','FontSize',16)
ylabel('W(t)','FontSize',16,'Rotation',0)

A higher simulation of Brownian path:
Average from about 1000 different paths. Sketch the average path and 5 of the brownian paths.
Here we evaluate the function u(W(t)) = exp(t + W(t)/2) along 1000 discretized Brownian paths.
The average of u(W(t)) over these paths is plotted with a solid blue line.
Five individual paths are also plotted using a dashed red line.

%BPATH3  Function along a Brownian path

randn('state',100)                                % set the state of randn
T = 1; N = 500; dt = T/N; t = [dt:dt:1];          % t = row vector, N components, from dt to 1, increments dt.

M = 1000;                                         % M paths simultaneously
dW = sqrt(dt)*randn(M,N);                         % increments for M paths. 1 path 1 row. Each row N increments ...
                                                % i.e., dW is an M-by-N array st dW(i,j) gives the increment dWj

W = cumsum(dW,2);                                 % cumulative sum accross the second dimension (column)
                                                % i.e., each row of W is the cumsum of the corresp. row of dW

U = exp(repmat(t,[M 1]) + 0.5*W);                 % repmat(t,[M 1]) form a "matrix" M row 1 column, each "entry" is t.
                                                % With t a N-dim row vector, this put M t's stack together,
                                                % forming a M by N matrix.
                                                % then add it to the matrix .5*W
                                                  % then take exp of this matrix
Umean = mean(U);                                  % compute columwise averages
plot([0,t],[1,Umean],'b-'), hold on               % plot mean over M paths
                                                 % x-axis is 0 and the vector t
                                                 % y-axis is 1. At time 0 the path begin at y=1. WHY???
plot([0,t],[ones(5,1),U(1:5,:)],'r--'), hold off  % plot 5 individual paths
                                                 % discrete "time" from 0 to 1
                                                 % over 5 paths, all begin at 1
                                                 % paths are first 5 rows of U
xlabel('t','FontSize',16)
ylabel('U(t)','FontSize',16,'Rotation',0,'HorizontalAlignment','right')
legend('mean of 1000 paths','5 individual paths',2)

averr = norm((Umean - exp(9*t/8)),'inf')          % sample error. Infinity norm of matrix
                                                 % largest row sum of A.
                                                 % w.o. 'inf' - 1-norm, largest column sum

To be explored....