Breakthrough in the theory of stochastic differential equations and their simulation

Amin · November 10th, 2023, 8:13 am

Friends, please notice the part of main program where I have specified MaxCount and MinCount in the code below. I noticed that this is relevant to half a million paths and have to be proportionally altered when total number of paths change. MaxCount caps the maximum paths that can be allocated to any volatility grid cell. In my implementation, Total number of paths, MaxCount and MinCount have to be multiple of each other (as I have used) but the grid can work for non-multiples as well(I have not tested those cases).
When I tried the program out to four years at the start after writing it, it was not working for lognormal type asset very well but when I decreased the MaxCount or alternatively increased sampling frequency, lognormal type asset was working far better.
This ad-hoc increase in sampling frequency might still not be perfectly optimal. But I learnt that sampling has to be proportional to local volatility of the asset in the subdivision. I still have not implemented this local volatility in design of the volatility grid cells. This is still a hypothesis but this hypothesis has enough weight.
I believe that sampling of the data by the polynomial grid has to be more carefully studied but I think that local polynomial sampling frequency has to be proportional to both probability of the SV and variance/volatility of the asset in the relevant grid cell. I must warn that this is just a hypothesis at this stage and needs to be more carefully checked.
My grid does not implement local asset volatility and is refined universally for asset volatility in right tail as I have used greater sampling frequency for lognormal type asset and lesser sampling frequency for asset volatility exponent equal to .5
I will try some better design for underlying grid to see if the hypothesis does hold and come back to friends. But it does make sense to think that sampling frequency should be proportional to how fast the underlying variable being modelled is changing locally.
There is another possible improvement point. In the current program, sampling frequency(which is not equi-probable and is not on a probability grid but is on a Z-grid) for calculation of call prices is totally independent of the sampling frequency elsewhere and only related to rest of the program through Z-series coefficients of conditional Coefficients of individual Z-series across the volatility grid. I also want to see how altering that sampling frequency of Z-grid for calculation of call prices increases or decreases the accuracy of call prices.

MaxCount=2000;
%For asset SDE volatility exponents close to lognormal, please choose
%MaxCount=500;
%For asset SDE volatility exponents close to .75, please choose
%MaxCount=1000;
%For asset SDE volatility exponents close to .5, please choose
%MaxCount=2000;

if(gammaX>.87)
    MaxCount=500;
elseif((gammaX<=.87)&&(gammaX>=.67))
    MaxCount=1000;
elseif(gammaX<.67)
    MaxCount=2000;
end
MinCount=100;

.
.

Amin · November 10th, 2023, 9:18 am

Friends, I want to share the details of my small algorithm to correct the moments of arrays of standard normals associated with equiprobable normal grid subdivisions. These arrays of standard normals would later be used for polynomial regression on the data to find parameters of Z-series from the data.

Suppose there are N data members and we want to do a hermite polynomial regression on this data to find Z-series parameters of data random variable. We will divide the standard normal density into N subdivisions with each subdivision sharing an equal probability mass of 1/N.
However discrete values of Zn when integrated on probability mass 1/N, can have moments different from moments of standard normal density which will create inaccuracies in our polynomial regression to find Z-series parameters of the data random variable.

I thought of multiplying the data array with an expression comprising sum of hermite polynomials with unspecified coefficients. We would write equations of moments and then find coefficients of hermite polynomials that pin down the moments of new data precisely equal to target moments.

One consideration is that our Zn array is perfectly symmetric around zero so all odd moments are zero. We want to add only even hermite polynomials in the expression of hermite polynomials with unspecified coefficients(that is multiplied to original Zn array to make its moments equal to target moments). We do not want odd hermite polynomials in this expression with unspecified coefficients as they would introduce skew in our array of equi-probable standard normal values(and then we would have to do unnecessary work to cancel the skew). In order to have zero odd moments(as in case of standard normals), we only multiply with a series in even hermite polynomials with unspecified coefficients that preserves perfect symmetry around zero so we would not have to bother with any odd moments and all odd moments will be zero by construction due to symmetry of our Zn array around zero.

Our original standard normal Zn array associated with equiprobable grid is multiplied with a series in hermite polynomials as

[$]Z_{n,Corrected} \, = \, Z_n \, (\alpha \, + \beta \, (Z_n^2-1) \, + \, \gamma \, ({Z_n}^4 \, + \, 6 \, {Z_n}^2 \, + \, 3 \,))\,[$]

We find coefficients [$]\alpha[$], [$]\beta[$], and [$]\gamma[$] so that our objective function matching moments of [$]Z_{n,Corrected} \, [$] with target moments is maximized.

[$]Moment2\, = \sum\limits_{n=1}^{n=N} \, {Z_{n,Corrected}}^2 \, \frac{1}{N} [$]
[$]= \, \sum\limits_{n=1}^{n=N} \, {Z_n}^2 \, {(\alpha \, + \beta \, (Z_n^2-1) \, + \, \gamma \, ({Z_n}^4 \, + \, 6 \, {Z_n}^2 \, + \, 3 \,))}^2\, \frac{1}{N} [$]
[$]Moment4\, = \sum\limits_{n=1}^{n=N} \, {Z_{n,Corrected}}^4 \, \frac{1}{N} [$]
[$]= \, \sum\limits_{n=1}^{n=N} \, {Z_n}^4 \, {(\alpha \, + \beta \, (Z_n^2-1) \, + \, \gamma \, ({Z_n}^4 \, + \, 6 \, {Z_n}^2 \, + \, 3 \,))}^4\, \frac{1}{N} [$]
[$]Moment6\, = \sum\limits_{n=1}^{n=N} \, {Z_{n,Corrected}}^6 \, \frac{1}{N} [$]
[$]= \, \sum\limits_{n=1}^{n=N} \, {Z_n}^6 \, {(\alpha \, + \beta \, (Z_n^2-1) \, + \, \gamma \, ({Z_n}^4 \, + \, 6 \, {Z_n}^2 \, + \, 3 \,))}^6\, \frac{1}{N} [$]

Here is an iterative optimization matlab program that finds [$]Z_{n,Corrected}[$] whose moments match the standard normal moments with great precision. I improved this program and now it takes only 3-7 seconds as opposed to 30-40 secs it was taking yesterday. So if great speed is not the concern, you can easily embed it in your programs.

.
.
.
Code: Select all
function [ZnNew] = MatchMomentsOfIdealZ(Zn)

Ndata=length(Zn);
if(Ndata<=100)
    Iterations=3500;
end
if((Ndata>100)&&(Ndata<=500))
    Iterations=2000;
end
if((Ndata>500)&&(Ndata<=1000))
    Iterations=1500;
end
if((Ndata>1000)&&(Ndata<=3000))
    Iterations=1000;
end
if((Ndata>3000)&&(Ndata<=6000))
    Iterations=750;
end
if((Ndata>6000)&&(Ndata<=10000))
    Iterations=550;
end
if(Ndata>10000)
    Iterations=50;
end
alpha=1;
beta=0;
gamma=0;

ZnNew(1:Ndata)=Zn(1:Ndata).*(alpha+beta*(Zn(1:Ndata).^2-1)+gamma*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));

Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata

w1=5000;
w2=200;

ObjFuncBest=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;


delta_alpha=.00025;
delta_beta=.0002;
delta_gamma=.0002;
alphaSuccessCount=0;
betaSuccessCount=0;
gammaSuccessCount=0;
SuccessPrevalpha=0;
SuccessPrevbeta=0;
SuccessPrevgamma=0;

mul1=.9;
mul2=1.1;

for iter=1:Iterations
    alphaNew=alpha+delta_alpha;
    betaNew=beta;
    gammaNew=gamma;
    ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
    Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
    Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
    Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;

    ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;

    
    if(ObjFuncNew<ObjFuncBest)
        alpha=alphaNew;
        ObjFuncBest=ObjFuncNew;
        if(SuccessPrevalpha==1)
            alphaSuccessCount=alphaSuccessCount+1;
        end
        SuccessPrevalpha=1;
    else
        alphaNew=alpha-delta_alpha;
        betaNew=beta;
        gammaNew=gamma;
        ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
        Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
        Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
        Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;

        ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;
        if(ObjFuncNew<ObjFuncBest)
            alpha=alphaNew;
            ObjFuncBest=ObjFuncNew;
            if(SuccessPrevalpha==1)
                alphaSuccessCount=alphaSuccessCount+1;
            end
            SuccessPrevalpha=1;
        else
            SuccessPrevalpha=0;
            alphaSuccessCount=0;
            delta_alpha=delta_alpha*mul1;
        end
    end
    
    
    alphaNew=alpha;
    betaNew=beta+delta_beta;
    gammaNew=gamma;
    ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
    Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
    Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
    Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;

    ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;

    if(ObjFuncNew<ObjFuncBest)
        beta=betaNew;
        ObjFuncBest=ObjFuncNew;
        if(SuccessPrevbeta==1)
            betaSuccessCount=betaSuccessCount+1;
        end
        SuccessPrevbeta=1;
    else
        alphaNew=alpha;
        betaNew=beta-delta_beta;
        gammaNew=gamma;
        ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
        Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
        Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
        Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;

        ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;
        if(ObjFuncNew<ObjFuncBest)
        beta=betaNew;
        ObjFuncBest=ObjFuncNew;
        if(SuccessPrevbeta==1)
            betaSuccessCount=betaSuccessCount+1;
        end
        SuccessPrevbeta=1;
        else
            SuccessPrevbeta=0;
            betaSuccessCount=0;
            delta_beta=delta_beta*mul1;
        end
    end
    
    alphaNew=alpha;
    betaNew=beta;
    gammaNew=gamma+delta_gamma;
    ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
    Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
    Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
    Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;

    ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;

    if(ObjFuncNew<ObjFuncBest)
        gamma=gammaNew;
        ObjFuncBest=ObjFuncNew;
        if(SuccessPrevgamma==1)
            gammaSuccessCount=gammaSuccessCount+1;
        end
        SuccessPrevgamma=1;
    else
        alphaNew=alpha;
        betaNew=beta;
        gammaNew=gamma-delta_gamma;
        ZnNew(1:Ndata)=Zn(1:Ndata).*(alphaNew+betaNew*(Zn(1:Ndata).^2-1)+gammaNew*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
        Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
        Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
        Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;

        ObjFuncNew=(Moment2-1).^2*w1+(Moment4-3).^2*w2+(Moment6-15).^2;
        if(ObjFuncNew<ObjFuncBest)
        gamma=gammaNew;
        ObjFuncBest=ObjFuncNew;
        if(SuccessPrevgamma==1)
            gammaSuccessCount=gammaSuccessCount+1;
        end
        SuccessPrevgamma=1;
        else
            delta_gamma=delta_gamma*mul1;
            SuccessPrevgamma=0;
            gammaSuccessCount=0;
        end
    end
    if(alphaSuccessCount>4)
        delta_alpha=delta_alpha*mul2;
    end
    if(betaSuccessCount>4)
        delta_beta=delta_beta*mul2;
    end
    if(gammaSuccessCount>4)
        delta_gamma=delta_gamma*mul2;
    end
end
ZnNew(1:Ndata)=Zn(1:Ndata).*(alpha+beta*(Zn(1:Ndata).^2-1)+gamma*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
        
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;        


Moment2
Moment4
Moment6

%str=input("Look at moment values")

end

.
.
.
I have in the copied post talked about moment matching for monte carlo where underlying probability of every path is equal and uniform. But this new post deals with cases when underlying probability is known with certainty and is not uniform or symmetric around median.

Another point under consideration is that I have used moment matching for regression that calculates Z-series of coefficients of conditional Z-series of asset density for every volatility grid cell. But our underlying stochastic volatility grid till now has been symmetric around median. Once we make it proportional to local change in parameters, the underlying volatility grid will not remain symmetric any more. For moment matching of a symmetric grid, we needed to match only three parameters(through hermite series multiplication) in our numerical iterative moment matching function but in case of a completely general grid, we will have to match six parameters in our grid so there will be some more work required to properly calculate everything.

Amin · November 11th, 2023, 11:52 am

Friends, I have been very careful about food and trying to avoid anything that can possibly be drugged and therefore have been able to survive and do decent research. But mind control agents are extremely upset that they failed to retard me despite drugging food in huge swathes of Lahore city and now they have resorted to other tactics.
I had told friends that they were targeting my teeth with electromagnetic waves. I would have pain in all teeth but one of the teeth that was slightly weaker would have huge pain. However they ended focusing EM waves on my teeth later and the pain ended in all of the teeth. But after failing to retard me in past few weeks, they have started targeting my teeth again with EM waves and there is huge pain in the teeth again. Mind control agents are extremely upset that I am trying to do good research and want to stop me at all cost.

Amin · November 11th, 2023, 6:28 pm

Friends, I updated the new grid generation algorithm and now it can generate systematic non-symmetric probability mass based grids. I have also modified the moment matching program to deal with moment matching of non-symmetric Z-based grids for regression. This new optimization program modifies the Z-array associated with the probability mass of un-symmetric grid by optimizing six parameters and works well.
I will be playing with the program for a day or two and then post the new program again on Wilmott.
We are getting closer to end of the first chapter of our research and will move on to more related and more interesting research work in about a week.
We have to tackle the correlated SV SDEs case and we also want to work with modelling the future distribution of asset conditional on both asset and volatility at some earlier point in time.

Amin · November 12th, 2023, 6:11 am

Friends, I have been thinking of modelling of correlated SV SDE simulated stochastic variables and other cases when two variables are correlated. I think I have good ideas about how to tackle the correlated case.
First of all, we have to find univariate Z-series of stochastic volatility and asset and represent them by respective hermite polynomial series. Then we find covariances between hermite polynomials of different order between the both univariate hermite polynomial series of asset and stochastic volatility. After finding covariace/correlations, we separate the asset into a Z-series part correlated with stochastic volatility and another Z-series part uncorrelated with stochastic volatility. We model the uncorrelated part of the asset just as we have been doing in our experiments for past few weeks. And model the correlated part of the asset as a pure function of stochastic volatility in terms of stochastic volatility Z-series.
We can apply the same general procedure for modelling of stochastic variables other than SDE simulated data.
I will try to explain the above procedure in terms of equations and in more detail later. I hope that we will have a working program ready for dealing with modelling of monte carlo simulated correlated stochastic variables in a few days.

Amin · November 16th, 2023, 7:07 am

Friends, for past few days, I tried to tackle the correlated SDEs case but I have not been successful yet. I am still working on it and hope that we will soon have a well defined strategy to deal with various types of data.
Meanwhile, I had other ideas about dealing with the 2D case. I believe that we can deal with two dimensional random variable case in a single weighted least squares regression. And similarly for three dimensional and higher dimensional random variables. I hope to present a program by the end of this weekend that solves the 2D density of asset and stochastic volatility problem with a single weighted least squares regression.
Sometimes, I think that we just need to increase dimensionality to deal with correlated case and I will have to check this possibility.
Please not that cross-sectional correlation that we have in our problem is different from correlation in time that is used to advance monte carlo simulations.

Meanwhile, my persecution by drugging the city continues unabated. They have distributed mind control drugs among a very large number of street restaurants and they keep adding them to their food. Day before Yesterday, I bought a very small amount of cooked food and had a complete hangover yesterday after having food.
Mind control agents are so desperate that they try to stop shopkeepers from selling me food when they notice that I am buying good food. They project their voice on shopkeepers through voice to skull and try to persuade the shopkeepers to not sell me good food. Some shopkeepers are nice and ignore the voice to skull messages but others become belligerent to not sell their food items. Day before yesterday, I went to a bakery in a remote far off area of the city. It had a brand of milk that has remained good for a long time. I asked the shopkeeper to give me a few packets of milk. He gave me the milk and I looked carefully at the production date of the milk which was in September 2023 and the expiry of milk was end December 2023. My eyes have become weak and it took me sometime to read the production date. In the meantime, mind control agents asked the shopkeeper over voice to skull to not sell me the milk. He had already placed milk tetra packs on the counter but he told me that he would not sell the milk as it was expired. I told him that he could verify on the tetra packs that milk would expire on 29th Dec 2023. He argued with me for a bit but then agreed that milk was not expired but still took the milk packets off the counter and placed them on the aisle saying that his partner would come in a few minutes and only he could give the milk to me. He totally refused to sell the milk to me until his partner would come.
And this is not a single isolated incident. Such incidents of people trying to not sell their food happen almost everyday even though some people really do not care about voice to skull and remain very nice.
Another similar thing. Army has already approached almost all pharmacies in Lahore and have asked them to obey whenever they would send voice to skull messages to pharmacy staff. I buy an antipsychotic injection after a month and also buy some vitamins and minerals after a month. And many of the pharmacy staff refuse that they have the vitamins/minerals I ask for when voice to skull asks them to not sell anything of the sort related to good vitamins/minerals. I try to take mineral supplements regularly since I think that many of the rare elements that are in our body in only trace amounts are used in our brain to make special neurotransmitters and I will also like to suggest to people who do a lot of mental and creative work that they should take some good mineral supplements for better neurotransmitters(even though there are many research reports that mineral and vitamin supplements might have no special benefit for everyone). And mind control agents have a very severe problem with my vitamin/mineral supplement. I would like to tell friends at this time that a large number of vitamin and mineral supplements in the pharmacies especially by big multinational drug companies and their associates in Pakistan are added with mind control chemicals at the time of manufacturing as ordered by army. But some lesser known brands of vitamin and mineral supplements remain good and do not have mind control chemicals in them and I have been able to buy them for past few months.

Amin · November 17th, 2023, 2:32 pm

Friends, when I returned home yesterday after buying food in the city, I was feeling very cold and my body has been breaking (has aches all over) for past few days. I never made a big deal and thought it is minor cold but it was more severe today. I slept relatively early yesterday and woke up at three in the morning to work but was too tired to do anything. In the morning, I checked my temperature and found that I had 101.2 F fever. I took two panadol tablets and had some rest and all fever was gone in a few hours and I felt really sturdy. I wrote core of the new program on matlab program to solve the 2D interpolation problem in one single regression in the afternoon. I had thought that my fever was gone but I started to feel cold now in the night and found out that I had 101 F fever again. I had hoped to complete the program and make it perfectly running tonight but now will have to wait till tomorrow at the earliest.
Meanwhile mind control agents continue to cause strong pain in my heart and right now I am reeling from this pain.
I am slightly afraid that I might have dengue fever since I had been working in the open at night and in the evening regularly. I had dengue fever about two/three years ago and it was extremely exhausting right from the start and it was not possible to do anything other try to rest in the bed but I really have not been feeling that bad at all despite that I have aches in the back and in the body for past few days. So I am not sure about whether it is just cold or dengue fever.
I hope I recover by tomorrow and post the new matlab program in Wilmott. I am not feeling extremely bad but I have aches all over and sometimes I feel cold.

Amin · November 17th, 2023, 3:18 pm

Even though I have fever and I am not doing anything and just trying to sleep, mind control agents continue to highly charge my both hands like electrodes and I have very awkward feelings in both hands.
I am lying on a cot in my room and I had placed a side table by the bed/cot with my computer on the side table and had made my previous post lying in the bed. I had a very slight sleep after that but woke up and felt that my both hands were getting charged like electrodes with pain in both hands. Since my laptop was by my side, I decided to make this post while still in the bed.

Amin · November 19th, 2023, 6:02 am

Friends, the new method I had talked before again required a 2D grid. It was supposed to solve the 2D interpolation problem in a single WLS regression. But even this method required to make a 2D grid first. But I have been thinking that creating a 2D/3D grid is a major difficulty and might not be easily possible when data is not too large. So I have been thinking hard about any possible methods in which we do not have to construct a grid. I have thought of some ideas and want to experiment with them. I am quite confident that I might be able to solve the 2D/3D problems without construction of an underlying 2D/3D grid and we might have to use joint data with 2D or 3D arrays. I hope to be coming back on wilmott in a few days with a lot of matlab stuff.

Regarding other activities. Yesterday, I went out to get food when it was just early afternoon and it was warm and sunny. I was thinking that I had been going out to get food in earlier days after sunset and it is much colder at that time and therefore I caught cold as I was taking it too easy and was not even wearing a sweater or jacket. Since I had just gotten out of fever, I decided to go out when it was warm and sunny. I also bought a very nice and surprisingly cheap sweater-jacket (It cost me slightly upwards of four dollars). Most of my life, I had been buying brand clothes but I know now how less affluent people buy very nice clothes.
Earlier, I would go out to get food everyday after sunset but for past ten days, I have been going out to get food every second day. It would take four to five hours everyday and it was a huge drain on my time that I needed to use towards my research. I also realized that I would spend less money when I go out every other day. So I remain mostly at home in between second day motor cycle trips other than when I have to take my mother somewhere for her class or for her shopping on the car which is considered my responsibility at home.
I have much less money left and I just spend money on food. My average expense on food is less than 14/15 dollars after every two days or seven dollars per day. And I hope that at this rate of spending, I can last till mid of February next year and then I would have to approach people with request for possible consulting projects in order to make money. Ironic that American defense is willing to spend hundreds of millions of dollars every year to drug the entire Lahore city to keep tabs on my research as racist pigs of American defense and some other insecure people get slighted when they see my doing good research. It is difficult to not get flattered given the grandiose scale of money spent by American defense to retard me.
What hell can I, who has a mere few hundred dollars, let loose that American defense has to spend hundreds of millions of dollars every year to drug Lahore city in order to retard me.
Even though I have a Muslim heritage, I do not believe in any truth about religion(s). Of course, due to my muslim background, I would like muslims to be educated, tolerant, liberal and embracing people of all other nations and religions with love. I would want them to have properties of good and nice human beings. I would want them to have a soft heart that cannot allow cruelty or wrong towards any other human being ever. But I wish the same for all people of other religions and nations and all human beings.
However, all other things apart, only important thing for me, at the moment, is to do as good research as possible till the mid of February and think of no other thing. I will reassess everything after that to see what to do but of course, my greatest priority is to continue doing good research and keep distributing it and I hope that I will get some decent consulting projects to support me while continuing my research into the future.

Amin · November 19th, 2023, 8:35 pm

Friends, I have been trying to make the single regression idea for modelling the 2D densities work but a single regression framework does not seem to work very well. I will continue to see if I can possibly alter some part to make it work but it seems a bit difficult at the moment.
I have not been able to work well for past few days. Earlier I had bad food about six days ago that was drugged and then I had body aches and other problems. I think later I caught cold and had fever also due to weakness and aches that were induced by bad drugged food. I am really good now finally but I could not properly work for many days and really lament the lost time.
Though I could not properly do the single regression for solution of 2D densities, I will quickly try to write a program that finds the 2D Z-series from 1st six and eight moments and all cross-moments of the same order. This program will not take the raw data as input but rather only moments and cross-moments. Given the case with single dimensional Z-series, similarly I do not expect that program will universally work without exception but I hope that it will be very useful for well defined densities. I hope that this new program will complement the analytic regression based framework for representing the 2D density.
Simultaneously, I also want to write a small function that calculates the graphs and densities of Z-series that have a negative first derivative dX/dZ of Z-series anywhere single or multiple times(when density folds back onto itself in previous graphing routine and the new program will just show one continuous density). I have worked on the logic of the program and would put it together quickly soon.
I also want to quickly write a routine for positive densities that have a small tail in negative X-axis to numerically convert them into an all positive densities with unit mass and as close to true positive densities as possible. This obviously will not be a Z-series representation and it will be a final numerical correction to Z-series density to convert it analytically into an alll positive density.
And I hope that I will continue the momentum of the 2D interpolation of density project and its applications to option pricing.

Amin · November 21st, 2023, 11:45 am

Friends, I had been trying to write optimization program for 2D Z-series but I wandered and had random thoughts and I realized that actual problem is much simpler. And I would probably not need the 2D optimization program that I was writing.
Once we want to find 2D Z-series of two random variables from their moments and cross-moments, We actually need to find their univariate Hermite series only. And we need to specify correlations between all hermite polynomials (of the same order across both random variables) up to a high degree(for example 20-30 or higher). Again the thesis is that once univariate Z-series of two random variables are known, all cross-moments can be pinned down if correlations between all hermite polynomials of two random variables till a high enough degree is known.
This is because powers of hermite polynomials can be re-written and replaced simply by a local hermite series(2nd power of 2nd hermite can be substituted by a hermite series of three terms in 4th, 2nd and zeroth hermite polynomials. Similarly for all other powers of various hermites) . We can expand cross moments from multiplication of powers of two univariate hermite series of both random variables and then substitute equivalent local hermite series in place of powers of hermite polynomials everywhere. These local hermite series have first powers of various hermite polynomials only and then we have an expression of cross-moments in terms of products of (first powers of) hermite polynomials of two random variables. Again we are trying to find cross moments between two random variables by taking product of appropriate powers of individual Hermite series of two random variables. Only free variables that can match this expression of product of powers of two univariate Z-series with given values of relevant cross-moments(found from data) are correlations across various hermite polynomials across two random variables since we have converted all the cross moment product expression into terms containing only products of first powers of different hermite polynomials(hence the correlations where applicable).
We can also more easily try to find these correlations from cross-moments data by optimization since correlations are only first order variables and are not found in the form of powers in our setup.
Again the thesis is that once individual Z-series/Hermite Series of two random variables are known, all cross-moments can be pinned down if correlations between hermite polynomials of two random variables till a high enough degree is known.
For friends, who have never followed previous discussions, correlations between hermite polynomials of two random variables are across the polynomials of same degree. first/nth hermite polynomial of one random variable can be correlated with first/nth hermite polynomial of some other random variable.
There might be wild inaccuracies in calculation of these correlations in practice but that is not the topic of this post and I hope we will continue to see what can be done.

Giving a simple example, If univariate Hermite-series of two different variables are expanded to fifth hermite term, the calcuation of E[X^3Y^3] can require correlations between up to fifteenth hermite polynomials between two random variables. When cross moment E[X^3Y^3] is calculated from data, we can possibly use it along with other cross moments to back out correlations between hermites by optimization since correlations are just first order variables.

I will try to explain this with Latex Equations later tonight. Obviously this result can be generalized to higher dimensions also.

Amin · November 21st, 2023, 12:41 pm

Again a 2D Hermite Series/Z Series of product of two random variables is merely a product of two univariate Hermite series of individual random variables but we have to do accounting for correlations between hemite polynomials of various orders. And in order to do that we have to represent the more complicated raw "product of two univariate Hermite-series expression" by a perfectly (algebraically) equivalent expression with terms containing only products of first powers of hermite polynomials of various orders and then we can apply the hermite correlations if known.

Similarly for addition of two random variables when their Hermite series are given, I expect that if all hermite polynomials across two random variables are uncorrelated we would add them by squared addition of coefficients and then taking square root. However if a correlation exists between two variables, they could be added by squared addition of uncorrelated parts of hermite polynomial coefficients while correlated parts of coefficients would be added linearly.

Amin · November 21st, 2023, 2:40 pm

Friends, I slept late yesterday night and woke up early this morning. I am feeling a bit off and therefore going to sleep in a bit and will make a formal post with latex equations tomorrow. I am sure friends who have been following my research would have surely found my previous posts helpful and already got the idea. Sorry for this delay.

Amin · November 22nd, 2023, 3:45 pm

Friends, I had started to believe that all magic is in correlations between various hermites. In order to test this idea, I wrote a small matlab program that evolves stochastic volatility and asset in a similar framework that we had been using. I calculated hermite series of both asset and SV and added them analytically (for that analytic part I calculated correlation between hermite polynomials of two variables from numerical data). Then I numerically added the SV and asset data from paths and made its Hermite series and plotted both analytical and numerical hermite series of the sum variable and the match between two densities was remarkable. I noticed that sometimes for extreme SDE parameters the method does not work well and I think the reason was that I was using only hermite series till 5th hermite and for accuracy in such cases we might need a 7th order hermite series.
When I used asset as base variable and added correlated part of stochastic volatility linearly to it before squaring, the method did not work well. But when I took volatility as base variable and added correlated part of asset linearly to it, it worked extremely well. So I believe some work has to be done in order to know which variable to add linearly to other variable. It was more obvious in our case since SV was independent variable(in our traditional sense).
I simply found correlations between hermite polynomials of both hermite series and added correlated part of asset hermite polynomials linearly to volatility hermite polynomials and then squared it and added it to square of uncorrelated part of the asset hermite polynomials.
Hermite polynomials are added as

[$]H_{S,0} \, = H_{V,0} \, + \, H_{X,0}[$]
[$]H_{S,1} \, = \sqrt{ {\big[H_{V,1} \, + \rho_1 \, H_{X,1}\big]}^2+ \, (1- \, {\rho_1}^2) \, {H_{X,1}}^2 }[$]
[$]H_{S,2} \, = \sqrt{ {\big[H_{V,2} \, + \rho_2 \, H_{X,2}\big]}^2+ \, (1- \, {\rho_2}^2) \, {H_{X,2}}^2 }[$]
[$]H_{S,3} \, = \sqrt{ {\big[H_{V,3} \, + \rho_3 \, H_{X,3}\big]}^2+ \, (1- \, {\rho_3}^2) \, {H_{X,3}}^2 }[$]
[$]H_{S,4} \, = \sqrt{ {\big[H_{V,4} \, + \rho_4 \, H_{X,4}\big]}^2+ \, (1- \, {\rho_4}^2) \, {H_{X,4}}^2 }[$]
[$]H_{S,5} \, = \sqrt{ {\big[H_{V,5} \, + \rho_5 \, H_{X,5}\big]}^2+ \, (1- \, {\rho_5}^2) \, {H_{X,5}}^2 }[$]

Here[$]H_{S,n}[$] is nth hermite of the sum random variable, [$]H_{V,n}[$] is nth hermite of stochastic volatility and [$]H_{X,n}[$] is the nth hermite of asset. It is SV and asset that are being added to find the sum variable.

I will post my simulation program for above within one and a half hour and hope that it will help in better understanding of the method.

Amin · November 22nd, 2023, 6:13 pm

Friends, I am posting the program in the hope that it would be helpful for friends and they could just try their version of equations for summation without having to go through writing all the basic program.
Sorry for breaking the excitement but it seems that more work needs to be done to settle the problem and I realized that a lot of scenarios are not working very well. I am trying to concentrate on this problem to solve the addition and products with correlations and will try to see if I can think of anything interesting.

Case for correlation between asset and SV needs to be more deeply considered. One reason for lack of good program might be that I am using almost lognormal asset and almost lognormal volatility and even fore moderate volatilities, they cannot be described by five hermite polynomials. Please try some smaller CEV indexes. Though I believe more serious work has to be done on the problem.

Here is the program. You might have to look for the thread for some basic functions that we have been using all the time.

function [] = AnalyticAddTwoSDEDensities()

X0=1.0;
alpha1=0;
alpha2=1;
kappaX=0;
thetaX=1;
a=kappaX*thetaX;
b=-kappaX;

rho=0;
sigma1=.65;
gammaV=.5;
gammaX=.95;
%%%%%%%%`%%%%%
V0=.950;
thetaV=.75;
kappaV=.650;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
gamma=.95;
sigma0=.65;


dt=.03125/2;
Tt=64*1;

T=Tt*dt;




%%%%%%%%%%%%%%%%%%%%%%%%5
rng(82139475, 'twister')
paths=500000;
V(1:paths)=V0;  %Original process monte carlo.
X(1:paths)=X0;
% paths0=paths/2;
% paths00=paths/4;
% Random1(1:paths00)=0;
% Random2(1:paths0)=0;


paths0=paths;
paths00=paths;
Random1(1:paths00)=0;
Random2(1:paths0)=0;

for tt=1:Tt
    
Random2=randn(size(Random2)); 
Random1=randn(size(Random1)); 

% sum(Random1)/paths
% sum(Random1.^2)/paths
% sum(Random1.^3)/paths
% sum(Random1.^4)/paths
% sum(Random1.^5)/paths
% sum(Random1.^6)/paths
% 
% tt
% 
% str=input("Look at moments of MC Normals");

%Random1(paths00+1:2*paths00)=-1*Random1(1:paths00);
%Random2(paths0+1:2*paths0)=-1*Random2(1:paths0);    
%Random1(paths0+1:2*paths0)=Random1(1:paths0);    
%antithetic variates above


X(1:paths)=X(1:paths)+ ...
    (a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dt + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dt) + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*sqrt(dt) + ...
    (a*alpha1* X(1:paths).^(alpha1-1)+b*alpha2* X(1:paths).^(alpha2-1)).* ...
    (((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dt^2/2)+ ...
    (sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) *(1-1/sqrt(3)).*dt^1.5+ ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*(1-1/sqrt(3)).*dt^1.5))+ ...
    .5*(a*alpha1*(alpha1-1)* X(1:paths).^(alpha1-2)+b*alpha2*(alpha2-1).* X(1:paths).^(alpha2-2)).* ...
    ( sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX)) *dt^2/2 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
    ((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random1(1:paths) * 1/sqrt(3).* dt^1.5 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dt/2 + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths)*dt/2)+ ...
    .5*sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
    (sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dt^1.5)+ ...
    sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
    ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random1(1:paths) * 1/sqrt(3).* dt^1.5 + ...
    sigma0*V(1:paths).^gamma.*Random1(1:paths).*Random2(1:paths)*dt/2)+ ...
    .5*sqrt(1-rho^2)* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
    (sigma0^2*V(1:paths).^(2*gamma).*Random1(1:paths)*1/sqrt(3)*dt^1.5)+ ...
    sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
    rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random1(1:paths)*1/sqrt(3)*dt^1.5+ ...
    rho* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
    ((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random2(1:paths) * 1/sqrt(3).* dt^1.5 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths) * dt/2 + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*(Random2(1:paths).^2-1)*dt/2)+ ...
    .5*rho* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
    (sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random2(1:paths).*1/sqrt(3).*dt^1.5)+ ...
    rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
    ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths) * 1/sqrt(3).* dt^1.5 + ...
    sigma0*V(1:paths).^gamma.*(Random2(1:paths).^2-1)*dt/2)+ ...
    .5*rho* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
    sigma0^2.*V(1:paths).^(2*gamma).*Random2(1:paths) * 1/sqrt(3).* dt^1.5+ ...
    rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
    rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random2(1:paths)*1/sqrt(3)*dt^1.5;
    

    V(1:paths)=V(1:paths)+ ...
        (mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dt + ...
        sigma0*V(1:paths).^gamma .*Random2(1:paths)*sqrt(dt) + ...
        (mu1.*beta1*V(1:paths).^(beta1-1) + mu2.*beta2.*V(1:paths).^(beta2-1)).* ...
        ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dt^2/2 + ...
        sigma0*V(1:paths).^gamma .*Random2(1:paths)*(1-1/sqrt(3))*dt^1.5) + ...
        .5*(mu1.*beta1.*(beta1-1).*V(1:paths).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths).^(beta2-2)).* ...    
        sigma0^2.*V(1:paths).^(2*gamma).*dt^2/2 + ...
        sigma0*gamma*V(1:paths).^(gamma-1) .* ...
        ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths).*1/sqrt(3)*dt^1.5 + ...
        sigma0.*V(1:paths).^gamma .*(Random2(1:paths).^2-1)*dt/2) + ...
        .5*sigma0*gamma*(gamma-1).*V(1:paths).^(gamma-2) .* ...
        sigma0^2.*V(1:paths).^(2*gamma) .*Random2(1:paths).*1/sqrt(3)*dt^1.5;


end
% X1=X;
% V1=V;
% 
% save('MonteCarloData1','X1','V1');


%str=input("Monte carlo data has been saved");
X(X<0)=0;
V(V<0)=0;
X=real(X);
V=real(V);

% SVolMeanAnalytic=thetaV+exp(-kappaV*T)*(V0-thetaV)
% SVolMean=sum(V(1:paths))/paths
% AssetMeanAnalytic=thetaX+exp(-kappaX*T)*(X0-thetaX)
% AssetMean=sum(X(1:paths))/paths

    W=X+V; %Add asset and volatility to find sum variable

    
    
    
    Ndata=paths;
    Mtable=readtable('C:\Users\Lenovo\Documents\MATLAB\MATLAB1\MomentMatchParams08.txt');
    [Xs,Ix]=sort(X);
    [ZI] = MatchMomentParametersOfIdealZhalf04(Ndata,Mtable);
    ZProb(1:Ndata)=1/Ndata;
    [a0,a] = CalculateZSeriesDensityByIdealMeasureMean05WithZI(Xs,ZI,ZProb);
    [ah0,ah] = ConvertZCoeffsToHCoeffs(a0,a,5);
    
    %Above ah0, ah are coefficients of hermite series of asset. a0,a are
    %coefficients of Z-series of assset.

    
    [Vs,Iv]=sort(V);
    %[ZI] = MatchMomentParametersOfIdealZhalf04(Ndata,Mtable);
    ZProb(1:Ndata)=1/Ndata;
    [b0,b] = CalculateZSeriesDensityByIdealMeasureMean05WithZI(Vs,ZI,ZProb);
    [bh0,bh] = ConvertZCoeffsToHCoeffs(b0,b,5);
    
     %Above bh0, bh are coefficients of hermite series of SV. b0,b are
    %coefficients of Z-series of SV.
    
    
    
    [Zx0] = CalculateZgivenXAndZSeriesC5Improved(X,a0,a);
    [Zv0] = CalculateZgivenXAndZSeriesC5Improved(V,b0,b);
    
    for pp=1:paths
       Zx1(Ix(pp))=norminv(pp/paths-.000001*pp/paths);
       Zv1(Iv(pp))=norminv(pp/paths-.000001*pp/paths);
    end
    
    Nhermites=5;
    [CorrH1] = CalculateCorrelationsBetweenTwoZSeriesRVs(Zx0,Zv0,paths,Nhermites);
    
    [CorrH] = CalculateCorrelationsBetweenTwoZSeriesRVs(Zx1,Zv1,paths,Nhermites);
    
    CorrH
    CorrH1
    str=input("Look at correlations");
    
    ch0=ah0+bh0;
    
    
%     chh(1)=sqrt((ah(1)+CorrH(1).*bh(1)).^2+(1-CorrH(1).^2).*bh(1).^2);
%     chh(2)=sqrt((ah(2)+CorrH(2).*bh(2)).^2+(1-CorrH(2).^2).*bh(2).^2);
%     chh(3)=sqrt((ah(3)+CorrH(3).*bh(3)).^2+(1-CorrH(3).^2).*bh(3).^2);
%     chh(4)=sqrt((ah(4)+CorrH(4).*bh(4)).^2+(1-CorrH(4).^2).*bh(4).^2);
%     chh(5)=sqrt((ah(5)+CorrH(5).*bh(5)).^2+(1-CorrH(5).^2).*bh(5).^2);
%     
%     

%Equations for Sum random variable as posted on wilmott.
     ch(1)=sqrt((bh(1)+CorrH(1).*ah(1)).^2+(1-CorrH(1).^2).*ah(1).^2);
     ch(2)=sqrt((bh(2)+CorrH(2).*ah(2)).^2+(1-CorrH(2).^2).*ah(2).^2);
     ch(3)=sqrt((bh(3)+CorrH(3).*ah(3)).^2+(1-CorrH(3).^2).*ah(3).^2);
     ch(4)=sqrt((bh(4)+CorrH(4).*ah(4)).^2+(1-CorrH(4).^2).*ah(4).^2);
     ch(5)=sqrt((bh(5)+CorrH(5).*ah(5)).^2+(1-CorrH(5).^2).*ah(5).^2);
     
     
     
%Works very well for uncorrelated(in traditional sense) but is seemingly wrong.     
%     ch(1)=sqrt((bh(1)+CorrH(1).*ah(1)).^2+(1-CorrH(1)).^2.*ah(1).^2);
%     ch(2)=sqrt((bh(2)+CorrH(2).*ah(2)).^2+(1-CorrH(2)).^2.*ah(2).^2);
%     ch(3)=sqrt((bh(3)+CorrH(3).*ah(3)).^2+(1-CorrH(3)).^2.*ah(3).^2);
%     ch(4)=sqrt((bh(4)+CorrH(4).*ah(4)).^2+(1-CorrH(4)).^2.*ah(4).^2);
%     ch(5)=sqrt((bh(5)+CorrH(5).*ah(5)).^2+(1-CorrH(5)).^2.*ah(5).^2);
%     
    
%     ch(1)=.5*ch(1)+.5*chh(1);
%     ch(2)=.5*ch(2)+.5*chh(2);
%     ch(3)=.5*ch(3)+.5*chh(3);
%     ch(4)=.5*ch(4)+.5*chh(4);
%     ch(5)=.5*ch(5)+.5*chh(5);
        
    
    [c0,c] = ConvertHCoeffsToZCoeffs(ch0,ch,5);
    
    [Ws,Iw]=sort(W);
    [d0,d] = CalculateZSeriesDensityByIdealMeasureMean05WithZI(Ws,ZI,ZProb);
    [dh0,dh] = ConvertHCoeffsToZCoeffs(d0,d,5);

    
    ah0
    ah
    bh0
    bh
    ch0
    ch
    dh0
    dh
    
    a0
    a
    b0
    b
    d0
    c0
    d
    c
    

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


XX(1:Nn)=c0;
for nn=1:5
    XX(1:Nn)=XX(1:Nn)+c(nn)*Z(1:Nn).^nn;
end

XX0(1:Nn)=d0;
for nn=1:5
    XX0(1:Nn)=XX0(1:Nn)+d(nn)*Z(1:Nn).^nn;
end

clf;

plot(Z(1:Nn),XX(1:Nn),'r',Z(1:Nn),XX0(1:Nn),'g')
    
    
    
    str=input("Look at coeffs of Added variable and Graph comparison")
    
    
    clf;
    PlotZSeriesDensity(d0,d,'g');
    hold on
    PlotZSeriesDensity(c0,c,'r');
    hold off
    str=input("Look at densities Graph comparison")
    
    
end

.
.
.

function [ZnNew] = MatchMomentParametersOfIdealZhalf04(Ndata,Mtable)
[Zn] = StandardNormalIdealDesity02(Ndata,1);


[alpha0,beta0,gamma0] = ReadParamsFile(Ndata,Mtable)



ZnNew(1:Ndata)=Zn(1:Ndata).*(alpha0+beta0*(Zn(1:Ndata).^2-1)+gamma0*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
        
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;        


Moment2
Moment4
Moment6
display("fitted Moments");
%str=input("Look at moment values")

end

.
.
.

function [a0,a] = CalculateZSeriesDensityByIdealMeasureMean05WithZI(data,ZI,ZProb)


Ndata=length(data);



Sdata=sort(data);


Y(1:Ndata,1)=Sdata(1:Ndata);
Mean=sum(Y(1:Ndata,1))/Ndata;
Y(1:Ndata,1)=Y(1:Ndata,1)-Mean; %subtract the mean. Variance standardization might not have an effect.
if (rem(Ndata,2)==0)
%Median1=sum(Y(Ndata/2,1)+Y(Ndata/2+1,1))/2;

x0=0;
x1=ZI(Ndata/2-3);
x2=ZI(Ndata/2-2);
x3=ZI(Ndata/2-1);
x4=ZI(Ndata/2+1);
x5=ZI(Ndata/2+2);
x6=ZI(Ndata/2+3);

y1=Y(Ndata/2-3);
y2=Y(Ndata/2-2);
y3=Y(Ndata/2-1);
y4=Y(Ndata/2+1);
y5=Y(Ndata/2+2);
y6=Y(Ndata/2+3);
N=6;
[y0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

Median=y0;

% Median1
% Median
% %str=input("Look at Medians");
else
Median=Y(floor(Ndata/2+1),1);
end

Ymu(1:Ndata,1)=Y(1:Ndata,1)-Median;%subtract the median
Zu=0;   %Since median is subtracted from mean-standardized data, we have to 
%subtract the Z=0 from ZI i.e subtract nothing.
X(1:Ndata,1:6)=0;
%Below I am regressing powers of Z as in Z-series but hermites could
%equally be regressed and I expect that results would probably be extremely
%close.
for nn=1:Ndata
    X(nn,1)=1;
    X(nn,2)=(ZI(nn)-Zu);
    X(nn,3)=(ZI(nn)-Zu).^2;
    X(nn,4)=(ZI(nn)-Zu).^3;
    X(nn,5)=(ZI(nn)-Zu).^4;
    X(nn,6)=(ZI(nn)-Zu).^5;
    %X(nn,7)=(ZI(nn)-Zu).^6;
    %X(nn,8)=(ZI(nn)-Zu).^7;
end

%coeff=inv(X'*X)*(X'*Ymu);
coeff=(X'*X)\(X'*Ymu);
c0=coeff(1,1);
c(1:5)=coeff(2:6,1);

a0=c0;
a=c;
a0=a0+Median+Mean;
end

.
.
.

function [CorrH] = CalculateCorrelationsBetweenTwoZSeriesRVs(Z1,Z2,paths,Nhermites)

CorrH(1:Nhermites)=0;
CovH(1:Nhermites)=0;
VarH1(1:Nhermites)=0;
VarH2(1:Nhermites)=0;
for pp=1:paths
    for hh=1:Nhermites
        if(hh==1)
            CovH(hh)=CovH(hh)+Z1(pp).*Z2(pp)/paths;
            VarH1(hh)=VarH1(hh)+Z1(pp).*Z1(pp)/paths;
            VarH2(hh)=VarH2(hh)+Z2(pp).*Z2(pp)/paths;
        end
        if(hh==2)
            CovH(hh)=CovH(hh)+(Z1(pp).^2-1).*(Z2(pp).^2-1)/paths;
            VarH1(hh)=VarH1(hh)+(Z1(pp).^2-1).*(Z1(pp).^2-1)/paths;
            VarH2(hh)=VarH2(hh)+(Z2(pp).^2-1).*(Z2(pp).^2-1)/paths;
        end
        if(hh==3)
            CovH(hh)=CovH(hh)+(Z1(pp).^3-3*Z1(pp)).*(Z2(pp).^3-3*Z2(pp))/paths;
            VarH1(hh)=VarH1(hh)+(Z1(pp).^3-3*Z1(pp)).*(Z1(pp).^3-3*Z1(pp))/paths;
            VarH2(hh)=VarH2(hh)+(Z2(pp).^3-3*Z2(pp)).*(Z2(pp).^3-3*Z2(pp))/paths;
        end
        if(hh==4)
            CovH(hh)=CovH(hh)+(Z1(pp).^4-6*Z1(pp).^2+3).*(Z2(pp).^4-6*Z2(pp).^2+3)/paths;
            VarH1(hh)=VarH1(hh)+(Z1(pp).^4-6*Z1(pp).^2+3).*(Z1(pp).^4-6*Z1(pp).^2+3)/paths;
            VarH2(hh)=VarH2(hh)+(Z2(pp).^4-6*Z2(pp).^2+3).*(Z2(pp).^4-6*Z2(pp).^2+3)/paths;
        end
        if(hh==5)
            CovH(hh)=CovH(hh)+(Z1(pp).^5-10*Z1(pp).^3+15*Z1(pp)).*(Z2(pp).^5-10*Z2(pp).^3+15*Z2(pp))/paths;
            VarH1(hh)=VarH1(hh)+(Z1(pp).^5-10*Z1(pp).^3+15*Z1(pp)).*(Z1(pp).^5-10*Z1(pp).^3+15*Z1(pp))/paths;
            VarH2(hh)=VarH2(hh)+(Z2(pp).^5-10*Z2(pp).^3+15*Z2(pp)).*(Z2(pp).^5-10*Z2(pp).^3+15*Z2(pp))/paths;
        end
    end
end


for hh=1:Nhermites
     CorrH(hh)=CovH(hh)./sqrt(VarH1(hh))./sqrt(VarH2(hh));
%      if(hh==1)
%          CorrH(hh)=CovH(hh);
%      end
%      if(hh==2)
%          CorrH(hh)=CovH(hh)/2;;
%      end
%     if(hh==3)
%          CorrH(hh)=CovH(hh)/6;
%     end
%     if(hh==4)
%          CorrH(hh)=CovH(hh)/24;
%     end
%     if(hh==5)
%          CorrH(hh)=CovH(hh)/120;
%      end 
     
end

VarH1
VarH2
CovH
CorrH

end

.
.
.

function [] = PlotZSeriesDensity(c0,c,Color)

Index=length(c);

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);
%Z
%str=input('Look at Z');
%Now calculate Y as a function of normal random variable.
Y(1:Nn)=c0;
for nn=1:Index
    Y(1:Nn)=Y(1:Nn)+c(nn)*Z(1:Nn).^nn;
end

%Now take change of densities derivative of Y with respect to normal
DfY(1:Nn)=0;
for nn=2:Nn-1
    DfY(nn) = (Y(nn + 1) - Y(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    %Change of variable derivative for densities
end
DfY(Nn)=DfY(Nn-1);
DfY(1)=DfY(2);

%Now calculate the density of Y from density of normal random variable
%using change of probability derivative.
pY(1:Nn)=0;
for nn = 1:Nn    
    pY(nn) = (normpdf(Z(nn),0, 1))/abs(DfY(nn));
end
Y;
%clf;
plot(Y(1:Nn),pY(1:Nn),Color);
title('Density of the Z-Series Variable with Fitted Moments');
str=input('Look at density');

% %clf;
% plot(Z(1:Nn),Y(1:Nn),Color);
% title('Graph of Z-Series Variable as a function of underlying standard normal');
% str=input('Look at Graph of Z-Series Variable as a function of underlying standard normal');
% 

end

.
.
.

Breakthrough in the theory of stochastic differential equations and their simulation

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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