Serving the Quantitative Finance Community

 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 2nd, 2024, 9:52 pm

Friends, I have decided that I will go ahead with this and first try to represent HO(Hydrogen Orbitals) and STO(Slater Type Orbitals) in the form of one common set of orthogonal polynomials (with one underlying density) and try to represent all of the different HO's and STO's in the form of a series in appropriate associated Laguerre or Laguerre polynomials(all with one underlying density). If it works well, I will try to see if we can do all sort of analytics with expectations of the appropriate generalized gamma random variables. 
I will try to write matlab programs for this and also post my matlab programs on the forum when they start to work well. I will start tomorrow but it can take several days before the first programs are ready.
For friends in mathematics who do not know quantum mechanics, linear combinations of HOs(Hydrogen Orbitals) and STOs(Slater Type Orbitals) are used (as basis) to approximate the wave functions of atoms and molecules with larger number of electrons with all sort of electron-electron interactions between them.
Atomic Orbital
Slater Type Orbitals
Linear Combination Of Atomic Orbitals
Hydrogen Like Atom
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 3rd, 2024, 8:40 pm

Friends, I could not do any work today. I also had fluanxol injection today. Fluanxol injections were cruelly and mercilessly drugged by Lundbeck in Pakistan. Till four or five months ago, injections were so bad that it was very difficult for me to even express myself properly to anyone for fifteen days after the injections because they would sap all neurotransmitters from the brain. And I would have all sorts of pains and aches in the body. I think a lot of good people protested and crooks at Lundbeck stopped drugging the injections. But recently crooks at Lundbeck started drugging the injections again and injections with manufacturing date of 08/23 are again being drugged. I think injections with many earlier dates are better but Lundbeck has very heavily distributed fluanxol injections with manufacturing date of 08/23 and more than 95% of the pharmacies have injections with this most recent manufacturing date. I tried to buy injection with an earlier date and it was extremely difficult but I was able to get one injection. I hope that injection I had would be better. 
Anyway, I hope to start experimenting tomorrow how we can faithfully represent functions/polynomials of one density in terms of functions/polynomials of the other density. We will try to see for which densities we could do that and for which densities we cannot do that. And then we would turn to one common basis for various Slater orbitals.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 4th, 2024, 8:21 am

Friends, I had already distributed a program in post 1837. This program constructs a series representation in terms of Laguerre Polynomials with gamma density as the base density in Laguerre Polynomials. Construction of this Laguerre density is done by moment matching with respect to input moments of the random variable. Then we convert this Laguerre series representation(with underlying gamma density) into Hermite series representation where underlying density is standard normal. 
As I recall that part about construction of Laguerre series density from exogenously specified moments was rather approximate but once the Laguerre Series density found and its coefficients determined, the part about converting the Laguerre Series to Hermite series coefficients was quite faithful and worked very well.
There was another time when I tried to convert functions of Hermite polynomials into Functions of Legendre polynomial (with underlying uniform density) but I had not thought clearly about the problem and I was being very goofy since uniform density does not take any derivatives and first derivative is zero so it was impossible to match any derivatives and I realized that conversion across Hermite and Legendre could not be done.
Matlab Program here about Conversion from Laguerre Series to Hermite Series in Post 1837: https://forum.wilmott.com/viewtopic.php?t=99702&start=1830#p875723
Please also read post 1838 and earlier posts (in conversion from Lognormal to Hermite Context) 1636, 1637, 1638. Every (polynomial or closed form) function of a lognormal density that takes continuous derivatives can be represented by a Hermite Series with underlying standard normal density but for functions that take large powers, hermite series expansion might be needed to a sufficiently large order.

Now in terms of quantum mechanics context, I will try to see if we can shift a series in one generalized gamma density into a different generalized gamma density(with completely general parameters). It will just find a plain polynomial series but since span of all the different polynomials of same order is usually the complete span so we can easily convert the plain polynomial series into an associated Laguerre Series (with respect to generalized gamma density which is the right inner product weight for the orthogonal associated Laguerre density) . I am working on it today and hope that my program will be ready in two days.
I am not sure how well it goes but just being curious I will also try to convert associated Laguerre series (generalized gamma densities) into Hermite series with underlying standard normal.

  
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 4th, 2024, 11:13 am

Friends, I have been thinking that there are so many densities that do not have any orthogonal polynomials and if there could be some possible way to construct orthogonal polynomials easily. And I had some basic thoughts about it that I want to share with friends.
It seems that we could possibly construct orthogonal polynomials for a large number of densities by solving a set of many linear and one non-linear equations recursively for each increasing polynomial.
Zeroth polynomial is chosen as a constant.
P0(X)= 1
The first polynomial can be  P1(X)=a1 X + a0 where a0 and a1 have to be found from equations
For P0 and P1 to be orthogonal the product of P0(X) and P1(X) when evaluated with respect to specified density should be zero.

E[P0(X) P1(X)]=a0 +  a1 E[X]=0;
and we can have one non-linear equation for normality that 
E[P1(X) P1(X)]=a0^2 + 2 a0 a1 E[X] + a1^2 E[X^2] = 1;

Once first polynomial is fixed, the next polynomial is P2(X)=b0+ b1 X + b2 X^2
For orthogonality, the product of second polynomial with zeroth polynomial and first polynomial with respect to underlying density should be zero. This product is

E[P0(X) P2(X)]= b0 + b1 E[X] + b2 E[X^2] =0
E[P1(X) P2(X)]= a0 b0 + a0 b1 E[X] + a0 b2 E[X^2] + a1 b0 E[X] + a1 b1 E[X^2] + a1 b2 E[X^3]=0
and one non-linear equation for normality given as
E[P2(X)^2] = b0^2 + 2 b0 b1 E[X] + 2 b0 b2 E[X^2]   + b1^2 E[X^2] +2 b1 b2 E[X^3] + b2^2 E[X^4] =0

This will give us three equations for three unknowns.

We repeat the procedure for third polynomial as
Now P3(X)= c0+ c1 X + c2 X^2 + c3 X^3
For orthogonality of this third polynomial to hold, expectation of its product with zeroth, first and second polynomials should be zero.
E[P0(X) P3(X)]= c0 + c1 E[X] + c2 E[X^2] + c3 E[X]^3 =0
E[P1(X) P3(X)]= a0 c0 + a0 c1 E[X] + a0 c2 E[X^2] + a0 c3 E[X^3] + a1 c0 E[X] + a1 c1 E[X^2] + a1 c2 E[X^3]+ a1 c3 E[X^4]=0
E[P2(X) P3(X)]= b0 c0 + b0 c1 E[X] + b0 c2 E[X^2] + b0 c3 E[X^3] + b1 c0 E[X] + b1 c1 E[X^2] + b1 c2 E[X^3]+ b1 c3 E[X^4]+ b2 c0 E[X^2] + b2 c1 E[X^3] + b2 c2 E[X^4]+ b2 c3 E[X^5]=0
and we would have to solve one non-linear equation for the product of third polynomial with itself and equate it to one for normality
E[P3(X) P3(X)]= c0^2 + 2 c0 c1 E[X] + 2 c0 c2 E[X^2] +2 c0 c3 E[X^3]  + c1^2 E[X^2] +2 c1 c2 E[X^3]+2 c1 c3 E[X^4] + c2^2 E[X^4] + c2 c3 E[X^5]+ c3^2 E[X^6]=0

We can recursively continue this process to find higher order orthonormalized polynomials 

I think most of the times the above equations should be solvable though I am sure for some densities they might not be solvable. But this seems like a simple common sense method of constructing orthonormalized polynomials for a given density.

Many times we would like to omit the normalizing step(with non-linear equation of coefficients) and try to impose the condition that coefficient on highest order power on every polynomial is one. This way we could match coefficients found from our method with coefficients of most of the known orthogonal polynomials.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 4th, 2024, 7:44 pm

Friends, I am writing a small post trying to explain the basics of how to convert functions and polynomials of one density with continuous derivatives into polynomials of another density with continuous derivatives. I will not go deep into mathematical calculations but try to give idea about logic and methods involved in such change of probability distributions.
Suppose we have one standard normal density and we write it as

[$]p_Z(z) \, = \, \frac{1}{\sqrt{2 \pi}} \, \exp(-\, \frac{z^2}{2}) [$]

and we have a gamma density and we write it as

[$]p_G(g) \, = \, \frac{1}{\Gamma(k)\, {\theta}^k}\, g^{(k-1)} \, \exp(-\frac{g}{\theta}) [$]

We know the change of variables equation for continuous univariate probability densities as

[$]\, p_G(g) \, |\frac{dg}{dz}|\, = \, p_Z(z) [$]  Eq(1)

We select a common CDF point for both densities where we find derivatives between gamma variable g and standard normal variable z. In my experiments I found that median of both densities was a very good choice. Median of many densities is not known analytically but can be calculated with very high precision using numerical methods if CDF is known. We suppose that pdf of both densities involved is known in closed from and we can analytically take all order of derivatives of each pdf.

From Eq(1) we can find first derivative between both densities at any common CDF point which is median in our special case. For that we just input values of z and g at their respective medians in the respective probability distribution functions in Eq(1) and back out the value of  [$]|\frac{dg}{dz}|\,[$]. We will have to be careful about absolute value but if both densities are increasing at median, [$]\frac{dg}{dz}\,[$] will always be positive and we would not need to worry about absolute sign. It turns out that we can find all higher derivatives [$]\frac{d^{n} g}{dz^n}\,[$] by differentiating the Eq(1) n-1 times on both sides and doing some algebra. Please read posts 1637, 1638, 1639 for details.
viewtopic.php?t=99702&start=1635#p873552
viewtopic.php?t=99702&start=1485#p871049
Alternatively we could also have found  [$]\frac{dz}{dg}\,[$] and all of the nth derivatives [$]\frac{d^{n} z}{dz^n}\,[$] by writing Eq(1) differently as

[$]\, p_G(g) \,  = \, p_Z(z) |\frac{dz}{dg}|\,[$]  Eq(2)

It turns out that differentiating many densities becomes very tedious after 9th or tenth derivatives but an algorithm can easily be written for respective densities that can differentiate them several tens of times for arbitrary precision.

Once we have found derivatives between both densities to sufficient order, we can use them to represent for example functions of standard normal as polynomials in gamma density variable or we could represent functions of gamma density as polynomials in standard normal random variable. Of course, some restrictions would apply but the conversions can be useful in a large number of settings.

Let us suppose we have a function or a polynomial in Gamma density variable written as [$]f(g)[$] and we want to represent it as polynomial function h(z) of standard normal variable z. Again we would do this polynomial expansion around common CDF points which we prefer to choose as median of both densities. In order to do that , we would have to write 

[$]h(z)=f(g(z))[$] Eq(3)

We will expand the right hand side of Eq(3) in a taylor series and it is given to fourth order as (please note that derivatives between gamma random variable g and normal random variable z of all orders [$]\frac{d^n g}{dz^n}[$] have previously been calculated at a common CDF point possibly median either from differentiating Eq(1) or Eq(2) on both sides repeatedly and using some algebra) Here [$]z_0[$] is median of standard normal density.

[$]h(z)\, = \, f(g(z)) \, = \, f(g(z_0))+(z-z_0) \, \frac{dg}{dz} \, f'(g(z_0)) \\
+\frac{1}{2} (z-z_0)^2 \big[ \, {(\frac{dg}{dz})}^2 f''(g(z_0))+ \, \frac{d^2g}{dz^2} \, f'(g(z_0))\big] \\
+\frac{1}{6} (z-z_0)^3 \big[f^{(3)}(g(z_0))  \, {( \frac{dg}{dz})}^3 \, +3  \, \frac{dg}{dz} \,  \, \frac{d^2g}{dz^2} \, f''(g(z_0))+ \, \frac{d^3g}{dz^3} \, f'(g(z_0))\big]\\
+\frac{1}{24} (z-z_0)^4 \big[f^{(4)}(g(z_0))  \, {(\frac{dg}{dz})}^4 \,+6 f^{(3)}(g(z_0)) \, {(\frac{dg}{dz})}^2 \, \, \frac{d^2g}{dz^2} \,+3 \, {(\frac{d^2g}{dz^2})}^2 \, f''(g(z_0))+4 \, \frac{d^3g}{dz^3} \, \, \frac{dg}{dz} \, f''(g(z_0))+\, {(\frac{d^4g}{dz^4})} \, f'(g(z_0))\big]\\
+O\left((z-z_0)^5\right)[$]

We can write algorithm for an automated program that expands the above series to arbitrary order as opposed to fourth order we have done in previous equation.
Similarly we could have just as easily converted many functions of standard normal into functions of gamma density using the same logic.

Once we have found a polynomial representation in variable of a particular density, we could convert it easily into a series in orthogonal polynomials of the same density. This is because span of all nth order polynomials is the same and everything that can be represented by an nth order polynomial, can also be represnted as a series in first n orthogonal polynomials and it just requires some algebra (we have always been converting all nth order polynomial functions of standard normal in the form of hermite polynomials series with first n polynomials).
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 6th, 2024, 10:41 am

Friends, I am writing a simple matlab program that will let friends expand normal, generalized gamma and lognormal type densities from original variable density to expansion in other densities. I have worked out the normal and generalized gamma part. 
Out of curiosity, I tried to create a standard normal density from an expansion in gamma density with parameters k=9, theta=.5 and did an expansion in first nine terms. The approximation of standard normal density was faithful but gamma-approximated tail decreased very slightly in probability mass between 2 SD and 3 SD but then increased again and then continued forever with the earlier mass defect with extremely small probability.
I will post these program later tonight or tomorrow and you could expand all sort of densities across each other with original density overlaid in different color.

Here standard normal density approximation with expansion in terms of  Gamma density gamma(9,.5) in the body. Green line in graph is true standard normal with analytical formula. Red line is from a series approximation of standard normal density in terms of gamma(9,.5) density variable.

This first graph concentrates on body of standard normal density from -5 to +5 which was range of green analytical density.
Image

This second graph shows all the tail with miniscule probability without concentrating on the body of the graph.


Image
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 6th, 2024, 11:43 am

Friends, I will post the new density programs tonight or later tomorrow. For a day or two, we will play with representation of various densities in terms of other densities. We will also try to work on various Slater Type orbitals as a series in a single type of orthogonal polynomials with a common base density and if that part works well, we will move on to attempts at solution of Stationary N-electron Schordinger Equation with our common orthogonal polynomials and try to find its solution.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 7th, 2024, 9:16 pm

Friends, I have worked out the program and written explanatory comments on main function. You can expand three of the standard normal, three parameter generalized gamma and two parameter lognormal density in terms of  either of the same three densities (with different values of parameters).  I want to check the program a bit more tomorrow therefore not posting it now. But will post it early tomorrow before friends start their new week in their offices. 
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 8th, 2024, 9:09 am

Friends, here is the new matlab program to convert random variables with one density in terms of an expansion in random variable of another different density. There are three input choices of random variables to be expanded. 1. Standard normal, 2. Generalized gamma(a2,d2,p2) and 3. Lognormal(mu2,sigma2).
One of these is also the choice of random variable in which expansion is done. 1. Standard normal, 2. Generalized gamma(a,d,p) and 3. Lognormal(mu,sigma).
We naturally assume that  a2, d2 and p2 are not the same as a, d, and p.

Here is the matlab program.
function [] = ConvertGenGammaDensityToOtherDensities()

%This program allows you to represent one of the three standard normal density, 
%generalized gamma density and lognormal density in terms of either of another 
%standard normal, generalized gamma density or lognormal density. You can
%for example choose generalized gamma density with certain parameters as
%input density to be modelled and expand in another generalized gamma
%density with different parameters. You can also try to expand normal
%density in terms of either of the two densities and so on. 
%In this program functions of the random variables associaated with 
%densities are not expanded. This program only expands densities with 
%respect to each other 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Input Densities Block. These are the densities that are to be modelled by
%expansions in another set of densities.There are three possibilities here
%1. Standard Normal density 2. Generalized gamma density
%3. Lognormal density
%Decomment appropriate density out of the three in this block


%Number 1. Uncomment next two lines in this block for Standard normal Density
%x2Median=0;
%[p2x2,dp2dx2] = CalculateStandardNormalDensityDerivatives(x2Median)
%VCase1=1;
% Nn2=201;
% dNn2=(Z2High-Z2low)/(Nn2-1);
% Z2(1:Nn2)=L2low+(0:Nn2-1)*dNn2;
% NormalPdf2(1:Nn2)=1/sqrt(2*pi).*exp(-.5*Z2(1:Nn2).^2);




%Number 2. Uncomment Below for generalized gamma density. 
%Feel free to change generalized gamma parameters
d2=2;
p2=1; %Keep it 1 for gamma density
a2=1;
[x2Median] = CalculateMedianOfGenGammaBisection(d2,p2,a2)
[p2x2,dp2dx2] = CalculateGeneralizedGammaDensityDerivatives(x2Median,d2,p2,a2)
VCase1=2;
%Below we analytically calculate the above Generalized Gamma density and 
%it would be overlaid on expanded density in green to represent the true density
[G2High] = InvertCDFOfGenGammaBisection(.99999,d2,p2,a2);
[G2low] = InvertCDFOfGenGammaBisection(.000001,d2,p2,a2); 

Nn2=200001;%Sometimes due to huge tail, this number has to be very large.
dNn2=(G2High-G2low)/(Nn2-1);
G2(1:Nn2)=G2low+(0:Nn2-1)*dNn2;
GenGammaPdf2(1:Nn2)=(p2/a2^d2)*G2(1:Nn2).^(d2-1).*exp(-(G2(1:Nn2)/a2).^p2)/gamma(d2/p2);




%Number 3. Uncomment for lognormal denstiy. 
% mu2=.5;
% sigma2=.4;
% x2Median=exp(mu2);
% [p2x2,dp2dx2] = CalculateLognormalDensityDerivatives(x2Median,mu2,sigma2)
%VCase1=3;
%Below we analytically calculate the above Lognormal density and 
%it would be overlaid on expanded density in green to represent the true density
% [L2High] = InvertCDFOfLognormalBisection(.99999,mu2,sigma2);
% [L2low] =  InvertCDFOfLognormalBisection(.00001,mu2,sigma2);
% 
% Nn2=201;
% dNn2=(L2High-L2low)/(Nn2-1);
% L2(1:Nn2)=L2low+(0:Nn2-1)*dNn2;
% LognormalPdf2(1:Nn2)=1/sqrt(2*pi)./sigma2./L2(1:Nn2).*exp(-.5*(log(L2(1:Nn2)-mu2)/sigma2).^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Expansion Densities Block.These are base densities used for expansion of other 
%variables.
%Uncomment appropriate density in this block.
%This is the density in whose varaible the expansion of the modelled
%density would be done.
%Uncomment Below to represent any density from previous code block in the
%form of variable of the esnity below.
% 

%Number 1. Uncomment next two lines in this block for if the other 
%density is being expanded in terms of standard normal random 
%variable.
x1Median=0;
[p1x1,dp1dx1] = CalculateStandardNormalDensityDerivatives(x1Median)
VCase2=1;

% %Uncomment next few lines in this block for if the other 
% %density is being expanded in terms of generalized gamma random 
% %variable. Also specify the parameters of generalized gamma variable in
% %whose form the other density is being expanded
% d=2;
% p=1;
% a=4;
% [x1Median] = CalculateMedianOfGenGammaBisection(d,p,a)
% [p1x1,dp1dx1] = CalculateGeneralizedGammaDensityDerivatives(x1Median,d,p,a)
% VCase2=2;

% % %Uncomment next few lines if the other density is being expanded in 
% %terms of lognormal random variable. You also have to specify
% %the parameters of lognormal random variable.
%  mu=0;
%  sigma=.3;
%  x1Median=exp(mu);
%  [p1x1,dp1dx1] = CalculateLognormalDensityDerivatives(x1Median,mu,sigma);
%  VCase2=3;
% %  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%The function below takes the values of density (p2x2) to be modelled by expansion
% and the values of derivatives of the density (dp2dx2) calculated at
% median. This is density being expanded.
%The function als takes the values (p1x1) and derivatives (dp1dx1)of the
%density in whose variable the other density would be expanded. This is
%density in which expansion is done.
 
if(((VCase1==1)||(VCase1==2)) && ((VCase2==1) || (VCase2==2)))
    OrderOfExpansion=9;
else
    OrderOfExpansion=7;
end
[dX2dX1] = CalculateDerivsBetweenRVarsFromDensityDerivs02(p1x1,p2x2,dp1dx1,dp2dx2,OrderOfExpansion);
%[dX2dX1b] = CalculateDerivsBetweenRVarsFromDensityDerivs(p1x1,p2x2,dp1dx1,dp2dx2);


%Below are coefficients of the expansion polynomial so that density to be
%modelled could be approximately expanded.
c0=x2Median;
for nn=1:OrderOfExpansion    
    c(nn)=dX2dX1(nn)/factorial(nn);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(VCase2==2)

[GHigh] = InvertCDFOfGenGammaBisection(.99999,d,p,a);
[Glow] = InvertCDFOfGenGammaBisection(.000001,d,p,a); 

Nn=200001;%Sometimes due to huge tail, this number has to be very large.
dNn=(GHigh-Glow)/(Nn-1);
G(1:Nn)=Glow+(0:Nn-1)*dNn;
GenGammaPdf(1:Nn)=(p/a^d)*G(1:Nn).^(d-1).*exp(-(G(1:Nn)/a).^p)/gamma(d/p);

%Now construct random variable from its Z-series coefficients
X(1:Nn)=c0;
for nn=1:OrderOfExpansion
    X(1:Nn)=X(1:Nn)+c(nn)*(G(1:Nn)-x1Median).^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))/(G(nn + 1) - G(nn - 1));
    
    %Change of variable derivative for densities
end
DfX(Nn)=DfX(Nn-1);
%Find lognormal density at the normal grid
pX(1:Nn)=0;
for nn = 1:Nn
    pX(nn) = (GenGammaPdf(nn))/abs(DfX(nn));
end

if(VCase1==1)
plot(X(1:Nn),pX(1:Nn),'r', Z2(1:Nn2),NormalPdf2(1:Nn2),'g');
end
%If Modelled density is Generalized gamma, uncomment the line below.
if(VCase1==2)
plot(X(1:Nn),pX(1:Nn),'r',G2(1:Nn2),GenGammaPdf2(1:Nn2),'g');
end
%If Modelled density is Lognormal, uncomment the line below.
if(VCase1==3)
plot(X(1:Nn),pX(1:Nn),'r',L2(1:Nn2),LognormalPdf2(1:Nn2),'g');

end


str=input("Look at the plot comparison of original density in green and its expansion in terms of Generalized Gamma density in red")

end


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(VCase2==3)


[LHigh] = InvertCDFOfLognormalBisection(.99999,mu,sigma);
[Llow] =  InvertCDFOfLognormalBisection(.00001,mu,sigma);

Nn=201;
dNn=(LHigh-Llow)/(Nn-1);
L(1:Nn)=Llow+(0:Nn-1)*dNn;
LognormalPdf(1:Nn)=1/sqrt(2*pi)./sigma./L(1:Nn).*exp(-.5*(log(L(1:Nn)-mu)/sigma).^2);
X(1:Nn)=c0;
for nn=1:OrderOfExpansion
    X(1:Nn)=X(1:Nn)+c(nn).*(L(1:Nn)-x1Median).^nn;
end

%Find change of density derivative with respect to lognormal for construction
%of density
DfX(1:Nn)=0;
for nn=2:Nn-1
    
    DfX(nn) = (X(nn + 1) - X(nn - 1))/(L(nn + 1) - L(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) = (LognormalPdf(nn))/abs(DfX(nn));
end



%If Modelled density is standard normal, uncomment the line below.
if(VCase1==1)
plot(X(1:Nn),pX(1:Nn),'r', Z2(1:Nn2),NormalPdf2(1:Nn2),'g');
end
%If Modelled density is Generalized gamma, uncomment the line below.
if(VCase1==2)
plot(X(1:Nn),pX(1:Nn),'r',G2(1:Nn2),GenGammaPdf2(1:Nn2),'g');
end
%If Modelled density is Lognormal, uncomment the line below.
if(VCase1==3)
plot(X(1:Nn),pX(1:Nn),'r',L2(1:Nn2),LognormalPdf2(1:Nn2),'g');

end
str=input("Look at the plot comparison of original density in green and its expansion in terms of lognormal density in red")

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(VCase2==1)
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 random 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



if(VCase1==1)
plot(X(1:Nn),pX(1:Nn),'r', Z2(1:Nn2),NormalPdf2(1:Nn2),'g');
end
%If Modelled density is Generalized gamma, uncomment the line below.
if(VCase1==2)
plot(X(1:Nn),pX(1:Nn),'r',G2(1:Nn2),GenGammaPdf2(1:Nn2),'g');
end
%If Modelled density is Lognormal, uncomment the line below.
if(VCase1==3)
plot(X(1:Nn),pX(1:Nn),'r',L2(1:Nn2),LognormalPdf2(1:Nn2),'g');

end
str=input("Look at the plot comparison of original density in green and its expansion in terms of normal density in red")
end
end

.
.
.
function [Median] = CalculateMedianOfGenGammaBisection(d,p,a)


Mu=a*gamma((d+1)/p)/gamma(d/p)


Med=0.576809914575588
MedCdf=(gamma(d/p)-igamma(d/p,(Med/a)^p))/gamma(d/p)
MuCdf=(gamma(d/p)-igamma(d/p,(Mu/a)^p))/gamma(d/p)
if (MuCdf<.5)
    BLow=Mu;
    BHigh=BLow*1.1;
    HCdf=(gamma(d/p)-igamma(d/p,(BHigh/a)^p))/gamma(d/p)
    while (HCdf<.5)
        BLow=BHigh;
        BHigh=BHigh*1.1;
        HCdf=(gamma(d/p)-igamma(d/p,(BHigh/a)^p))/gamma(d/p)
    end
    
end
if (MuCdf>.5)
    BHigh=Mu;
    BLow=BHigh/1.1;
    LCdf=(gamma(d/p)-igamma(d/p,(BLow/a)^p))/gamma(d/p);
    while (LCdf>.5)
        BHigh=BLow;
        BLow=BHigh/1.1;
        LCdf=(gamma(d/p)-igamma(d/p,(BLow/a)^p))/gamma(d/p);
    end
end

for nn=1:35
    BMid=(BLow+BHigh)/2;
    MidCdf=(gamma(d/p)-igamma(d/p,(BMid/a)^p))/gamma(d/p);
    if(MidCdf>.5)
        BHigh=BMid;
    end
    if(MidCdf<.5)
        BLow=BMid;
    end
end

Median=(BHigh+BLow)/2;

end

,
,
,
function [pg,dpdg] = CalculateGeneralizedGammaDensityDerivatives(gMedian,d,p,a)


g=gMedian;

g1=g^(d-1);
g2=exp(-(g/a)^p);


pg=g1*g2*(p/a^d)/gamma(d/p);

dpdg1(1)=(d-1)*g^(d-2);
for nn=2:9
    dpdg1(nn)=dpdg1(nn-1)*(d-nn)/g;
end
    

dpdg2(1)=-((exp(-(g/a)^p)* (g/a)^p *p)/g);
dpdg2(2)=(exp(-(g/a)^p)* (g/a)^p* p *(1+(-1+(g/a)^p)* p))/g^2;
dpdg2(3)=-((exp(-(g/a)^p)* (g/a)^p* p* (2+3 *(-1+(g/a)^p)* p+ ...
    (1-3* (g/a)^p+(g/a)^(2* p)) *p^2))/g^3);
dpdg2(4)=(1/(g^4))*exp(-(g/a)^p) * (g/a)^p *p *(6+11* (-1+(g/a)^p)* p+ ...
    6 *(1-3 *(g/a)^p+(g/a)^(2* p))* p^2+ ...
    (-1+7 *(g/a)^p- 6 *(g/a)^(2* p)+(g/a)^(3* p))* p^3);
dpdg2(5)=-(1/(g^5))*exp(-(g/a)^p) *(g/a)^p* p *(24+50* (-1+(g/a)^p)* p+  ...
    35 *(1-3* (g/a)^p+(g/a)^(2* p))* p^2+ ...
    10* (-1+7 *(g/a)^p-6* (g/a)^(2* p)+(g/a)^(3* p))* p^3+ ...
    (1-15 *(g/a)^p+25 *(g/a)^(2* p)-10 *(g/a)^(3* p)+(g/a)^(4* p))* p^4);
dpdg2(6)=(1/(g^6))*exp(-(g/a)^p) *(g/a)^p* p *(120+274* (-1+(g/a)^p)* p+ ...
    225 *(1-3 *(g/a)^p+(g/a)^(2 *p))* p^2+ ...
    85 *(-1+7 *(g/a)^p-6 *(g/a)^(2* p)+(g/a)^(3* p))* p^3+ ...
    15 *(1-15 *(g/a)^p+25 *(g/a)^(2* p)-10* (g/a)^(3* p)+(g/a)^(4 *p))* p^4+ ...
    (-1+31* (g/a)^p-90* (g/a)^(2 *p)+65 *(g/a)^(3* p)-15* (g/a)^(4* p)+(g/a)^(5 *p))* p^5);
dpdg2(7)=-(1/(g^7))*exp(-(g/a)^p)* (g/a)^p* p* (720+1764* (-1+(g/a)^p) *p+1624 *(1-3 *(g/a)^p+(g/a)^(2 *p)) *p^2+ ...
    735* (-1+7* (g/a)^p-6* (g/a)^(2* p)+(g/a)^(3* p))* p^3+175* (1-15 *(g/a)^p+25* (g/a)^(2 *p)- ...
    10* (g/a)^(3 *p)+(g/a)^(4 *p))* p^4+21 *(-1+31* (g/a)^p-90* (g/a)^(2* p)+65 *(g/a)^(3* p)- ...
    15* (g/a)^(4* p)+(g/a)^(5* p))* p^5+(1-63* (g/a)^p+301* (g/a)^(2* p)-350* (g/a)^(3 *p)+140 *(g/a)^(4* p)- ...
    21 *(g/a)^(5* p)+(g/a)^(6* p))* p^6);

dpdg2(8)=(1/(g^8))*exp(-(g/a)^p)* (g/a)^p* p* (5040+13068* (-1+(g/a)^p)* p+13132 *(1-3* (g/a)^p+(g/a)^(2* p))* p^2+ ...
    6769* (-1+7* (g/a)^p-6* (g/a)^(2 *p)+(g/a)^(3* p))* p^3+1960 *(1-15 *(g/a)^p+25* (g/a)^(2* p)-10 *(g/a)^(3* p)+(g/a)^(4* p))* p^4+322* (-1+31* (g/a)^p-90* (g/a)^(2* p)+ ...
    65* (g/a)^(3* p)-15 *(g/a)^(4* p)+(g/a)^(5* p))* p^5+28* (1-63* (g/a)^p+301* (g/a)^(2* p)-350* (g/a)^(3* p)+ ...
    140* (g/a)^(4* p)-21 *(g/a)^(5* p)+(g/a)^(6* p))* p^6+(-1+127* (g/a)^p-966* (g/a)^(2* p)+1701* (g/a)^(3 *p)- ...
    1050* (g/a)^(4* p)+266* (g/a)^(5* p)-28* (g/a)^(6* p)+(g/a)^(7* p))* p^7);

dpdg2(9)=-(1/(g^9)).*exp(-(g/a)^p)* (g/a)^p* p* (40320+109584* (-1+(g/a)^p)* p+118124 .*(1-3 *(g/a)^p+(g/a)^(2* p))* p^2+67284* (-1+7 *(g/a)^p-6* (g/a)^(2* p)+(g/a)^(3* p))* p^3+ ...
    22449 *(1-15* (g/a)^p+25* (g/a)^(2* p)-10* (g/a)^(3 *p)+(g/a)^(4* p))* p^4+4536* (-1+31 *(g/a)^p-90* (g/a)^(2 *p)+65* (g/a)^(3* p)-15* (g/a)^(4* p)+(g/a)^(5 *p))* p^5+ ...
    546* (1-63* (g/a)^p+301* (g/a)^(2* p)-350* (g/a)^(3* p)+140* (g/a)^(4* p)-21* (g/a)^(5* p)+(g/a)^(6* p))* p^6+ ...
    36* (-1+127 *(g/a)^p-966 *(g/a)^(2* p)+1701* (g/a)^(3* p)-1050* (g/a)^(4* p)+266* (g/a)^(5* p)-28 *(g/a)^(6* p)+(g/a)^(7 *p))* p^7+ ...
    (1-255* (g/a)^p+3025 *(g/a)^(2 *p)-7770 *(g/a)^(3 *p)+6951* (g/a)^(4* p)-2646 *(g/a)^(5* p)+462* (g/a)^(6 *p)-36 *(g/a)^(7* p)+(g/a)^(8* p))* p^8);



for nn=1:9
    dpdg(nn)=0;
    dpdg(nn)=dpdg(nn)+g1*dpdg2(nn)+dpdg1(nn).*g2;
    for mm=2:nn
        kk=mm-1;
        ll=nn-mm+1;
        dpdg(nn)=dpdg(nn)+factorial(nn)/factorial(nn-mm+1)/factorial(mm-1).*dpdg1(kk).*dpdg2(ll);
    end
end

dpdg=dpdg.*(p/a^d)/gamma(d/p);


end

.
.
.
function [XGamma] = InvertCDFOfGenGammaBisection(CDF,d,p,a)


Mu=a*gamma((d+1)/p)/gamma(d/p)


%Med=0.576809914575588
%MedCdf=(gamma(d/p)-igamma(d/p,(Med/a)^p))/gamma(d/p)
MuCdf=(gamma(d/p)-igamma(d/p,(Mu/a)^p))/gamma(d/p)
if (MuCdf<CDF)
    BLow=Mu;
    BHigh=BLow*1.1;
    HCdf=(gamma(d/p)-igamma(d/p,(BHigh/a)^p))/gamma(d/p)
    while (HCdf<CDF)
        BLow=BHigh;
        BHigh=BHigh*1.1;
        HCdf=(gamma(d/p)-igamma(d/p,(BHigh/a)^p))/gamma(d/p)
    end
    
end
if (MuCdf>CDF)
    BHigh=Mu;
    BLow=BHigh/1.1;
    LCdf=(gamma(d/p)-igamma(d/p,(BLow/a)^p))/gamma(d/p)
    while (LCdf>CDF)
        BHigh=BLow;
        BLow=BHigh/1.1;
        LCdf=(gamma(d/p)-igamma(d/p,(BLow/a)^p))/gamma(d/p)
    end
end

for nn=1:35
    BMid=(BLow+BHigh)/2;
    MidCdf=(gamma(d/p)-igamma(d/p,(BMid/a)^p))/gamma(d/p);
    if(MidCdf>CDF)
        BHigh=BMid;
    end
    if(MidCdf<CDF)
        BLow=BMid;
    end
end

XGamma=(BHigh+BLow)/2;

end

.
.
.
function [px,dpdx] = CalculateLognormalDensityDerivatives(Xmedian,mu,sigma)


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);      


  
  dpdx(1)=dpxdX;
  dpdx(2)=d2pxdX2;
  dpdx(3)=d3pxdX3;
  dpdx(4)=d4pxdX4;
  dpdx(5)=d5pxdX5;
  dpdx(6)=d6pxdX6;
  dpdx(7)=d7pxdX7;
  
  
end

,
,
,
function [XLognormal] = InvertCDFOfLognormalBisection(CDF,mu,sigma)


Mu=exp(mu);


%Med=0.576809914575588
%MedCdf=(gamma(d/p)-igamma(d/p,(Med/a)^p))/gamma(d/p)
MuCdf=normcdf((log(Mu)-mu)/sigma);
if (MuCdf<CDF)
    BLow=Mu;
    BHigh=BLow*1.1;
    HCdf=normcdf((log(BHigh)-mu)/sigma)
    while (HCdf<CDF)
        BLow=BHigh;
        BHigh=BHigh*1.1;
        HCdf=normcdf((log(BHigh)-mu)/sigma)
    end
    
end
if (MuCdf>CDF)
    BHigh=Mu;
    BLow=BHigh/1.1;
    LCdf=normcdf((log(BLow)-mu)/sigma)
    while (LCdf>CDF)
        BHigh=BLow;
        BLow=BHigh/1.1;
        LCdf=normcdf((log(BLow)-mu)/sigma)
    end
end

for nn=1:35
    BMid=(BLow+BHigh)/2;
    MidCdf=normcdf((log(BMid)-mu)/sigma);
    if(MidCdf>CDF)
        BHigh=BMid;
    end
    if(MidCdf<CDF)
        BLow=BMid;
    end
end

XLognormal=(BHigh+BLow)/2;

end
.
.
.
function [pz,dpdz] = CalculateStandardNormalDensityDerivatives(Zm)


Hz=HermitePoly(10);
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);

for nn=1:10
    dpdz(nn)=Hz(nn+1,1);
    for mm=1:nn
        dpdz(nn)=dpdz(nn)+Hz(nn+1,mm+1).*Zm^(mm);
    end
    dpdz(nn)=(-1)^nn.*dpdz(nn).*exp(-0.5* (Zm.^2))/sqrt(2*pi);
end



end

.
.
.
function [dX2dX1] = CalculateDerivsBetweenRVarsFromDensityDerivs02(px1,px2,dp1dx1,dp2dx2,OrderOfExpansion)



% px1=px1;
% dp1dx1(1)=dp1dx1(1);
% dp1dx1(2)=dp1dx1(2);
% dp1dx1(3)=dp1dx1(3);
% dp1dx1(4)=dp1dx1(4);
% dp1dx1(5)=dp1dx1(5);
% dp1dx1(6)=dp1dx1(6);
% dp1dx1(7)=dp1dx1(7);
% dp1dx1(8)=dp1dx1(8);
% dp1dx1(9)=dp1dx1(9);
% 
% px2=px2;
% dp2dx2(1)=dp2dx2(1);
% dp2dx2(2)=dp2dx2(2);
% dp2dx2(3)=dp2dx2(3);
% dp2dx2(4)=dp2dx2(4);
% dp2dx2(5)=dp2dx2(5);
% dp2dx2(6)=dp2dx2(6);
% dp2dx2(7)=dp2dx2(7);
% dp2dx2(8)=dp2dx2(8);
% dp2dx2(9)=dp2dx2(9);

if(OrderOfExpansion>=1)
dX2dX1(1)=px1/px2;
end
if(OrderOfExpansion>=2)
dX2dX1(2)=(dp1dx1(1)-(dp2dx2(1) *dX2dX1(1)^2))/px2; 
end

if(OrderOfExpansion>=3)
dX2dX1(3)=(dp1dx1(2)- (2* dp2dx2(1)* dX2dX1(1) *dX2dX1(2)+dX2dX1(1) *(dX2dX1(1)^2* dp2dx2(2)+dp2dx2(1) *dX2dX1(2))))/px2 ;
end

if(OrderOfExpansion>=4)
dX2dX1(4)=(dp1dx1(3)-(3* dX2dX1(2)* (dX2dX1(1)^2* dp2dx2(2)+dp2dx2(1) *dX2dX1(2))+3* dp2dx2(1) *dX2dX1(1) *dX2dX1(3)+dX2dX1(1)* (3 *dX2dX1(1)* dp2dx2(2)* dX2dX1(2)+dX2dX1(1)^3* dp2dx2(3)+dp2dx2(1)* dX2dX1(3))))/px2 ;
end

if(OrderOfExpansion>=5)
dX2dX1(5)=(dp1dx1(4)- (6 *(dX2dX1(1)^2 *dp2dx2(2)+dp2dx2(1) *dX2dX1(2)) *dX2dX1(3)+4 *dX2dX1(2)* (3* dX2dX1(1) *dp2dx2(2) *dX2dX1(2)+dX2dX1(1)^3* dp2dx2(3)+dp2dx2(1)* dX2dX1(3))+4* dp2dx2(1)* dX2dX1(1)* dX2dX1(4)+dX2dX1(1) *(3* dp2dx2(2)* dX2dX1(2)^2+6* dX2dX1(1)^2* dX2dX1(2)* dp2dx2(3)+4 *dX2dX1(1)* dp2dx2(2)* dX2dX1(3)+dX2dX1(1)^4* dp2dx2(4)+dp2dx2(1)* dX2dX1(4))))/px2 ;
end

if(OrderOfExpansion>=6)
dX2dX1(6)=(dp1dx1(5)- (10* dX2dX1(3) *(3 *dX2dX1(1)* dp2dx2(2)* dX2dX1(2)+dX2dX1(1)^3* dp2dx2(3)+dp2dx2(1) *dX2dX1(3))+10* (dX2dX1(1)^2* dp2dx2(2)+dp2dx2(1) *dX2dX1(2))* dX2dX1(4)+5* dX2dX1(2)* (3* dp2dx2(2)* dX2dX1(2)^2+6* dX2dX1(1)^2* dX2dX1(2)* dp2dx2(3)+4 *dX2dX1(1)* dp2dx2(2)* dX2dX1(3)+dX2dX1(1)^4* dp2dx2(4)+dp2dx2(1)* dX2dX1(4))+5 *dp2dx2(1)* dX2dX1(1)* dX2dX1(5)+dX2dX1(1) *(15* dX2dX1(1)* dX2dX1(2)^2* dp2dx2(3)+10* dp2dx2(2)* dX2dX1(2)* dX2dX1(3)+10* dX2dX1(1)^2 *dp2dx2(3)* dX2dX1(3)+10* dX2dX1(1)^3* dX2dX1(2)* dp2dx2(4)+5* dX2dX1(1)* dp2dx2(2)* dX2dX1(4)+dX2dX1(1)^5* dp2dx2(5)+dp2dx2(1)* dX2dX1(5))))/px2 ;
end


if(OrderOfExpansion>=7)
dX2dX1(7)=(dp1dx1(6)- (20 *(3* dX2dX1(1)* dp2dx2(2)* dX2dX1(2)+dX2dX1(1)^3 *dp2dx2(3)+dp2dx2(1) *dX2dX1(3))* dX2dX1(4)+15* dX2dX1(3)* (3 *dp2dx2(2) *dX2dX1(2)^2+6* dX2dX1(1)^2* dX2dX1(2)* dp2dx2(3)+4 *dX2dX1(1)* dp2dx2(2)* dX2dX1(3)+dX2dX1(1)^4* dp2dx2(4)+dp2dx2(1)* dX2dX1(4))+15* (dX2dX1(1)^2* dp2dx2(2)+dp2dx2(1)* dX2dX1(2))* dX2dX1(5)+6* dX2dX1(2)* (15* dX2dX1(1)* dX2dX1(2)^2* dp2dx2(3)+10* dp2dx2(2)* dX2dX1(2) *dX2dX1(3)+10 *dX2dX1(1)^2 *dp2dx2(3) *dX2dX1(3)+10 *dX2dX1(1)^3 *dX2dX1(2)* dp2dx2(4)+5 *dX2dX1(1) *dp2dx2(2) *dX2dX1(4)+dX2dX1(1)^5 *dp2dx2(5)+dp2dx2(1)* dX2dX1(5))+6* dp2dx2(1)* dX2dX1(1) *dX2dX1(6)+dX2dX1(1)* (15 *dX2dX1(2)^3 *dp2dx2(3)+60 *dX2dX1(1) *dX2dX1(2) *dp2dx2(3)* dX2dX1(3)+10 *dp2dx2(2) *dX2dX1(3)^2+45 *dX2dX1(1)^2* dX2dX1(2)^2 *dp2dx2(4)+20* dX2dX1(1)^3 *dX2dX1(3)* dp2dx2(4)+15 *dp2dx2(2)* dX2dX1(2) *dX2dX1(4)+15 *dX2dX1(1)^2* dp2dx2(3) *dX2dX1(4)+15 *dX2dX1(1)^4* dX2dX1(2)* dp2dx2(5)+6 *dX2dX1(1) *dp2dx2(2)* dX2dX1(5)+dX2dX1(1)^6* dp2dx2(6)+dp2dx2(1) *dX2dX1(6))))/px2 ;
end

if(OrderOfExpansion>=8)
dX2dX1(8)=(dp1dx1(7)- (35 *dX2dX1(4) *(3 *dp2dx2(2)* dX2dX1(2)^2+6* dX2dX1(1)^2 *dX2dX1(2) *dp2dx2(3)+4* dX2dX1(1)* dp2dx2(2)* dX2dX1(3)+dX2dX1(1)^4 *dp2dx2(4)+dp2dx2(1) *dX2dX1(4))+35 *(3* dX2dX1(1) *dp2dx2(2) *dX2dX1(2)+dX2dX1(1)^3 *dp2dx2(3)+dp2dx2(1)* dX2dX1(3))* dX2dX1(5)+21* dX2dX1(3) *(15 *dX2dX1(1)* dX2dX1(2)^2 *dp2dx2(3)+10 *dp2dx2(2) *dX2dX1(2) *dX2dX1(3)+10* dX2dX1(1)^2 *dp2dx2(3) *dX2dX1(3)+10 *dX2dX1(1)^3 *dX2dX1(2) *dp2dx2(4)+5* dX2dX1(1) *dp2dx2(2)* dX2dX1(4)+dX2dX1(1)^5 *dp2dx2(5)+dp2dx2(1) *dX2dX1(5))+21 *(dX2dX1(1)^2 *dp2dx2(2)+dp2dx2(1)* dX2dX1(2)) *dX2dX1(6)+7 *dX2dX1(2)* (15 *dX2dX1(2)^3 *dp2dx2(3)+60 *dX2dX1(1) *dX2dX1(2) *dp2dx2(3) *dX2dX1(3)+10 *dp2dx2(2)* dX2dX1(3)^2+45* dX2dX1(1)^2 *dX2dX1(2)^2 *dp2dx2(4)+20* dX2dX1(1)^3* dX2dX1(3) *dp2dx2(4)+15 *dp2dx2(2) *dX2dX1(2) *dX2dX1(4)+15 *dX2dX1(1)^2 *dp2dx2(3) *dX2dX1(4)+15* dX2dX1(1)^4 *dX2dX1(2) *dp2dx2(5)+6 *dX2dX1(1)* dp2dx2(2) *dX2dX1(5)+dX2dX1(1)^6 *dp2dx2(6)+dp2dx2(1) *dX2dX1(6))+7 *dp2dx2(1)* dX2dX1(1)* dX2dX1(7)+dX2dX1(1) *(105 *dX2dX1(2)^2 *dp2dx2(3) *dX2dX1(3)+70 *dX2dX1(1)* dp2dx2(3) *dX2dX1(3)^2+105 *dX2dX1(1) *dX2dX1(2)^3 *dp2dx2(4)+210* dX2dX1(1)^2 *dX2dX1(2) *dX2dX1(3) *dp2dx2(4)+105* dX2dX1(1) *dX2dX1(2) *dp2dx2(3) *dX2dX1(4)+35 *dp2dx2(2)* dX2dX1(3) *dX2dX1(4)+35 *dX2dX1(1)^3* dp2dx2(4) *dX2dX1(4)+105* dX2dX1(1)^3 *dX2dX1(2)^2* dp2dx2(5)+35* dX2dX1(1)^4* dX2dX1(3) *dp2dx2(5)+21 *dp2dx2(2) *dX2dX1(2) *dX2dX1(5)+21 *dX2dX1(1)^2 *dp2dx2(3) *dX2dX1(5)+21 *dX2dX1(1)^5 *dX2dX1(2) *dp2dx2(6)+7 *dX2dX1(1)* dp2dx2(2)* dX2dX1(6)+dX2dX1(1)^7* dp2dx2(7)+dp2dx2(1) *dX2dX1(7))))/px2 ;
end

if(OrderOfExpansion>=9)
dX2dX1(9)=(dp1dx1(8)- (70* (3 *dp2dx2(2)* dX2dX1(2)^2+6* dX2dX1(1)^2* dX2dX1(2)* dp2dx2(3)+4* dX2dX1(1)* dp2dx2(2)* dX2dX1(3)+dX2dX1(1)^4 *dp2dx2(4)+dp2dx2(1)* dX2dX1(4)) *dX2dX1(5)+56* dX2dX1(4)* (15* dX2dX1(1)* dX2dX1(2)^2 *dp2dx2(3)+10 *dp2dx2(2) *dX2dX1(2) *dX2dX1(3)+10 *dX2dX1(1)^2 *dp2dx2(3) *dX2dX1(3)+10 *dX2dX1(1)^3 *dX2dX1(2) *dp2dx2(4)+5 *dX2dX1(1) *dp2dx2(2)* dX2dX1(4)+dX2dX1(1)^5 *dp2dx2(5)+dp2dx2(1) *dX2dX1(5))+56 *(3 *dX2dX1(1) *dp2dx2(2) *dX2dX1(2)+dX2dX1(1)^3* dp2dx2(3)+dp2dx2(1) *dX2dX1(3)) *dX2dX1(6)+28* dX2dX1(3) *(15 *dX2dX1(2)^3 *dp2dx2(3)+60 *dX2dX1(1) *dX2dX1(2) *dp2dx2(3) *dX2dX1(3)+10 *dp2dx2(2)* dX2dX1(3)^2+45 *dX2dX1(1)^2 *dX2dX1(2)^2 *dp2dx2(4)+20* dX2dX1(1)^3 *dX2dX1(3) *dp2dx2(4)+15* dp2dx2(2) *dX2dX1(2) *dX2dX1(4)+15 *dX2dX1(1)^2 *dp2dx2(3) *dX2dX1(4)+15 *dX2dX1(1)^4 *dX2dX1(2) *dp2dx2(5)+6 *dX2dX1(1) *dp2dx2(2)* dX2dX1(5)+dX2dX1(1)^6 *dp2dx2(6)+dp2dx2(1) *dX2dX1(6))+28* (dX2dX1(1)^2 *dp2dx2(2)+dp2dx2(1) *dX2dX1(2)) *dX2dX1(7)+8 *dX2dX1(2) *(105 *dX2dX1(2)^2 *dp2dx2(3) *dX2dX1(3)+70* dX2dX1(1)* dp2dx2(3) *dX2dX1(3)^2+105 *dX2dX1(1)* dX2dX1(2)^3 *dp2dx2(4)+210 *dX2dX1(1)^2 *dX2dX1(2)* dX2dX1(3)* dp2dx2(4)+105 *dX2dX1(1)* dX2dX1(2) *dp2dx2(3) *dX2dX1(4)+35 *dp2dx2(2) *dX2dX1(3) *dX2dX1(4)+35* dX2dX1(1)^3 *dp2dx2(4) *dX2dX1(4)+105* dX2dX1(1)^3 *dX2dX1(2)^2 *dp2dx2(5)+35 *dX2dX1(1)^4 *dX2dX1(3) *dp2dx2(5)+21* dp2dx2(2)* dX2dX1(2) *dX2dX1(5)+21 *dX2dX1(1)^2 *dp2dx2(3)* dX2dX1(5)+21* dX2dX1(1)^5* dX2dX1(2) *dp2dx2(6)+7 *dX2dX1(1)* dp2dx2(2) *dX2dX1(6)+dX2dX1(1)^7 *dp2dx2(7)+dp2dx2(1) *dX2dX1(7))+8 *dp2dx2(1)* dX2dX1(1)* dX2dX1(8)+dX2dX1(1) *(280* dX2dX1(2) *dp2dx2(3) *dX2dX1(3)^2+105 *dX2dX1(2)^4 *dp2dx2(4)+840 *dX2dX1(1)* dX2dX1(2)^2 *dX2dX1(3) *dp2dx2(4)+280 *dX2dX1(1)^2 *dX2dX1(3)^2 *dp2dx2(4)+210 *dX2dX1(2)^2 *dp2dx2(3) *dX2dX1(4)+280 *dX2dX1(1)* dp2dx2(3) *dX2dX1(3) *dX2dX1(4)+420 *dX2dX1(1)^2 *dX2dX1(2) *dp2dx2(4) *dX2dX1(4)+35* dp2dx2(2) *dX2dX1(4)^2+420 *dX2dX1(1)^2 *dX2dX1(2)^3* dp2dx2(5)+560* dX2dX1(1)^3 *dX2dX1(2)* dX2dX1(3)* dp2dx2(5)+70* dX2dX1(1)^4 *dX2dX1(4) *dp2dx2(5)+168 *dX2dX1(1) *dX2dX1(2) *dp2dx2(3) *dX2dX1(5)+56 *dp2dx2(2) *dX2dX1(3) *dX2dX1(5)+56* dX2dX1(1)^3 *dp2dx2(4) *dX2dX1(5)+210 *dX2dX1(1)^4 *dX2dX1(2)^2 *dp2dx2(6)+56 *dX2dX1(1)^5 *dX2dX1(3) *dp2dx2(6)+28 *dp2dx2(2) *dX2dX1(2) *dX2dX1(6)+28* dX2dX1(1)^2* dp2dx2(3) *dX2dX1(6)+28 *dX2dX1(1)^6 *dX2dX1(2) *dp2dx2(7)+8 *dX2dX1(1) *dp2dx2(2) *dX2dX1(7)+dX2dX1(1)^8 *dp2dx2(8)+dp2dx2(1) *dX2dX1(8))))/px2 ;
end

if(OrderOfExpansion>=10)
dX2dX1(10)=(dp1dx1(9)- (126*dX2dX1(5)*(15*dX2dX1(1)*dX2dX1(2)^2*dp2dx2(3)+10*dp2dx2(2)*dX2dX1(2)*dX2dX1(3)+10*dX2dX1(1)^2*dp2dx2(3)*dX2dX1(3)+10*dX2dX1(1)^3*dX2dX1(2)*dp2dx2(4)+5*dX2dX1(1)*dp2dx2(2)*dX2dX1(4)+dX2dX1(1)^5*dp2dx2(5)+dp2dx2(1)*dX2dX1(5))+126*(3*dp2dx2(2)*dX2dX1(2)^2+6*dX2dX1(1)^2*dX2dX1(2)*dp2dx2(3)+4*dX2dX1(1)*dp2dx2(2)*dX2dX1(3)+dX2dX1(1)^4*dp2dx2(4)+dp2dx2(1)*dX2dX1(4))*dX2dX1(6)+84*dX2dX1(4)*(15*dX2dX1(2)^3*dp2dx2(3)+60*dX2dX1(1)*dX2dX1(2)*dp2dx2(3)*dX2dX1(3)+10*dp2dx2(2)*dX2dX1(3)^2+45*dX2dX1(1)^2*dX2dX1(2)^2*dp2dx2(4)+20*dX2dX1(1)^3*dX2dX1(3)*dp2dx2(4)+15*dp2dx2(2)*dX2dX1(2)*dX2dX1(4)+15*dX2dX1(1)^2*dp2dx2(3)*dX2dX1(4)+15*dX2dX1(1)^4*dX2dX1(2)*dp2dx2(5)+6*dX2dX1(1)*dp2dx2(2)*dX2dX1(5)+dX2dX1(1)^6*dp2dx2(6)+dp2dx2(1)*dX2dX1(6))+84*(3*dX2dX1(1)*dp2dx2(2)*dX2dX1(2)+dX2dX1(1)^3*dp2dx2(3)+dp2dx2(1)*dX2dX1(3))*dX2dX1(7)+36*dX2dX1(3)*(105*dX2dX1(2)^2*dp2dx2(3)*dX2dX1(3)+70*dX2dX1(1)*dp2dx2(3)*dX2dX1(3)^2+105*dX2dX1(1)*dX2dX1(2)^3*dp2dx2(4)+210*dX2dX1(1)^2*dX2dX1(2)*dX2dX1(3)*dp2dx2(4)+105*dX2dX1(1)*dX2dX1(2)*dp2dx2(3)*dX2dX1(4)+35*dp2dx2(2)*dX2dX1(3)*dX2dX1(4)+35*dX2dX1(1)^3*dp2dx2(4)*dX2dX1(4)+105*dX2dX1(1)^3*dX2dX1(2)^2*dp2dx2(5)+35*dX2dX1(1)^4*dX2dX1(3)*dp2dx2(5)+21*dp2dx2(2)*dX2dX1(2)*dX2dX1(5)+21*dX2dX1(1)^2*dp2dx2(3)*dX2dX1(5)+21*dX2dX1(1)^5*dX2dX1(2)*dp2dx2(6)+7*dX2dX1(1)*dp2dx2(2)*dX2dX1(6)+dX2dX1(1)^7*dp2dx2(7)+dp2dx2(1)*dX2dX1(7))+36*(dX2dX1(1)^2*dp2dx2(2)+dp2dx2(1)*dX2dX1(2))*dX2dX1(8)+9*dX2dX1(2)*(280*dX2dX1(2)*dp2dx2(3)*dX2dX1(3)^2+105*dX2dX1(2)^4*dp2dx2(4)+840*dX2dX1(1)*dX2dX1(2)^2*dX2dX1(3)*dp2dx2(4)+280*dX2dX1(1)^2*dX2dX1(3)^2*dp2dx2(4)+210*dX2dX1(2)^2*dp2dx2(3)*dX2dX1(4)+280*dX2dX1(1)*dp2dx2(3)*dX2dX1(3)*dX2dX1(4)+420*dX2dX1(1)^2*dX2dX1(2)*dp2dx2(4)*dX2dX1(4)+35*dp2dx2(2)*dX2dX1(4)^2+420*dX2dX1(1)^2*dX2dX1(2)^3*dp2dx2(5)+560*dX2dX1(1)^3*dX2dX1(2)*dX2dX1(3)*dp2dx2(5)+70*dX2dX1(1)^4*dX2dX1(4)*dp2dx2(5)+168*dX2dX1(1)*dX2dX1(2)*dp2dx2(3)*dX2dX1(5)+56*dp2dx2(2)*dX2dX1(3)*dX2dX1(5)+56*dX2dX1(1)^3*dp2dx2(4)*dX2dX1(5)+210*dX2dX1(1)^4*dX2dX1(2)^2*dp2dx2(6)+56*dX2dX1(1)^5*dX2dX1(3)*dp2dx2(6)+28*dp2dx2(2)*dX2dX1(2)*dX2dX1(6)+28*dX2dX1(1)^2*dp2dx2(3)*dX2dX1(6)+28*dX2dX1(1)^6*dX2dX1(2)*dp2dx2(7)+8*dX2dX1(1)*dp2dx2(2)*dX2dX1(7)+dX2dX1(1)^8*dp2dx2(8)+dp2dx2(1)*dX2dX1(8))+9*dp2dx2(1)*dX2dX1(1)*dX2dX1(9)+dX2dX1(1)*(280*dp2dx2(3)*dX2dX1(3)^3+1260*dX2dX1(2)^3*dX2dX1(3)*dp2dx2(4)+2520*dX2dX1(1)*dX2dX1(2)*dX2dX1(3)^2*dp2dx2(4)+1260*dX2dX1(2)*dp2dx2(3)*dX2dX1(3)*dX2dX1(4)+1890*dX2dX1(1)*dX2dX1(2)^2*dp2dx2(4)*dX2dX1(4)+1260*dX2dX1(1)^2*dX2dX1(3)*dp2dx2(4)*dX2dX1(4)+315*dX2dX1(1)*dp2dx2(3)*dX2dX1(4)^2+945*dX2dX1(1)*dX2dX1(2)^4*dp2dx2(5)+3780*dX2dX1(1)^2*dX2dX1(2)^2*dX2dX1(3)*dp2dx2(5)+840*dX2dX1(1)^3*dX2dX1(3)^2*dp2dx2(5)+1260*dX2dX1(1)^3*dX2dX1(2)*dX2dX1(4)*dp2dx2(5)+378*dX2dX1(2)^2*dp2dx2(3)*dX2dX1(5)+504*dX2dX1(1)*dp2dx2(3)*dX2dX1(3)*dX2dX1(5)+756*dX2dX1(1)^2*dX2dX1(2)*dp2dx2(4)*dX2dX1(5)+126*dp2dx2(2)*dX2dX1(4)*dX2dX1(5)+126*dX2dX1(1)^4*dp2dx2(5)*dX2dX1(5)+1260*dX2dX1(1)^3*dX2dX1(2)^3*dp2dx2(6)+1260*dX2dX1(1)^4*dX2dX1(2)*dX2dX1(3)*dp2dx2(6)+126*dX2dX1(1)^5*dX2dX1(4)*dp2dx2(6)+252*dX2dX1(1)*dX2dX1(2)*dp2dx2(3)*dX2dX1(6)+84*dp2dx2(2)*dX2dX1(3)*dX2dX1(6)+84*dX2dX1(1)^3*dp2dx2(4)*dX2dX1(6)+378*dX2dX1(1)^5*dX2dX1(2)^2*dp2dx2(7)+84*dX2dX1(1)^6*dX2dX1(3)*dp2dx2(7)+36*dp2dx2(2)*dX2dX1(2)*dX2dX1(7)+36*dX2dX1(1)^2*dp2dx2(3)*dX2dX1(7)+36*dX2dX1(1)^7*dX2dX1(2)*dp2dx2(8)+9*dX2dX1(1)*dp2dx2(2)*dX2dX1(8)+dX2dX1(1)^9*dp2dx2(9)+dp2dx2(1)*dX2dX1(9))))/px2;
end


% dX2dX1(1)=dX2dX1(1);
% dX2dX1(2)=dX2dX1(2);
% dX2dX1(3)=dX2dX1(3);
% dX2dX1(4)=dX2dX1(4);
% dX2dX1(5)=dX2dX1(5);
% dX2dX1(6)=dX2dX1(6);
% dX2dX1(7)=dX2dX1(7);
% dX2dX1(8)=dX2dX1(8);
% dX2dX1(9)=dX2dX1(9);
% dX2dX1(10)=dX2dX1(10);





end

.
.
.
function [Hz] = HermitePoly(Order)
Hz(1:Order+1,1:Order+1)=0.0;
Hz(1,1)=1.0;
for k=1:Order
    for m=1:k
       Hz(k+1,m+1)=Hz(k,m);
       if (m>1)
            Hz(k+1,m-1)=Hz(k+1,m-1)-(m-1)*Hz(k,m);
       end
    end
end

end

.
.
.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 8th, 2024, 9:31 am

Here is what you see in graph when you run the previous program. Later today, I will try to expand Slater type orbitals in terms of a common Gamma density and also in terms of standard normal density.
Please make sure you rescale the matlab graph from edit menu and look at the body of the graph since many times tails of Generalized gamma might be extremely long to millions while most body of the graph might only be only a few tens. Green is true density so please make sure you rescale the graph so that you can focus on the body of green density.

Image
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 2586
Joined: July 14th, 2002, 3:00 am

Re: Breakthrough in the theory of stochastic differential equations and their simulation

April 10th, 2024, 7:13 pm

Friends, American defense and their backers are thoroughly adamant to continue my persecution at all cost and they try to undo anything that gives me relief from mind control or allows me to work creatively. We had Eid-ul-fitr today which is celebrated by everyone in our society so today (and yesterday) there was some lull in the storm since most agents would be celebrating Eid now. But till two days ago, they were trying their best to drug my food and water all over the city. And a few days ago, I was hit again and again by bad food. 
They continue to use all sorts of gases and charges to remove my neurotransmitters. Gases have decreased somewhat than earlier times but they still use gas on me several times a day. However they did not use charges on me regularly(gas filled with charges to charge and static on my body). They would use a large number of gases earlier but they almost never used gas charges on me(to the extent that I could feel them like needles light or sharp). They somehow use potential energy and the charges from the air flow all the way from bottom of the pants where they get inside the pants and then flow upwards inside the pants towards the testicles bag and settle on that. When they hit the body they feel like light needles(depending upon the intensity of charges and highly charged particles can feel literally like sharp needles.) These charges are very disturbing and thoroughly discomforting and it becomes impossible to concentrate on anything when they release these charges in the air. Previously mind control agents used to be very high on gases and they were low on charges but now they are high on charges and relatively low on torture gases that cause sickening feeling in lungs and the stomach.
I washed my clothes last night and I believe that they released charges in tap water and when I wore some of the cloths that I had washed yesterday that I realized that my clothes were easily getting charged and I could feel needles in  my body after wearing those clothes today. I realized that they had added mind control chemicals in tap water since they had anticipated that I would wash clothes yesterday.
During my sleep, they still continue to force mind control chemicals inside my mouth that settles inside the mouth. 
When I go to the city on motorbike to get food, they continue to force waves into my back using devices hidden inside the seat and sometimes it   would become very very painful to sit on the motorbike for hours. I recently started using a plastic shopping bag and filled it with several aluminium foils from snacks and also with some cloth and I would place it over the seat and then drive the motorbike but it only helps very partially.
Finding good water in the city is very difficult. I told friends that a water brand from Murree Brewery was good water but Pak army approached Murree Brewery and the water supplied in recent dates is again mixed with mind control chemicals. And Nestle and Dasani water etc. have not turned better at all. If I get good ground water from some portion of the city, it is made sure by Pak army that ground water in that part of the city is thoroughly drugged and whenever I retry to get good water after a day or two I cannot get good water from any place in that particular area.
Lundbeck has supplied antipsychotic Fluanxol injections all over the city that carry mind control drugs and I was lucky this time to be able to get a better injection but I am very afraid that getting a good injection would be extremely difficult next month. 
 
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal