Analytic method exploits change of derivatives formula for densities to calculate derivatives of lognormal random variable with respect to standard normal variable.

Theoretically you could expand normal density as a Series in lognormal density meaning you can expand one random variable whose density has higher analytic derivatives as an expansion with respect to any other random variable whose higher derivatives also exist. I am sure there could be some restrictions but I am sure you can expand many random variables with different densities as a series in random variables with other densities. It is not necessary that expansion would be in the form of standard normal as we have been doing. though normal random variable is ideal for expanding other random variables but many other random variable types could just equally be used for expansion.

I will give details of the method in next post.

Here is the simple program with brief comments that expands a lognormal variable as a Z-series in standard normals and then compares the density found from Z-series to density from analytic closed form formula of the lognormal density.

.

.

.

Code: Select all

```
function [] = LognormalZSeries()
%This Z-series is for lognormal density and not for lognormal diffusions.
%For lognormal diffusions, you would have to subtract half of squared
%volatility from mean to take effective mean.
%Please specify vol and mean term for lognormal density.
sigma=1;
mu=.5;
%Our expansion is around median of densities so we find the median of
%lognormal density.
Xmedian=exp(mu);
%Below we calculate the value of density and its derivatives at median.
px=1/(sqrt(2* pi)*sigma* Xmedian)* exp(-.5* (log(Xmedian)-mu)^2/sigma^2);
pexp=exp(-.5* (log(Xmedian)-mu)^2/sigma^2);
dpxdX=pexp *(0.398942* mu-0.398942* sigma^2-0.398942 *log(Xmedian))/(sigma^3* Xmedian^2);
d2pxdX2=(pexp* (0.398942* mu^2-0.398942* sigma^2-1.19683* mu *sigma^2+0.797885* sigma^4+(-0.797885* mu+1.19683 *sigma^2) *log(Xmedian)+0.398942* log(Xmedian).^2))/(sigma^5 *Xmedian^3);
d3pxdX3=(1/(sigma^7* Xmedian^4))*pexp* (0.398942 *mu^3-1.19683* mu* sigma^2-2.39365* mu^2* sigma^2+2.39365* sigma^4+4.38837* mu* sigma^4-2.39365* sigma^6+(-1.19683 *mu^2+1.19683 *sigma^2+4.78731* mu* sigma^2-4.38837* sigma^4)* log(Xmedian)+(1.19683 *mu-2.39365* sigma^2)* (log(Xmedian)).^2-0.398942* (log(Xmedian)).^3)
d4pxdX4=(1/(sigma^9* Xmedian^5))*pexp* (0.398942* mu^4 - 3.98942* mu^3* sigma^2 + 1.19683* sigma^4 - 13.963* sigma^6 + 9.57461 *sigma^8 + mu^2* (-2.39365* sigma^2 + 13.963* sigma^4) + ...
mu *(11.9683 *sigma^4 - 19.9471 *sigma^6) + (-1.59577* mu^3 + 4.78731 *mu *sigma^2 + 11.9683* mu^2* sigma^2 - 11.9683* sigma^4 - 27.926* mu* sigma^4 + 19.9471* sigma^6)* log(Xmedian) + (2.39365* mu^2 - 2.39365* sigma^2 - 11.9683* mu* sigma^2 + 13.963* sigma^4)* (log(Xmedian)).^2 + ...
(-1.59577* mu + 3.98942* sigma^2)* (log(Xmedian)).^3 + 0.398942* (log(Xmedian))^4);
d5pxdX5=(1/(sigma^11* Xmedian^6))*pexp* (0.398942* mu^5-5.98413 *mu^4* sigma^2+sigma^6* (-17.9524+89.762* sigma^2-47.8731* sigma^4)+ ...
mu^3* (-3.98942* sigma^2+33.9101* sigma^4)+mu^2 *(35.9048* sigma^4-89.762* sigma^6)+ mu* (5.98413* sigma^4-101.73 *sigma^6+109.31* sigma^8)+ ...
(-1.99471* mu^4+23.9365* mu^3* sigma^2-5.98413* sigma^4+101.73* sigma^6-109.31* sigma^8+mu^2* (11.9683* sigma^2-101.73* sigma^4)+ ...
mu* (-71.8096* sigma^4+179.524* sigma^6))* log(Xmedian)+ ...
(3.98942* mu^3-11.9683* mu* sigma^2-35.9048* mu^2* sigma^2+35.9048* sigma^4+101.73* mu* sigma^4-89.762* sigma^6)* (log(Xmedian)).^2+ ...
(-3.98942* mu^2+3.98942* sigma^2+23.9365* mu* sigma^2-33.9101* sigma^4)* (log(Xmedian)).^3+ ...
(1.99471* mu-5.98413* sigma^2)* (log(Xmedian)).^4-0.398942* (log(Xmedian)).^5);
d6pxdX6=1/(sigma^13* Xmedian^7)*pexp* (0.398942* mu^6 - 8.37779 *mu^5* sigma^2 + mu* sigma^6* (-125.667 + 879.668 *sigma^2 - 703.734* sigma^4) + ...
mu^4* (-5.98413* sigma^2 + 69.8149* sigma^4) + mu^3* (83.7779 *sigma^4 - 293.223* sigma^6) + sigma^6 *(-5.98413 + 209.445* sigma^2 - 647.882* sigma^4 + ...
287.238 *sigma^6) + mu^2 *(17.9524 *sigma^4 - 418.889* sigma^6 + 647.882 *sigma^8) + (-2.39365* mu^5 + 41.8889 *mu^4* sigma^2 + ...
mu^3 *(23.9365 *sigma^2 - 279.26 *sigma^4) + sigma^6* (125.667 - 879.668* sigma^2 + 703.734* sigma^4) + mu^2 *(-251.334 *sigma^4 + 879.668 *sigma^6) + ...
mu *(-35.9048 *sigma^4 + 837.779 *sigma^6 - 1295.76 *sigma^8))* log(Xmedian) + ...
(5.98413* mu^4 - 83.7779 *mu^3 *sigma^2 + 17.9524* sigma^4 - 418.889 *sigma^6 + 647.882* sigma^8 + mu^2* (-35.9048 *sigma^2 + 418.889 *sigma^4) + ...
mu *(251.334 *sigma^4 - 879.668* sigma^6))* (log(Xmedian)).^2 + ...
(-7.97885* mu^3 + 23.9365* mu *sigma^2 + 83.7779* mu^2* sigma^2 - 83.7779 *sigma^4 - 279.26 *mu* sigma^4 + 293.223* sigma^6)* (log(Xmedian)).^3 + ...
(5.98413* mu^2 - 5.98413* sigma^2 - 41.8889* mu* sigma^2 + 69.8149* sigma^4) *(log(Xmedian)).^4 + ...
(-2.39365* mu + 8.37779* sigma^2) *(log(Xmedian)).^5 + 0.398942 *(log(Xmedian)).^6);
% d7pxdX7=1/(sigma^15* Xmedian^8)*pexp* (0.398942* mu^7 - 11.1704 *mu^6* sigma^2 + mu^2* sigma^6 *(-502.667 + 4691.56* sigma^2 - 5238.91* sigma^4) + ...
% mu^5 *(-8.37779 *sigma^2 + 128.459 *sigma^4) + sigma^8* (167.556 - 2345.78* sigma^2 + 5238.91* sigma^4 - 2010.67 *sigma^6) + mu^4* (167.556* sigma^4 - 781.927* sigma^6) + ...
% mu* sigma^6* (-41.8889 + 1926.89 *sigma^2 - 8101.32 *sigma^4 + 5213.38* sigma^6) + mu^3* (41.8889 *sigma^4 - 1284.59 *sigma^6 + ...
% 2700.44 *sigma^8) + (-2.7926* mu^6 + 67.0223* mu^5* sigma^2 + mu^4 *(41.8889* sigma^2 - 642.297* sigma^4) + ...
% mu *sigma^6 *(1005.33 - 9383.12 *sigma^2 + 10477.8 *sigma^4) + sigma^6 *(41.8889 - 1926.89* sigma^2 + 8101.32* sigma^4 - 5213.38* sigma^6) + ...
% mu^3 *(-670.223* sigma^4 + 3127.71 *sigma^6) + mu^2 *(-125.667* sigma^4 + 3853.78 *sigma^6 - 8101.32 *sigma^8)) *log(Xmedian) + ...
% (8.37779 *mu^5 - 167.556* mu^4 *sigma^2 + sigma^6* (-502.667 + 4691.56* sigma^2 - 5238.91* sigma^4) + mu^3* (-83.7779 *sigma^2 + 1284.59 *sigma^4) + ...
% mu^2* (1005.33* sigma^4 - 4691.56* sigma^6) + mu *(125.667* sigma^4 - 3853.78* sigma^6 + 8101.32 *sigma^8))* (log(Xmedian))^2 + ...
% (-13.963* mu^4 + 223.408* mu^3 *sigma^2 - 41.8889 *sigma^4 + 1284.59 *sigma^6 - 2700.44* sigma^8 + mu^2* (83.7779 *sigma^2 - 1284.59 *sigma^4) + ...
% mu* (-670.223* sigma^4 + 3127.71 *sigma^6))* (log(Xmedian)).^3 + ...
% (13.963 *mu^3 - 41.8889* mu *sigma^2 - 167.556* mu^2* sigma^2 + 167.556 *sigma^4 + 642.297 *mu *sigma^4 - 781.927 *sigma^6)* (log(Xmedian)).^4 + ...
% (-8.37779* mu^2 + 8.37779* sigma^2 + 67.0223* mu *sigma^2 - 128.459* sigma^4)* (log(Xmedian)).^5 +...
% (2.7926* mu - 11.1704 *sigma^2)* (log(Xmedian)).^6 - 0.398942* (log(Xmedian)).^7);
%Now we take median of standard normal density and find its point density
%and derivatives at median
Zm=0;
pz=exp(-0.5* (Zm.^2))/sqrt(2*pi);
dpzdZ=-Zm*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d2pzdZ2=(Zm^2-1)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d3pzdZ3=-(Zm^3-3*Zm)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d4pzdZ4=(Zm^4-6*Zm^2+3)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d5pzdZ5=-(Zm^5-10*Zm^3+15*Zm)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d6pzdZ6=(Zm^6-15*Zm^4+45*Zm.^2-15)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
%d7pzdZ7=-(Zm^7-21*Zm^5+105*Zm.^3-105*Zm)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
%d8pzdZ8=(Zm^8-28*Zm^6+210*Zm.^4-420*Zm.^2+105)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
%d9pzdZ9=-(Zm^9-36*Zm^7+378*Zm.^5-1260*Zm.^3+945*Zm)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
%From relationships between densities through change of density derivative,
%we determine the following derivatives. See post 1492 on this thread for
%details.
dXdZ=pz/px;
d2XdZ2=(dpzdZ-(dpxdX *dXdZ^2))/px
d3XdZ3=(d2pzdZ2- (2* dpxdX* dXdZ *d2XdZ2+dXdZ *(dXdZ^2* d2pxdX2+dpxdX *d2XdZ2)))/px
d4XdZ4=(d3pzdZ3-(3* d2XdZ2* (dXdZ^2* d2pxdX2+dpxdX *d2XdZ2)+3* dpxdX *dXdZ *d3XdZ3+dXdZ* (3 *dXdZ* d2pxdX2* d2XdZ2+dXdZ^3* d3pxdX3+dpxdX* d3XdZ3)))/px
d5XdZ5=(d4pzdZ4- (6 *(dXdZ^2 *d2pxdX2+dpxdX *d2XdZ2) *d3XdZ3+4 *d2XdZ2* (3* dXdZ *d2pxdX2 *d2XdZ2+dXdZ^3* d3pxdX3+dpxdX* d3XdZ3)+4* dpxdX* dXdZ* d4XdZ4+dXdZ *(3* d2pxdX2* d2XdZ2^2+6* dXdZ^2* d2XdZ2* d3pxdX3+4 *dXdZ* d2pxdX2* d3XdZ3+dXdZ^4* d4pxdX4+dpxdX* d4XdZ4)))/px
d6XdZ6=(d5pzdZ5- (10* d3XdZ3 *(3 *dXdZ* d2pxdX2* d2XdZ2+dXdZ^3* d3pxdX3+dpxdX *d3XdZ3)+10* (dXdZ^2* d2pxdX2+dpxdX *d2XdZ2)* d4XdZ4+5* d2XdZ2* (3* d2pxdX2* d2XdZ2^2+6* dXdZ^2* d2XdZ2* d3pxdX3+4 *dXdZ* d2pxdX2* d3XdZ3+dXdZ^4* d4pxdX4+dpxdX* d4XdZ4)+5 *dpxdX* dXdZ* d5XdZ5+dXdZ *(15* dXdZ* d2XdZ2^2* d3pxdX3+10* d2pxdX2* d2XdZ2* d3XdZ3+10* dXdZ^2 *d3pxdX3* d3XdZ3+10* dXdZ^3* d2XdZ2* d4pxdX4+5* dXdZ* d2pxdX2* d4XdZ4+dXdZ^5* d5pxdX5+dpxdX* d5XdZ5)))/px
d7XdZ7=(d6pzdZ6- (20 *(3* dXdZ* d2pxdX2* d2XdZ2+dXdZ^3 *d3pxdX3+dpxdX *d3XdZ3)* d4XdZ4+15* d3XdZ3* (3 *d2pxdX2 *d2XdZ2^2+6* dXdZ^2* d2XdZ2* d3pxdX3+4 *dXdZ* d2pxdX2* d3XdZ3+dXdZ^4* d4pxdX4+dpxdX* d4XdZ4)+15* (dXdZ^2* d2pxdX2+dpxdX* d2XdZ2)* d5XdZ5+6* d2XdZ2* (15* dXdZ* d2XdZ2^2* d3pxdX3+10* d2pxdX2* d2XdZ2 *d3XdZ3+10 *dXdZ^2 *d3pxdX3 *d3XdZ3+10 *dXdZ^3 *d2XdZ2* d4pxdX4+5 *dXdZ *d2pxdX2 *d4XdZ4+dXdZ^5 *d5pxdX5+dpxdX* d5XdZ5)+6* dpxdX* dXdZ *d6XdZ6+dXdZ* (15 *d2XdZ2^3 *d3pxdX3+60 *dXdZ *d2XdZ2 *d3pxdX3* d3XdZ3+10 *d2pxdX2 *d3XdZ3^2+45 *dXdZ^2* d2XdZ2^2 *d4pxdX4+20* dXdZ^3 *d3XdZ3* d4pxdX4+15 *d2pxdX2* d2XdZ2 *d4XdZ4+15 *dXdZ^2* d3pxdX3 *d4XdZ4+15 *dXdZ^4* d2XdZ2* d5pxdX5+6 *dXdZ *d2pxdX2* d5XdZ5+dXdZ^6* d6pxdX6+dpxdX *d6XdZ6)))/px
%d8XdZ8=(d7pzdZ7- (35 *d4XdZ4 *(3 *d2pxdX2* d2XdZ2^2+6* dXdZ^2 *d2XdZ2 *d3pxdX3+4* dXdZ* d2pxdX2* d3XdZ3+dXdZ^4 *d4pxdX4+dpxdX *d4XdZ4)+35 *(3* dXdZ *d2pxdX2 *d2XdZ2+dXdZ^3 *d3pxdX3+dpxdX* d3XdZ3)* d5XdZ5+21* d3XdZ3 *(15 *dXdZ* d2XdZ2^2 *d3pxdX3+10 *d2pxdX2 *d2XdZ2 *d3XdZ3+10* dXdZ^2 *d3pxdX3 *d3XdZ3+10 *dXdZ^3 *d2XdZ2 *d4pxdX4+5* dXdZ *d2pxdX2* d4XdZ4+dXdZ^5 *d5pxdX5+dpxdX *d5XdZ5)+21 *(dXdZ^2 *d2pxdX2+dpxdX* d2XdZ2) *d6XdZ6+7 *d2XdZ2* (15 *d2XdZ2^3 *d3pxdX3+60 *dXdZ *d2XdZ2 *d3pxdX3 *d3XdZ3+10 *d2pxdX2* d3XdZ3^2+45* dXdZ^2 *d2XdZ2^2 *d4pxdX4+20* dXdZ^3* d3XdZ3 *d4pxdX4+15 *d2pxdX2 *d2XdZ2 *d4XdZ4+15 *dXdZ^2 *d3pxdX3 *d4XdZ4+15* dXdZ^4 *d2XdZ2 *d5pxdX5+6 *dXdZ* d2pxdX2 *d5XdZ5+dXdZ^6 *d6pxdX6+dpxdX *d6XdZ6)+7 *dpxdX* dXdZ* d7XdZ7+dXdZ *(105 *d2XdZ2^2 *d3pxdX3 *d3XdZ3+70 *dXdZ* d3pxdX3 *d3XdZ3^2+105 *dXdZ *d2XdZ2^3 *d4pxdX4+210* dXdZ^2 *d2XdZ2 *d3XdZ3 *d4pxdX4+105* dXdZ *d2XdZ2 *d3pxdX3 *d4XdZ4+35 *d2pxdX2* d3XdZ3 *d4XdZ4+35 *dXdZ^3* d4pxdX4 *d4XdZ4+105* dXdZ^3 *d2XdZ2^2* d5pxdX5+35* dXdZ^4* d3XdZ3 *d5pxdX5+21 *d2pxdX2 *d2XdZ2 *d5XdZ5+21 *dXdZ^2 *d3pxdX3 *d5XdZ5+21 *dXdZ^5 *d2XdZ2 *d6pxdX6+7 *dXdZ* d2pxdX2* d6XdZ6+dXdZ^7* d7pxdX7+dpxdX *d7XdZ7)))/px
%d9XdZ9=(d8pzdZ8- (70* (3 *d2pxdX2* d2XdZ2^2+6* dXdZ^2* d2XdZ2* d3pxdX3+4* dXdZ* d2pxdX2* d3XdZ3+dXdZ^4 *d4pxdX4+dpxdX* d4XdZ4) *d5XdZ5+56* d4XdZ4* (15* dXdZ* d2XdZ2^2 *d3pxdX3+10 *d2pxdX2 *d2XdZ2 *d3XdZ3+10 *dXdZ^2 *d3pxdX3 *d3XdZ3+10 *dXdZ^3 *d2XdZ2 *d4pxdX4+5 *dXdZ *d2pxdX2* d4XdZ4+dXdZ^5 *d5pxdX5+dpxdX *d5XdZ5)+56 *(3 *dXdZ *d2pxdX2 *d2XdZ2+dXdZ^3* d3pxdX3+dpxdX *d3XdZ3) *d6XdZ6+28* d3XdZ3 *(15 *d2XdZ2^3 *d3pxdX3+60 *dXdZ *d2XdZ2 *d3pxdX3 *d3XdZ3+10 *d2pxdX2* d3XdZ3^2+45 *dXdZ^2 *d2XdZ2^2 *d4pxdX4+20* dXdZ^3 *d3XdZ3 *d4pxdX4+15* d2pxdX2 *d2XdZ2 *d4XdZ4+15 *dXdZ^2 *d3pxdX3 *d4XdZ4+15 *dXdZ^4 *d2XdZ2 *d5pxdX5+6 *dXdZ *d2pxdX2* d5XdZ5+dXdZ^6 *d6pxdX6+dpxdX *d6XdZ6)+28* (dXdZ^2 *d2pxdX2+dpxdX *d2XdZ2) *d7XdZ7+8 *d2XdZ2 *(105 *d2XdZ2^2 *d3pxdX3 *d3XdZ3+70* dXdZ* d3pxdX3 *d3XdZ3^2+105 *dXdZ* d2XdZ2^3 *d4pxdX4+210 *dXdZ^2 *d2XdZ2* d3XdZ3* d4pxdX4+105 *dXdZ* d2XdZ2 *d3pxdX3 *d4XdZ4+35 *d2pxdX2 *d3XdZ3 *d4XdZ4+35* dXdZ^3 *d4pxdX4 *d4XdZ4+105* dXdZ^3 *d2XdZ2^2 *d5pxdX5+35 *dXdZ^4 *d3XdZ3 *d5pxdX5+21* d2pxdX2* d2XdZ2 *d5XdZ5+21 *dXdZ^2 *d3pxdX3* d5XdZ5+21* dXdZ^5* d2XdZ2 *d6pxdX6+7 *dXdZ* d2pxdX2 *d6XdZ6+dXdZ^7 *d7pxdX7+dpxdX *d7XdZ7)+8 *dpxdX* dXdZ* d8XdZ8+dXdZ *(280* d2XdZ2 *d3pxdX3 *d3XdZ3^2+105 *d2XdZ2^4 *d4pxdX4+840 *dXdZ* d2XdZ2^2 *d3XdZ3 *d4pxdX4+280 *dXdZ^2 *d3XdZ3^2 *d4pxdX4+210 *d2XdZ2^2 *d3pxdX3 *d4XdZ4+280 *dXdZ* d3pxdX3 *d3XdZ3 *d4XdZ4+420 *dXdZ^2 *d2XdZ2 *d4pxdX4 *d4XdZ4+35* d2pxdX2 *d4XdZ4^2+420 *dXdZ^2 *d2XdZ2^3* d5pxdX5+560* dXdZ^3 *d2XdZ2* d3XdZ3* d5pxdX5+70* dXdZ^4 *d4XdZ4 *d5pxdX5+168 *dXdZ *d2XdZ2 *d3pxdX3 *d5XdZ5+56 *d2pxdX2 *d3XdZ3 *d5XdZ5+56* dXdZ^3 *d4pxdX4 *d5XdZ5+210 *dXdZ^4 *d2XdZ2^2 *d6pxdX6+56 *dXdZ^5 *d3XdZ3 *d6pxdX6+28 *d2pxdX2 *d2XdZ2 *d6XdZ6+28* dXdZ^2* d3pxdX3 *d6XdZ6+28 *dXdZ^6 *d2XdZ2 *d7pxdX7+8 *dXdZ *d2pxdX2 *d7XdZ7+dXdZ^8 *d8pxdX8+dpxdX *d8XdZ8)))/px
%d10XdZ10=(d9pzdZ9- (126*d5XdZ5*(15*dXdZ*d2XdZ2^2*d3pxdX3+10*d2pxdX2*d2XdZ2*d3XdZ3+10*dXdZ^2*d3pxdX3*d3XdZ3+10*dXdZ^3*d2XdZ2*d4pxdX4+5*dXdZ*d2pxdX2*d4XdZ4+dXdZ^5*d5pxdX5+dpxdX*d5XdZ5)+126*(3*d2pxdX2*d2XdZ2^2+6*dXdZ^2*d2XdZ2*d3pxdX3+4*dXdZ*d2pxdX2*d3XdZ3+dXdZ^4*d4pxdX4+dpxdX*d4XdZ4)*d6XdZ6+84*d4XdZ4*(15*d2XdZ2^3*d3pxdX3+60*dXdZ*d2XdZ2*d3pxdX3*d3XdZ3+10*d2pxdX2*d3XdZ3^2+45*dXdZ^2*d2XdZ2^2*d4pxdX4+20*dXdZ^3*d3XdZ3*d4pxdX4+15*d2pxdX2*d2XdZ2*d4XdZ4+15*dXdZ^2*d3pxdX3*d4XdZ4+15*dXdZ^4*d2XdZ2*d5pxdX5+6*dXdZ*d2pxdX2*d5XdZ5+dXdZ^6*d6pxdX6+dpxdX*d6XdZ6)+84*(3*dXdZ*d2pxdX2*d2XdZ2+dXdZ^3*d3pxdX3+dpxdX*d3XdZ3)*d7XdZ7+36*d3XdZ3*(105*d2XdZ2^2*d3pxdX3*d3XdZ3+70*dXdZ*d3pxdX3*d3XdZ3^2+105*dXdZ*d2XdZ2^3*d4pxdX4+210*dXdZ^2*d2XdZ2*d3XdZ3*d4pxdX4+105*dXdZ*d2XdZ2*d3pxdX3*d4XdZ4+35*d2pxdX2*d3XdZ3*d4XdZ4+35*dXdZ^3*d4pxdX4*d4XdZ4+105*dXdZ^3*d2XdZ2^2*d5pxdX5+35*dXdZ^4*d3XdZ3*d5pxdX5+21*d2pxdX2*d2XdZ2*d5XdZ5+21*dXdZ^2*d3pxdX3*d5XdZ5+21*dXdZ^5*d2XdZ2*d6pxdX6+7*dXdZ*d2pxdX2*d6XdZ6+dXdZ^7*d7pxdX7+dpxdX*d7XdZ7)+36*(dXdZ^2*d2pxdX2+dpxdX*d2XdZ2)*d8XdZ8+9*d2XdZ2*(280*d2XdZ2*d3pxdX3*d3XdZ3^2+105*d2XdZ2^4*d4pxdX4+840*dXdZ*d2XdZ2^2*d3XdZ3*d4pxdX4+280*dXdZ^2*d3XdZ3^2*d4pxdX4+210*d2XdZ2^2*d3pxdX3*d4XdZ4+280*dXdZ*d3pxdX3*d3XdZ3*d4XdZ4+420*dXdZ^2*d2XdZ2*d4pxdX4*d4XdZ4+35*d2pxdX2*d4XdZ4^2+420*dXdZ^2*d2XdZ2^3*d5pxdX5+560*dXdZ^3*d2XdZ2*d3XdZ3*d5pxdX5+70*dXdZ^4*d4XdZ4*d5pxdX5+168*dXdZ*d2XdZ2*d3pxdX3*d5XdZ5+56*d2pxdX2*d3XdZ3*d5XdZ5+56*dXdZ^3*d4pxdX4*d5XdZ5+210*dXdZ^4*d2XdZ2^2*d6pxdX6+56*dXdZ^5*d3XdZ3*d6pxdX6+28*d2pxdX2*d2XdZ2*d6XdZ6+28*dXdZ^2*d3pxdX3*d6XdZ6+28*dXdZ^6*d2XdZ2*d7pxdX7+8*dXdZ*d2pxdX2*d7XdZ7+dXdZ^8*d8pxdX8+dpxdX*d8XdZ8)+9*dpxdX*dXdZ*d9XdZ9+dXdZ*(280*d3pxdX3*d3XdZ3^3+1260*d2XdZ2^3*d3XdZ3*d4pxdX4+2520*dXdZ*d2XdZ2*d3XdZ3^2*d4pxdX4+1260*d2XdZ2*d3pxdX3*d3XdZ3*d4XdZ4+1890*dXdZ*d2XdZ2^2*d4pxdX4*d4XdZ4+1260*dXdZ^2*d3XdZ3*d4pxdX4*d4XdZ4+315*dXdZ*d3pxdX3*d4XdZ4^2+945*dXdZ*d2XdZ2^4*d5pxdX5+3780*dXdZ^2*d2XdZ2^2*d3XdZ3*d5pxdX5+840*dXdZ^3*d3XdZ3^2*d5pxdX5+1260*dXdZ^3*d2XdZ2*d4XdZ4*d5pxdX5+378*d2XdZ2^2*d3pxdX3*d5XdZ5+504*dXdZ*d3pxdX3*d3XdZ3*d5XdZ5+756*dXdZ^2*d2XdZ2*d4pxdX4*d5XdZ5+126*d2pxdX2*d4XdZ4*d5XdZ5+126*dXdZ^4*d5pxdX5*d5XdZ5+1260*dXdZ^3*d2XdZ2^3*d6pxdX6+1260*dXdZ^4*d2XdZ2*d3XdZ3*d6pxdX6+126*dXdZ^5*d4XdZ4*d6pxdX6+252*dXdZ*d2XdZ2*d3pxdX3*d6XdZ6+84*d2pxdX2*d3XdZ3*d6XdZ6+84*dXdZ^3*d4pxdX4*d6XdZ6+378*dXdZ^5*d2XdZ2^2*d7pxdX7+84*dXdZ^6*d3XdZ3*d7pxdX7+36*d2pxdX2*d2XdZ2*d7XdZ7+36*dXdZ^2*d3pxdX3*d7XdZ7+36*dXdZ^7*d2XdZ2*d8pxdX8+9*dXdZ*d2pxdX2*d8XdZ8+dXdZ^9*d9pxdX9+dpxdX*d9XdZ9)))/px
%Since we have found the derivatives of lognormal density variable with
%respect to standard normal variable, we can convert them as coefficients
%of Z-series.
c0=Xmedian
c(1)=dXdZ
c(2)=1/2*1*d2XdZ2
c(3)=1/6*1*d3XdZ3
c(4)=1/24*1*d4XdZ4
c(5)=1/120*1*d5XdZ5
c(6)=1/720*1*d6XdZ6
c(7)=1/720/7*1*d7XdZ7
%Below construct lognormal density from Z-series coefficients
%First find discretization of standard normal
dNn=.1/2; % Normal density subdivisions width. would change with number of subdivisions
Nn=45*4; % No of normal density subdivisions
NnMid=((1+Nn)/2)*dNn;
Z(1:Nn)=(((1:Nn)*dNn)-NnMid);
%Now construct lognormal variable from its Z-series coefficients
X(1:Nn)=c0;
for nn=1:7
X(1:Nn)=X(1:Nn)+c(nn)*Z(1:Nn).^nn;
end
%Find change of density derivative with respect to normal for construction
%of density
DfX(1:Nn)=0;
for nn=2:Nn-1
DfX(nn) = (X(nn + 1) - X(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
DfX(Nn)=DfX(Nn-1);
DfX(1)=DfX(2);
%Find lognormal density at the normal grid
pX(1:Nn)=0;
for nn = 1:Nn
pX(nn) = (normpdf(Z(nn),0, 1))/abs(DfX(nn));
end
%Below directly find lognormal density from analytic formula
Y(1:15000)=.00001+.01*(1:15000);
pY=1./(sqrt(2* pi)*sigma.* Y).* exp(-.5.* (log(Y)-mu).^2./sigma^2);
%Compare the graph of same lognormal density derived two different ways from
%Z-series coefficients --red and from analytic formula as green.
plot(X(1:Nn),pX(1:Nn),'r',Y(1:15000),pY(1:15000),'g');
title('Comparison of Lognormal Density from Z-series VS from analytic closed form formula');%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
legend({'Lognormal density From Z-series','Lognormal density from analytic formula of probability density'},'Location','northeast')
str=input('red line is lognormal density from Z-series, green is analytic formula.');
end
```

.

.

Here is the graph comparison from the above program.

Here is the comparison of tails of the lognormal density from Z-series to exact density. Please note that you can arbitrarily improve the tail by increasing the order of Z-series which only requires simple algebra as opposed to cumbersome optimizations that we have to do when we find the Z-series from moments of random variables.