Breakthrough in the theory of stochastic differential equations and their simulation

Amin · May 27th, 2021, 10:18 am

Friends, I had been thinking about the filtering model for statistical trading. I had a lot of interesting ideas and I want to share them with friends.
First of all the model we want to develop could be used in for vanilla trading on a sub-minute scale to a few minutes scale(or much more as we will later see.) The model could also be used in an excellent fashion in pairs and relative value trading(in scalar and vector form for a group of several pairs simultaneously with covariances involved).
1. We will start with a simple linear model, gain intuition, learn stylized facts a and then see if we need to extend it to non-linear filtering models. Here I describe the model with SDEs(or an equivalent form in discrete state space version). First we take volatility constant and suppose that drift has its own dynamics dictated by another SDE.
2. Since we are dealing with a level of minutes at maximum, we will not consider the asset dynamics as lognormal and just consider its noise as normal scaled by volatility. Our discrete time spacing can range from five seconds to twenty seconds for vanilla trading. At this level, considering a lognormal SDE would be a bad choice.
3. since we are working/trading on a level of few minutes at max, units of time will not be years but rather possibly in minutes or hours.
4. In order to decrease noise due to bid-ask crossing, we will work with mid-prices.
5. drift cannot continue to increase so it must be slowly mean-reverting.

Here is the SDE of financial asset or pair ratio
[$]X(t+1)=X(t) \,+ \mu(t) \, dt + \sigma \, dz(t)[$]
while [$]\mu(t)[$] has its SDE given as
[$]\mu(t+1)=\mu(t) \, - \theta \, \mu(t) \, dt+ \sigma_{mu} dw(t)[$]

Keeping the SDE for drift the same, We can write the SDE for financial asset in returns form for the asset as
[$]dX(t)=\, \mu(t) \, dt + \sigma \, dz(t)[$]

If the model(SDEs and parameters) is well-specified, its variance with drift will be smaller than its variance without drift and if model is poorly specified, its variance will be larger as compared to without drift variance.
We notice that asset price is fully observable while drift is a hidden variable and the two SDEs are exactly in Gaussian Kalman filter form. I write the two SDEs again
[$]dX(t)= \, \mu(t) \, dt + \sigma \, dz(t)[$]
while [$]\mu(t)[$] has its SDE given as
[$]\mu(t+1)=\mu(t) \, - \theta \, \mu(t) \, dt+ \sigma_{mu} dw(t)[$]

So we can use the huge literature and machinery developed for Kalman filter to solve this problem.
Returns of financial are highly notorious to defy attempts at finding structure but I am sure the above simple model would work successfully in a large number of cases in reality. With properly optimized parameters, The above model would continue to tell us when the drift is significantly different from zero and help us in trading decisions.
As I stated in my previous post, with a properly optimized model, if the level of drift constantly stays insignificant and far smaller than volatility, we can say that we are in a mean-reverting regime and we need to trade accordingly. On the other hand, when the drift continues to remain in significant range, we can make directional trades according to short-term trend dictated by direction of drift.
Since we are in a Kalman filter mode, we can make our model very interesting by making it large dimensional and by adding covariances across different financial assets and their drifts and the model will still remain robust and tractable.
I will be coming with worked out programs and more equations in a few days.
In the above post, I have hinted at high dimension Kalman filter with a hidden vector process of drift but even if we make the drift stochastic, volatility would still not be constant and would continue to remain stochastic and I was thinking of simpler ways to factor in the effect of changing volatility.
In the equation below upper case bold letters represent matrices, lower case bold letters represent vectors and non-bold lower case represent scalars. Writing the observation equation for the financial asset, we have
[$]\textbf{dx(t)} \, = \, \boldsymbol{\mu(t)} \, dt \, + \eta \, \boldsymbol{R \, z_1(t)} \, +  \, \boldsymbol{\Sigma \, z_2(t)} [$]
In hte above equation scalar [$]\eta[$] is a single dimensional hidden variable with its own hidden state space dynamics and total noise variance that enters the observation equation now becomes [$]{\eta}^2 \boldsymbol{RR' \,+\, \Sigma {\Sigma}' }[$]. The matrices [$]\boldsymbol{R}[$] and [$]\boldsymbol{\Sigma}[$] are kept constant and have been pre-calculated in the optimization phase. So we have a constant variance and there is a time changing variance that is scaled by a single hidden parameter [$]\eta[$] that has its own dynamics.
So we can write the filtering equations as
[$]\textbf{dx(t)} \, = \, \boldsymbol{\mu(t)} \, dt \, + \eta \, \boldsymbol{R \, z_1(t)} \, +  \, \boldsymbol{\Sigma \, z_2(t)} [$]

[$]\boldsymbol{\mu(t+1)}=\boldsymbol{\mu(t)} \, - \boldsymbol{\Theta \, \mu(t) }\, dt+ \boldsymbol{\sigma_{mu} \, w(t)}[$]

and scalar [$]\eta[$] process follows its own simple dynamics as

[$]\eta(t+1)=\eta(t) \, - \theta \, \eta(t) \, dt+ \sigma_{eta} v(t)[$]

due to [$]\eta[$] this filter will not be a kalman filter and will require some hack to solve this intelligently somehow. This may not be extremely difficult since we kept it as a single dimensional variable in a large vector-matrix system.

We could try other recipes for variance by choosing [$]eta[$] as the scalar magnitude of largest stochastic eigenvalue with a system of the sort like
[$]\textbf{dx(t)} \, = \, \boldsymbol{\mu(t)} \, dt \, + \eta_1 \, \boldsymbol{r_1} \, z_1(t) \,+ \eta_2 \, \boldsymbol{r_2} \, z_2(t) \, +  \, \boldsymbol{\Sigma \, z_3(t)} [$]

..
.
I could not complete the previous post so I am writing it again.
In the above post, I have hinted at high dimension Kalman filter with a hidden vector process of drift but even if we make the drift stochastic, volatility would still not be constant and would continue to remain stochastic and I was thinking of simpler ways to factor in the effect of changing volatility.
In the equation below upper case bold letters represent matrices, lower case bold letters represent vectors and non-bold lower case represent scalars. Writing the observation equation for the financial asset, we have
[$]\textbf{dx(t)} \, = \, \boldsymbol{\mu(t)} \, dt \, + \eta \, \boldsymbol{R \, z_1(t)} \, +  \, \boldsymbol{\Sigma \, z_2(t)} [$]
In hte above equation scalar [$]\eta[$] is a single dimensional hidden variable with its own hidden state space dynamics and total noise variance that enters the observation equation now becomes [$]{\eta}^2 \boldsymbol{RR' \,+\, \Sigma {\Sigma}' }[$]. The matrices [$]\boldsymbol{R}[$] and [$]\boldsymbol{\Sigma}[$] are kept constant and have been pre-calculated in the optimization phase. So we have a constant variance and there is a time changing variance that is scaled by a single hidden parameter [$]\eta[$] that has its own dynamics.
So we can write the filtering equations as
[$]\textbf{dx(t)} \, = \, \boldsymbol{\mu(t)} \, dt \, + \eta \, \boldsymbol{R \, z_1(t)} \, +  \, \boldsymbol{\Sigma \, z_2(t)} [$]

[$]\boldsymbol{\mu(t+1)}=\boldsymbol{\mu(t)} \, - \boldsymbol{\Theta \, \mu(t) }\, dt+ \boldsymbol{\sigma_{mu} \, w(t)}[$]

and scalar [$]\eta[$] process follows its own simple dynamics as

[$]\eta(t+1)=\eta(t) \, - \theta \, \eta(t) \, dt+ \sigma_{eta} v(t)[$]

due to [$]\eta[$] this filter will not be a kalman filter and will require some hack to solve this intelligently somehow. This may not be extremely difficult since we kept it as a single dimensional variable in a large vector-matrix system.

We could try other recipes for variance by choosing [$]\eta[$] as the scalar magnitude of largest stochastic eigenvalue with a system of the sort like
[$]\textbf{dx(t)} \, = \, \boldsymbol{\mu(t)} \, dt \, + \eta_1 \, \boldsymbol{r_1} \, z_1(t) \,+ \eta_2 \, \boldsymbol{r_2} \, z_2(t) \, +  \, \boldsymbol{\Sigma \, z_3(t)} [$]

here [$]\eta_1[$] is the largest eigenvalue and [$]\boldsymbol{r_1}[$] is the eigenvector associated with the variable variance that is calculated after taking into account constant variance and variance associated with the drift and similarly [$]\eta_2[$] is the second largest eigenvalue and associated second largest eigenvector is [$]\boldsymbol{r_2}[$] and The matrix [$]\boldsymbol{\Sigma}[$] is a pre-optimized matrix of coefficients associated with constant variance component.
So in both of the above scenarios, we have divided the volatility into a constant component and a stochastic component with parsimonious but hidden parameters that have to be additionally filtered. Due to hidden stochastic volatility parameters our filter is not a traditional Kalman filter but I expect that we can remain close to Kalman framework and solve for volatility using some hack.
We have to understand that our estimates of  [$]\eta_1[$], [$]\boldsymbol{r_1}[$] , [$]\eta_2[$], [$]\boldsymbol{r_2}[$] would be very different from traditional estimates of first two eigenvectors and their eigenvalues since stochastic drift would possibly take up a large part of explainable structural variance that is usually assigned to significant eigenvectors in traditional analysis without proper specification of drift vector and its dynamics. So we would have to optimize for these vector parameters and others in our filtering framework after taking into account the stochastic drift..
I have just given these examples to include stochastic volatility component in filtering in a parsimonious way but they are just meant to stimulate thought and I am sure there can be better specifications and I am sure friends would like to think about interesting models and would have many more innovative ideas.

Amin · May 28th, 2021, 11:01 am

My phone was stolen today from my car. I have written this post in the off-topic about it. Please protest to mind control agencies and ask them to end my persecution and stop any further machinations.
https://forum.wilmott.com/viewtopic.php?f=15&t=94796&p=866345#p866345

Amin · May 29th, 2021, 4:24 am

I have subscribed to CNN's "Meanwhile in America" and I was reading that new President wants to honor the memory of dead in a massacre of American blacks in Tulsa a hundred years ago. Though it is a great gesture by the new President, there are many times when I wonder why American establishment and government are so late to acknowledge the wrongs in American society and redress them.
We can recall when Martin Luther King Jr. was trying to lead the civil rights movement towards success, he was called a great traitor by American government and American intelligence agencies were trying to do everything in their power to sabotage him. Though Martin Luther King made a great impact on American society, the government apparatus remained antagonistic during his life and he was finally murdered. But It is a great thing that American society today acknowledges his contribution to American society and Martin Luther king is known as one of the great leaders in American history.
They say that the new President owes his presidency to African American voters. And I must say gain that it is a great symbolic gesture by the new President to visit the site of massacre of Black Americans. But I would like to say that is far more important (and probably less politically convenient) to bring a change in the life of living and end torture on them than honor the dead. I am trying to point towards mind control torture on thousands of intelligent blacks who have been victimized by racist hardliners in American army and mind control agencies.
Because racism has greatly decreased in American society, it is very difficult to openly continue racist practices in US society. But there is still a substantial segment of hardliners in American society who would do everything in their power to fail people of color especially when they could recognize brilliant people very early due to their neurotransmitters long before these people have distinguished themselves in American society. When many intelligent people(of color) would not have done anything extraordinary in their early life and their talent would be recognized by crooks in American army and mind control agencies, these racist crooks are quick to put such people of color and foreigners on mind control on all sort of totally false pretexts. Many large American universities proudly engage in this practice.
I want to request the new president if he really means to bring a great and rightful change in the life of American blacks, please end the mind control of thousands of intelligent blacks who are running around all over the country to avoid mind control torture. While it may not be politically expedient, it may be the perfectly right thing to do.
I had subscribed to a mind control group on yahoo (Yahoo has ended groups now) and I used to receive heartbreaking accounts of torture on this extremely nice and intelligent black woman who would post all the stories of torture every week or two on yahoo mind control group. I always knew that there is very little humanity in mind control crooks and even women of color are not safe from their inhuman animal torture since racist crooks are very sure they can get away with it.
A very large number of Americans and people all over the world rightfully praise American society for its great diversity but there are still a significant segment of population that does not like this diversity and when these racist crooks become powerful in American army and mind control agencies, they do everything in their power to fail foreigners, blacks and other people of color.
As history tell us again and again, I am very sure that evil of mind control will definitely end in American society but I will request Americans to end it before a large number of people, who are a victim of mind control torture for decades, finally have died. I am sure most people would agree that it is far better to end the wrong on living humans than celebrate the memory of dead who had been wronged.

Amin · May 29th, 2021, 12:44 pm

I want to tell friends about another machinatory thing that mind control crooks have done. I had bought ten years tick data for S&P500 two years ago till the end of 2019. I checked four or five days ago when I made the first post on this forum about kalman filtering for drift and the ten year S&P500 data was there on one of my external hard drives. But when I checked today for the data to plan everything in tackling the data step by step, I found that all data had been deleted from my external hard drive. Crooks of Pakistan army openly come to my house in my absence and do whatever they want with my stuff especially lowly things like stealing my under-wears or jeans that do not allow mind control frequencies to be effective (once they stole some of my books that I was using for research) and my family fully cooperates with them. When I complain about any such thing to my mother, she starts laughing in a goofy style and blames everything on my schizophrenia.
I still have a lot of Nasdaq ITCH file data and I have my C++ programs to construct limit order book from ITCH data so I might be able to play with Nasdaq data with my filters for now if they do not delete my ITCH files.
I want to tell friends that I have a meeting with my psychiatrist over wahtsapp around 2:00 pm on Monday. I do not have a phone at the moment so I would talk using my mother's phone or my sister's phone who is already here from Islamabad. I want to request American friends and European embassies to please record my conversation with the psychiatrist who will be very bitter and desperate to put me on daily medication because fortnightly injections are not being successful in retarding me.
I want to request good Americans to please stop these racist mind control crooks from stealing my things and continuing my persecution forever.

Amin · May 29th, 2021, 1:52 pm

Friends, I would soon be sharing my C++ program to construct the limit order book from Nasdaq ITCH protocol data. Though I have quickly saved a copy of the C++ program on my USB, I am afraid that crooks will still try to alter my program. I will be running the program several time in next 2-3 days to make sure everything is alright and output of the program is exact and then share it probably on Tuesday after writing explanatory comments on the program. Please stop mind control crooks from altering my programs.

Amin · June 5th, 2021, 4:43 pm

Friends, sorry for being away so long. Lots of things were happening and I had to spend some time away from my research. The update is that I have checked my C++ programs to construct the limit order book from Nasdaq ITCH data and they seem to work perfectly fine. I expect very similar book construction programs to work well for NYSE tick data format as well. I need another day to write detailed comments on LOB construction programs so that friends can easily understand the logic and then I will post the C++ program with comments in another day or two.

Amin · June 7th, 2021, 1:07 pm

Friends, I still could not add comments therefore I have given a brief overview of program logic and I am posting the programs. There are just two files in my program. A header file and a .cpp file.

My algorithm for limit order book construction is very different. There are two deques for "bid limit order book volume" and "ask limit order book volume". These deques start from top of each side of the book and every element of deques stores volume of corresponding limit order book based on its distance(no of levels) from top of the book on that side. These deques only keep track of volume in both bid and ask books. You can limit the size of deque to some large number in order to avoid limit order book being too large due to orders at zero or other other orders much away from the top of the book.
As new limit orders come, they are stored in a vector(of structure with all the characteristics like volume, price, time etc. in the structure) in a serial fashion with length of the "limit order vector" increasing with each new limit order.
When a new limit order arrives in the market, we update the limit order book at the order price by adding appropriate volume and also append the limit order vector with push_back of the structure containing all order information.
The key to algorithm is a map that stores Nasdaq order number when it arrives and maps it to the end of the limit order array where the new order is appended in the limit order vector so whenever we are given NAsdaq Order number, we use the map to look up the appropriate entry in limit order vector where the limit order(with Nasdaq Order number) was stored.
When there is a cancel or delete order, we use the Nasdaq order number and map the position of the Nasdaq order in the limit order vector and limit order vector iterator is pointed to appropriate element of limit order vector and then we retrieve the order price and order size from the limit order vector and use it to decrease the volume at appropriate place in the bid or ask volume deque.
When we have an execute or Fill order with Nasdaq order number, we again look at corresponding element of the limit order vector(using the map) to find out order size and decrease the volume on relevant element of the appropriate limit order volume book. We also check if the Fill or execute order does not lie on top of the book by matching the price information from the appropriate element of limit order vector and in that case algorithm flags an error that book has not been constructed properly.(This error is not flagged in my programs since construction of limit order book remain good.)
After an order has been completely filled or completely deleted, we know that it will not be recalled again and so we delete that entry from the map that points each Nasdaq limit order number to corresponding limit order entry in limit order vector.

#ifndef ConstructLOB02_H
#define ConstructLOB02_H
#pragma once
#include <iostream>
#include <deque>
#include <vector>
#include <list>
#include <fstream>
#include <sstream>
#include<iterator>
#include <string>
#include <iomanip>
#include <cstdint>
#include <math.h>
#include <map>
#include <algorithm>

struct OrderFormat {
	uint32_t TimeStamp;
	std::string OrderSymbol;
	uint64_t OrderNumber;
	char OrderAction;
	uint32_t OrderSize;
	uint32_t OrderPrice;
};

struct LimitOrderRecord {

	uint64_t OrderNumber;
	uint32_t OrderArrivalTime;//Original Time Stamp;
	char OrderType;//B or S;
	char OrderStatus;//U=unfilled,P=partiallyFilled,E=filled,C=Cancelled;
	uint32_t OrderPrice;//Bid or Ask Price;
	uint32_t OrderBookLevel;//BidBookLevel or AskBookLevel;
	uint32_t OrderVolume;//
	uint32_t OrderExecuteOrCancelTime;//Execution Time Stamp;

};

struct LimitOrderLevelDataRecord {


	uint32_t DataTime0;//Original Time Stamp;
	uint32_t BidLevel3;
	uint32_t BidLevel3Volume;
	uint32_t BidLevel2;
	uint32_t BidLevel2Volume;
	uint32_t BidLevel1;
	uint32_t BidLevel1Volume;
	uint32_t AskLevel1;
	uint32_t AskLevel1Volume;
	uint32_t AskLevel2;
	uint32_t AskLevel2Volume;
	uint32_t AskLevel3;
	uint32_t AskLevel3Volume;
};




class LimitOrderBook {


public:
	LimitOrderBook::LimitOrderBook(void);
	void ShowTheOrderFlowRecord(std::string ReasonForJump, uint32_t TimeStamp, uint32_t MBOArrivalRAte,
		uint32_t MSOArrivalRate, int32_t LOBBArrivaRate, int32_t LOBAArivalRate);
	void AskLimitOrderArrivalActionOnLimitOrderBook(uint32_t OrderPrice, uint32_t OrderSize);
	void AskLimitOrderArrivalActionOnLimitOrderVector(uint64_t OrderNumber, uint32_t TimeStamp, uint32_t OrderPrice,
		uint32_t OrderSize);
	void LimitOrderArrivalActionOnOrder_struct(uint64_t OrderNumber);

	void BidLimitOrderArrivalActionOnLimitOrderBook(uint32_t OrderPrice, uint32_t OrderSize);
	void BidLimitOrderArrivalActionOnLimitOrderVector(uint64_t OrderNumber, uint32_t TimeStamp, uint32_t OrderPrice,
		uint32_t OrderSize);
	void BidLimitOrderDeleteOrCancelActionOnLimitOrderBook(uint32_t OrderPrice, uint32_t OrderSize);
	//void BidLimitOrderDeleteActionOnLimitOrderBook(uint64_t OrderNumber);
	//void LimitOrderDeleteActionOnLimitOrderList(uint64_t OrderNumber, uint32_t TimeStamp);
	//void AskLimitOrderDeleteActionOnLimitOrderBook(uint64_t OrderNumber);
	void AskLimitOrderDeleteOrCancelActionOnLimitOrderBook(uint32_t OrderPrice, uint32_t OrderSize);
	//void BidLimitOrderCancelActionOnLimitOrderBook(uint64_t OrderNumber, uint32_t OrderSize);
	//void AskLimitOrderCancelActionOnLimitOrderBook(uint64_t OrderNumber, uint32_t OrderSize);
	//void LimitOrderCancelActionOnLimitOrderList(uint64_t OrderNumber, uint32_t OrderSize, uint32_t TimeStamp);
	void LimitOrderCancelActionOnLimitOrderList(uint32_t DeletedOrderSerialNumber, uint32_t OrderSize, uint32_t TimeStamp);

	void BidLimitOrderFillOrExecuteActionOnLimitOrderBook(uint32_t OrderSize);
	void AskLimitOrderFillOrExecuteActionOnLimitOrderBook(uint32_t OrderSize);
	void LimitOrderExecuteActionOnLimitOrderList(uint32_t ExecutedOrderSerialNumber, uint32_t OrderSize, uint32_t TimeStamp);

	LimitOrderLevelDataRecord LimitOrderLevelData;
	//void ConstructLimitOrderLevelDataRecord(uint32_t DataTime);
	LimitOrderLevelDataRecord ConstructLimitOrderLevelDataRecord(uint32_t DataTime);
	
	std::vector <LimitOrderLevelDataRecord> LimitOrderLevelDataList;
	//LimitOrderLevelDataRecord LimitOrderLevelData;

	uint32_t TopOfAskBook;
	uint32_t TopOfBidBook;

	uint32_t BottomOfAskBook;
	uint32_t BottomOfBidBook;

	std::deque<uint32_t>BidLimitOrderBookVolume;
	std::deque<uint32_t>AskLimitOrderBookVolume;
	std::map<uint64_t, uint32_t> Order_struct;
	uint32_t OrderSerialNumber;
	std::vector <LimitOrderRecord> LimitOrderList;
private:
	std::deque<uint32_t>::iterator Dit1;

};


//Line::Line(void) {
//	cout << "Object is being created" << endl;
//}

LimitOrderBook::LimitOrderBook(void) {
	TopOfAskBook = 0;
	TopOfBidBook = 0;
	BottomOfAskBook = 0;
	BottomOfBidBook = 0;
	OrderSerialNumber = 0;
}


void LimitOrderBook::ShowTheOrderFlowRecord(std::string ReasonForJump, uint32_t TimeStamp, uint32_t MBOArrivalRate,
	uint32_t MSOArrivalRate, int32_t LOBBArrivalRate, int32_t LOBAArrivalRate)
{
	std::cout << "TimeStamp of Order  " << TimeStamp << std::endl;
	std::cout << ReasonForJump << std::endl;
	//std::cout << "OrderFlowRecord.StartTime  " << OrderFlowRecord.StartTime << std::endl;
	//std::cout << "OrderFlowRecord.EndTime  " << OrderFlowRecord.EndTime << std::endl;

	std::cout << "AskLimitOrderExistingVolumeL[nn]   ";
	std::cout << AskLimitOrderBookVolume[0];

	std::cout << std::endl;
	std::cout << "AskLimitNetOrderArrivalVolumeL[nn]   ";

	std::cout << LOBAArrivalRate;
	std::cout << ";";

	std::cout << std::endl;

	std::cout << "MarketBuyOrderArrivalVolume   ";
	std::cout << MBOArrivalRate;
	std::cout << ";";

	std::cout << std::endl;
	std::cout << "TopOfAskBook   ";
	std::cout << TopOfAskBook;
	std::cout << ";";
	std::cout << std::endl;
	std::cout << "TopOfBidBook   ";
	std::cout << TopOfBidBook;
	std::cout << ";";
	std::cout << std::endl;
	std::cout << "MarketSellOrderArrivalVolume   ";
	std::cout << MSOArrivalRate;
	std::cout << ";";
	std::cout << std::endl;
	std::cout << "BuyLimitOrderNetArrivalVolumeL[nn]  ";

	std::cout << LOBBArrivalRate;
	std::cout << ";";
	std::cout << "BuyLimitOrderExistingVolumeL[nn]  ";

	std::cout << BidLimitOrderBookVolume[0];
	std::cout << ";";

	std::cout << std::endl;


	std::cout << std::endl;



	std::cout << std::endl;
}



LimitOrderLevelDataRecord LimitOrderBook::ConstructLimitOrderLevelDataRecord(uint32_t DataTime)
{
	
	LimitOrderLevelData.DataTime0 = DataTime;
	LimitOrderLevelData.BidLevel3 = TopOfBidBook - 2;
	LimitOrderLevelData.BidLevel3Volume = (BidLimitOrderBookVolume.size() >= 3) ? (BidLimitOrderBookVolume[2]) : (0);
	LimitOrderLevelData.BidLevel2 = TopOfBidBook - 1;
	LimitOrderLevelData.BidLevel2Volume = (BidLimitOrderBookVolume.size() >= 2) ? (BidLimitOrderBookVolume[1]) : (0);
	LimitOrderLevelData.BidLevel1 = TopOfBidBook;
	LimitOrderLevelData.BidLevel1Volume = (BidLimitOrderBookVolume.size() >= 1) ? (BidLimitOrderBookVolume[0]) : (0);
	LimitOrderLevelData.AskLevel1 = TopOfAskBook;
	LimitOrderLevelData.AskLevel1Volume = (AskLimitOrderBookVolume.size() >= 1) ? (AskLimitOrderBookVolume[0]) : (0);
	LimitOrderLevelData.AskLevel2 = TopOfAskBook + 1;
	LimitOrderLevelData.AskLevel2Volume = (AskLimitOrderBookVolume.size() >= 2) ? (AskLimitOrderBookVolume[1]) : (0);
	LimitOrderLevelData.AskLevel3 = TopOfAskBook + 2;
	LimitOrderLevelData.AskLevel3Volume = (AskLimitOrderBookVolume.size() >= 3) ? (AskLimitOrderBookVolume[2]) : (0);

	return(LimitOrderLevelData);
}




void LimitOrderBook::BidLimitOrderArrivalActionOnLimitOrderBook(uint32_t OrderPrice, uint32_t OrderSize)
{
	//OrderPrice = (uint32_t)std::round(OrderPrice / 100);


	if (BidLimitOrderBookVolume.size() > 0)
	{
		if ((OrderPrice <= TopOfBidBook) && (OrderPrice >= TopOfBidBook - BidLimitOrderBookVolume.size() + 1))
		{
			uint32_t index = TopOfBidBook - OrderPrice;
			Dit1 = BidLimitOrderBookVolume.begin() + index;
			//*Dit1 = *Dit1 + it1->OrderSize;
			*Dit1 = *Dit1 + OrderSize;


		}
		else if (OrderPrice > TopOfBidBook)
		{

			uint32_t IndicesToAdd = OrderPrice - TopOfBidBook;
			TopOfBidBook = OrderPrice;
			for (uint32_t ii = 0; ii < IndicesToAdd; ii++)
			{
				BidLimitOrderBookVolume.push_front(0);
			}
			Dit1 = BidLimitOrderBookVolume.begin();
			*Dit1 = OrderSize;


		}

		else if (OrderPrice < TopOfBidBook - BidLimitOrderBookVolume.size() + 1)//check to make sure that ends are aligned.
		{
			BottomOfBidBook = OrderPrice;
			uint32_t IndicesToAdd = TopOfBidBook - BidLimitOrderBookVolume.size() + 1 - OrderPrice;
			if (IndicesToAdd < 50000)
			{
				Dit1 = BidLimitOrderBookVolume.end();
				BidLimitOrderBookVolume.insert(Dit1, IndicesToAdd, 0);
				Dit1 = BidLimitOrderBookVolume.end() - 1;
				*Dit1 = OrderSize;

			}
		}

	//	mydeque.insert(it, 2, 20);
    //
	//	else if (OrderPrice < TopOfBidBook - BidLimitOrderBookVolume.size() + 1)//check to make sure that ends are aligned.
	//	{
	//		uint32_t IndicesToAdd = TopOfBidBook - BidLimitOrderBookVolume.size() + 1 - OrderPrice;
//
	//		for (uint32_t ii = 0; ii < IndicesToAdd; ii++)
	//		{
	//			BidLimitOrderBookVolume.push_back(0);
	//		}
	//		Dit1 = BidLimitOrderBookVolume.end() - 1;
	//		*Dit1 = OrderSize;
	//	}
	}
	else
	{
		TopOfBidBook = OrderPrice;
		BottomOfBidBook = OrderPrice;
		//Dit1 = BidLimitOrderBookVolume.begin();
		//BidLimitOrderBookVolume.insert(Dit1,it1->OrderSize);
		BidLimitOrderBookVolume.push_back(OrderSize);

	}
}

void LimitOrderBook::BidLimitOrderArrivalActionOnLimitOrderVector(uint64_t OrderNumber, uint32_t TimeStamp, uint32_t OrderPrice,
	uint32_t OrderSize)
{
	LimitOrderRecord NewLimitOrder;
	NewLimitOrder.OrderNumber = OrderNumber;
	NewLimitOrder.OrderArrivalTime = TimeStamp;
	NewLimitOrder.OrderPrice = OrderPrice;
	NewLimitOrder.OrderStatus = 'U';
	NewLimitOrder.OrderType = 'B';
	NewLimitOrder.OrderVolume = OrderSize;
	NewLimitOrder.OrderBookLevel = TopOfBidBook - (NewLimitOrder.OrderPrice);
	NewLimitOrder.OrderExecuteOrCancelTime = 0;

	LimitOrderList.push_back(NewLimitOrder);

}

void LimitOrderBook::LimitOrderArrivalActionOnOrder_struct(uint64_t OrderNumber)
{
	//Order_struct.insert({ it1->OrderNumber, OrderSerialNumber });
	Order_struct.insert({ OrderNumber, OrderSerialNumber });
	OrderSerialNumber++;
}


void LimitOrderBook::AskLimitOrderArrivalActionOnLimitOrderBook(uint32_t OrderPrice, uint32_t OrderSize)
{
	//OrderPrice = (uint32_t)std::round(OrderPrice / 100);

	if (AskLimitOrderBookVolume.size() > 0)
	{
		if ((OrderPrice >= TopOfAskBook) && (OrderPrice <= TopOfAskBook + AskLimitOrderBookVolume.size() - 1))
		{
			uint32_t index = OrderPrice - TopOfAskBook;
			Dit1 = AskLimitOrderBookVolume.begin() + index;
			*Dit1 = *Dit1 + OrderSize;
		}
		else if (OrderPrice < TopOfAskBook)
		{
			uint32_t IndicesToAdd = TopOfAskBook - OrderPrice;
			TopOfAskBook = OrderPrice;
			for (uint32_t ii = 0; ii < IndicesToAdd; ii++)
			{
				AskLimitOrderBookVolume.push_front(0);
			}
			Dit1 = AskLimitOrderBookVolume.begin();
			*Dit1 = OrderSize;

		}
		//else if (OrderPrice > TopOfAskBook + AskLimitOrderBookVolume.size() - 1)//check to make sure that ends are aligned.
		//{
		//	uint32_t IndicesToAdd = OrderPrice - (TopOfAskBook + AskLimitOrderBookVolume.size() - 1);
		//	for (uint32_t ii = 0; ii < IndicesToAdd; ii++)
		//	{
		//		AskLimitOrderBookVolume.push_back(0);
		//	}
		//	Dit1 = AskLimitOrderBookVolume.end() - 1;
		//	*Dit1 = OrderSize;
		//
		//}
		else if (OrderPrice > TopOfAskBook + AskLimitOrderBookVolume.size() - 1)//check to make sure that ends are aligned.
		{
			BottomOfAskBook = OrderPrice;
			uint32_t IndicesToAdd = OrderPrice - (TopOfAskBook + AskLimitOrderBookVolume.size() - 1);
			if (IndicesToAdd < 50000)
			{
				Dit1 = AskLimitOrderBookVolume.end();
				AskLimitOrderBookVolume.insert(Dit1, IndicesToAdd, 0);
				Dit1 = AskLimitOrderBookVolume.end() - 1;
				*Dit1 = OrderSize;
			}
		}
	}
	else
	{

		TopOfAskBook = OrderPrice;
		BottomOfAskBook = OrderPrice;
		AskLimitOrderBookVolume.push_back(OrderSize);


	}
}
void LimitOrderBook::AskLimitOrderArrivalActionOnLimitOrderVector(uint64_t OrderNumber, uint32_t TimeStamp, uint32_t OrderPrice,
	uint32_t OrderSize)
{
	LimitOrderRecord NewLimitOrder;

	NewLimitOrder.OrderNumber = OrderNumber;
	NewLimitOrder.OrderArrivalTime = TimeStamp;
	NewLimitOrder.OrderPrice = OrderPrice;

	//NewLimitOrder.OrderPrice = it1->OrderPrice; //has already been assigned value
	NewLimitOrder.OrderStatus = 'U';
	NewLimitOrder.OrderType = 'S';
	NewLimitOrder.OrderVolume = OrderSize;
	NewLimitOrder.OrderBookLevel = (NewLimitOrder.OrderPrice) - TopOfAskBook;
	NewLimitOrder.OrderExecuteOrCancelTime = 0;

	LimitOrderList.push_back(NewLimitOrder);
}



//void LimitOrderBook::BidLimitOrderDeleteActionOnLimitOrderBook(uint64_t OrderNumber)
//{
//	std::map<uint64_t, uint32_t>::iterator findIter = Order_struct.find(OrderNumber);
//	uint32_t DeletedOrderSerialNumber = findIter->second;
//	Order_struct.erase(findIter);
//	
//	uint32_t OrderPrice = LimitOrderList[DeletedOrderSerialNumber].OrderPrice;
//	uint32_t OrderSize = LimitOrderList[DeletedOrderSerialNumber].OrderVolume;
void LimitOrderBook::BidLimitOrderDeleteOrCancelActionOnLimitOrderBook(uint32_t OrderPrice, uint32_t OrderSize)
{
	if ((OrderPrice < TopOfBidBook) && (OrderPrice > TopOfBidBook - BidLimitOrderBookVolume.size() + 1))
	{
		uint32_t index = TopOfBidBook - OrderPrice;
		Dit1 = BidLimitOrderBookVolume.begin() + index;
		//	std::cout << "existing volume  " << *Dit1 << std::endl;
		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);

	}
	else if (OrderPrice == TopOfBidBook)
	{

		Dit1 = BidLimitOrderBookVolume.begin();
		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
		if (*Dit1 == 0)
		{
			bool DeleteFlag = 1;
			while (DeleteFlag)
			{
				BidLimitOrderBookVolume.pop_front();
				TopOfBidBook = TopOfBidBook - 1;
				if (!BidLimitOrderBookVolume.empty())
				{
					Dit1 = BidLimitOrderBookVolume.begin();
					if (*Dit1 == 0)
					{
						DeleteFlag = 1;

					}
					else
					{
						DeleteFlag = 0;
					}
				}
				else
				{
					DeleteFlag = 0;
					TopOfBidBook = 0;
					BottomOfBidBook = 0;
				}

			}
		}

	}
	else if (OrderPrice == TopOfBidBook - BidLimitOrderBookVolume.size() + 1)
	{
		uint32_t index = BidLimitOrderBookVolume.size() - 1;

		Dit1 = BidLimitOrderBookVolume.end() - 1;
		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
		//*Dit1 = *Dit1 - LimitOrderDeleted.OrderVolume;
		if (*Dit1 == 0)
		{
			bool DeleteFlag = 1;
			while (DeleteFlag)
			{
				BidLimitOrderBookVolume.pop_back();
				BottomOfBidBook = BottomOfBidBook + 1;
				if (!BidLimitOrderBookVolume.empty())
				{
					Dit1 = BidLimitOrderBookVolume.end() - 1;
					if (*Dit1 == 0)
					{
						DeleteFlag = 1;

					}
					else
					{
						DeleteFlag = 0;
					}
				}
				else
				{
					DeleteFlag = 0;
					TopOfBidBook = 0;
					BottomOfBidBook = 0;
				}

			}
		}

	}
}

//void LimitOrderBook::LimitOrderDeletionActionOnLimitOrderList(uint64_t OrderNumber, uint32_t TimeStamp)
//{
//	std::map<uint64_t, uint32_t>::iterator findIter = Order_struct.find(OrderNumber);
//	uint32_t DeletedOrderSerialNumber = findIter->second;
//	Order_struct.erase(findIter);
//
//	LimitOrderList[DeletedOrderSerialNumber].OrderExecuteOrCancelTime = TimeStamp;
//	LimitOrderList[DeletedOrderSerialNumber].OrderStatus = 'D';
//	//For OrderAction D, volume of deleted Limit Order is not changed.
//}


//void LimitOrderBook::AskLimitOrderDeleteActionOnLimitOrderBook(uint64_t OrderNumber)
//{
//	std::map<uint64_t, uint32_t>::iterator findIter = Order_struct.find(OrderNumber);
//	uint32_t DeletedOrderSerialNumber = findIter->second;
//	Order_struct.erase(findIter);
//
//	uint32_t OrderPrice = LimitOrderList[DeletedOrderSerialNumber].OrderPrice;
//	uint32_t OrderSize = LimitOrderList[DeletedOrderSerialNumber].OrderVolume;

void LimitOrderBook::AskLimitOrderDeleteOrCancelActionOnLimitOrderBook(uint32_t OrderPrice, uint32_t OrderSize)
{
	if ((OrderPrice > TopOfAskBook) && (OrderPrice < TopOfAskBook + AskLimitOrderBookVolume.size() - 1))
	{
		uint32_t index = OrderPrice - TopOfAskBook;
		Dit1 = AskLimitOrderBookVolume.begin() + index;
		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
		//*Dit1 = *Dit1 - LimitOrderDeleted.OrderVolume;
	}
	else if (OrderPrice == TopOfAskBook)
	{

		Dit1 = AskLimitOrderBookVolume.begin();
		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
		if (*Dit1 == 0)
		{
			bool DeleteFlag = 1;
			while (DeleteFlag)
			{
				AskLimitOrderBookVolume.pop_front();
				TopOfAskBook = TopOfAskBook + 1;
				if (!AskLimitOrderBookVolume.empty())
				{
					Dit1 = AskLimitOrderBookVolume.begin();
					if (*Dit1 == 0)
					{
						DeleteFlag = 1;

					}
					else
					{
						DeleteFlag = 0;
					}
				}
				else
				{
					DeleteFlag = 0;
					TopOfAskBook = 0;
					BottomOfAskBook = 0;
				}

			}
		}
	}
	else if (OrderPrice == TopOfAskBook + AskLimitOrderBookVolume.size() - 1)
	{
		uint32_t index = AskLimitOrderBookVolume.size() - 1;
		Dit1 = AskLimitOrderBookVolume.end() - 1;
		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
		//*Dit1 = *Dit1 - LimitOrderDeleted.OrderVolume;
		if (*Dit1 == 0)
		{
			bool DeleteFlag = 1;
			while (DeleteFlag)
			{
				AskLimitOrderBookVolume.pop_back();
				BottomOfAskBook = BottomOfAskBook - 1;
				if (!AskLimitOrderBookVolume.empty())
				{
					Dit1 = AskLimitOrderBookVolume.end() - 1;
					if (*Dit1 == 0)
					{
						DeleteFlag = 1;
					}
					else
					{
						DeleteFlag = 0;
					}
				}
				else
				{
					DeleteFlag = 0;
					TopOfAskBook = 0;
					BottomOfAskBook = 0;
				}
			}
		}
	}
}

//void LimitOrderBook::BidLimitOrderCancelActionOnLimitOrderBook(uint64_t OrderPrice, uint32_t OrderSize)
//{
////	std::map<uint64_t, uint32_t>::iterator findIter = Order_struct.find(OrderNumber);
////	uint32_t DeletedOrderSerialNumber = findIter->second;
////	//Order_struct.erase(findIter);//must not be erased in partial cancel/delete.
//
//	uint32_t OrderPrice = LimitOrderList[DeletedOrderSerialNumber].OrderPrice;
//
//	if ((OrderPrice < TopOfBidBook) && (OrderPrice > TopOfBidBook - BidLimitOrderBookVolume.size() + 1))
//	{
//		uint32_t index = TopOfBidBook - OrderPrice;
//		Dit1 = BidLimitOrderBookVolume.begin() + index;
//		//	std::cout << "existing volume  " << *Dit1 << std::endl;
//		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
//
//	}
//	else if (OrderPrice == TopOfBidBook)
//	{
//
//		Dit1 = BidLimitOrderBookVolume.begin();
//		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
//		if (*Dit1 == 0)
//		{
//			bool DeleteFlag = 1;
//			while (DeleteFlag)
//			{
//				BidLimitOrderBookVolume.pop_front();
//				TopOfBidBook = TopOfBidBook - 1;
//				if (!BidLimitOrderBookVolume.empty())
//				{
//					Dit1 = BidLimitOrderBookVolume.begin();
//					if (*Dit1 == 0)
//					{
//						DeleteFlag = 1;
//
//					}
//					else
//					{
//						DeleteFlag = 0;
//					}
//				}
//				else
//				{
//					DeleteFlag = 0;
//					TopOfBidBook = 0;
//				}
//
//			}
//		}
//
//	}
//	else if (OrderPrice == TopOfBidBook - BidLimitOrderBookVolume.size() + 1)
//	{
//		uint32_t index = BidLimitOrderBookVolume.size() - 1;
//
//		Dit1 = BidLimitOrderBookVolume.end() - 1;
//		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
//		//*Dit1 = *Dit1 - LimitOrderDeleted.OrderVolume;
//		if (*Dit1 == 0)
//		{
//			bool DeleteFlag = 1;
//			while (DeleteFlag)
//			{
//				BidLimitOrderBookVolume.pop_back();
//				if (!BidLimitOrderBookVolume.empty())
//				{
//					Dit1 = BidLimitOrderBookVolume.end() - 1;
//					if (*Dit1 == 0)
//					{
//						DeleteFlag = 1;
//
//					}
//					else
//					{
//						DeleteFlag = 0;
//					}
//				}
//				else
//				{
//					DeleteFlag = 0;
//					TopOfBidBook = 0;
//				}
//
//			}
//		}
//
//	}
//}
//
////void LimitOrderBook::AskLimitOrderCancelActionOnLimitOrderBook(uint64_t OrderNumber, uint32_t OrderSize);
//
//void LimitOrderBook::AskLimitOrderCancelActionOnLimitOrderBook(uint64_t OrderNumber, uint32_t OrderSize)
//{
//	std::map<uint64_t, uint32_t>::iterator findIter = Order_struct.find(OrderNumber);
//	uint32_t DeletedOrderSerialNumber = findIter->second;
//	//Order_struct.erase(findIter);//must not be erased in partial cancel/delete.
//
//	uint32_t OrderPrice = LimitOrderList[DeletedOrderSerialNumber].OrderPrice;
//	//uint32_t OrderSize = LimitOrderList[DeletedOrderSerialNumber].OrderVolume;
//
//	if ((OrderPrice > TopOfAskBook) && (OrderPrice < TopOfAskBook + AskLimitOrderBookVolume.size() - 1))
//	{
//		uint32_t index = OrderPrice - TopOfAskBook;
//		Dit1 = AskLimitOrderBookVolume.begin() + index;
//		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
//		//*Dit1 = *Dit1 - LimitOrderDeleted.OrderVolume;
//	}
//	else if (OrderPrice == TopOfAskBook)
//	{
//
//		Dit1 = AskLimitOrderBookVolume.begin();
//		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
//		if (*Dit1 == 0)
//		{
//			bool DeleteFlag = 1;
//			while (DeleteFlag)
//			{
//				AskLimitOrderBookVolume.pop_front();
//				TopOfAskBook = TopOfAskBook + 1;
//				if (!AskLimitOrderBookVolume.empty())
//				{
//					Dit1 = AskLimitOrderBookVolume.begin();
//					if (*Dit1 == 0)
//					{
//						DeleteFlag = 1;
//
//					}
//					else
//					{
//						DeleteFlag = 0;
//					}
//				}
//				else
//				{
//					DeleteFlag = 0;
//					TopOfAskBook = 0;
//				}
//
//			}
//		}
//	}
//	else if (OrderPrice == TopOfAskBook + AskLimitOrderBookVolume.size() - 1)
//	{
//		uint32_t index = AskLimitOrderBookVolume.size() - 1;
//		Dit1 = AskLimitOrderBookVolume.end() - 1;
//		*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
//		//*Dit1 = *Dit1 - LimitOrderDeleted.OrderVolume;
//		if (*Dit1 == 0)
//		{
//			bool DeleteFlag = 1;
//			while (DeleteFlag)
//			{
//				AskLimitOrderBookVolume.pop_back();
//				if (!AskLimitOrderBookVolume.empty())
//				{
//					Dit1 = AskLimitOrderBookVolume.end() - 1;
//					if (*Dit1 == 0)
//					{
//						DeleteFlag = 1;
//					}
//					else
//					{
//						DeleteFlag = 0;
//					}
//				}
//				else
//				{
//					DeleteFlag = 0;
//					TopOfAskBook = 0;
//				}
//			}
//		}
//	}
//}

void LimitOrderBook::LimitOrderCancelActionOnLimitOrderList(uint32_t DeletedOrderSerialNumber, uint32_t OrderSize, uint32_t TimeStamp)
{
	//std::map<uint64_t, uint32_t>::iterator findIter = Order_struct.find(OrderNumber);
	//uint32_t DeletedOrderSerialNumber = findIter->second;
	////Order_struct.erase(findIter);

	LimitOrderList[DeletedOrderSerialNumber].OrderExecuteOrCancelTime = TimeStamp;
	LimitOrderList[DeletedOrderSerialNumber].OrderStatus = 'C';
	if (OrderSize <= LimitOrderList[DeletedOrderSerialNumber].OrderVolume)
	{
		LimitOrderList[DeletedOrderSerialNumber].OrderVolume =
			LimitOrderList[DeletedOrderSerialNumber].OrderVolume - OrderSize;
	}
	else
	{
		LimitOrderList[DeletedOrderSerialNumber].OrderVolume = 0;
	}
}


//
//
////FilledOrderIn.OrderNumber = it1->OrderNumber;
////FilledOrderIn.OrderSerialNumber = 0;
//FilledOrderIn = it1->OrderNumber;
//
////goodliffe::skip_list<OrderStruct>::iterator findIter = Order_struct.find(FilledOrderIn);
//std::map<uint64_t, uint32_t>::iterator findIter = Order_struct.find(FilledOrderIn);
//FilledOrderOut = findIter->second;
//Order_struct.erase(findIter);
//
//LimitOrderFilled = LimitOrderList[FilledOrderOut];
//
//void LimitOrderBook::BidLimitOrderFillActionOnLimitOrderBook(uint64_t OrderNumber, uint32_t OrderSize);
//{
//	//NewMarketOrder.OrderSerialNumber = MarketOrderSerialNumber;
//	//NewMarketOrder.OrderArrivalTime = it1->TimeStamp;//Original Time Stamp;
//	//NewMarketOrder.OrderType = 'S';//B or S;
//	//NewMarketOrder.OrderPrice1stLevel = LimitOrderFilled.OrderPrice;//Bid or Ask Price;
//	MarketOrderList[MarketOrderSerialNumber - 1].OrderPriceLastLevel = LimitOrderFilled.OrderPrice;//Bid or Ask Price;
//	MarketOrderList[MarketOrderSerialNumber - 1].LimitOrdersFilledCount++;
//	MarketOrderList[MarketOrderSerialNumber - 1].OrderVolume =
//		MarketOrderList[MarketOrderSerialNumber - 1].OrderVolume + LimitOrderFilled.OrderVolume;//
//


//if ((LimitOrderFilled.OrderType == 'B') && (BPrevTimeStamp != it1->TimeStamp))
//{
//	NewMarketOrder.OrderSerialNumber = MarketOrderSerialNumber;
//	NewMarketOrder.OrderArrivalTime = it1->TimeStamp;//Original Time Stamp;
//	NewMarketOrder.OrderType = 'S';//B or S;
//	NewMarketOrder.OrderPrice1stLevel = LimitOrderFilled.OrderPrice;//Bid or Ask Price;
//	NewMarketOrder.OrderPriceLastLevel = LimitOrderFilled.OrderPrice;//Bid or Ask Price;
//	NewMarketOrder.LimitOrdersFilledCount = 1;
//	NewMarketOrder.OrderVolume = LimitOrderFilled.OrderVolume;//
//
//
//
//	MarketOrderList.push_back(NewMarketOrder);
//	MarketOrderSerialNumber++;
//	BPrevTimeStamp = it1->TimeStamp;
//	}
//}



void LimitOrderBook::BidLimitOrderFillOrExecuteActionOnLimitOrderBook(uint32_t OrderSize)
{

	Dit1 = BidLimitOrderBookVolume.begin();
	*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
	//*Dit1 = *Dit1 - LimitOrderFilled.OrderVolume;
	if (*Dit1 == 0)
	{
		bool DeleteFlag = 1;
		while (DeleteFlag)
		{
			BidLimitOrderBookVolume.pop_front();
			TopOfBidBook = TopOfBidBook - 1;

			if (!BidLimitOrderBookVolume.empty())
			{
				Dit1 = BidLimitOrderBookVolume.begin();
				if (*Dit1 == 0)
				{
					DeleteFlag = 1;

				}
				else
				{
					DeleteFlag = 0;
				}
			}
			else
			{
				DeleteFlag = 0;
				TopOfBidBook = 0;
				BottomOfBidBook = 0;
			}

		}
	}

}


void LimitOrderBook::AskLimitOrderFillOrExecuteActionOnLimitOrderBook(uint32_t OrderSize)
{

	Dit1 = AskLimitOrderBookVolume.begin();
	*Dit1 = *Dit1 - ((*Dit1 >= OrderSize) ? OrderSize : *Dit1);
	//*Dit1 = *Dit1 - LimitOrderFilled.OrderVolume;
	if (*Dit1 == 0)
	{
		bool DeleteFlag = 1;
		while (DeleteFlag)
		{
			AskLimitOrderBookVolume.pop_front();
			TopOfAskBook = TopOfAskBook + 1;
			if (!AskLimitOrderBookVolume.empty())
			{
				Dit1 = AskLimitOrderBookVolume.begin();
				if (*Dit1 == 0)
				{
					DeleteFlag = 1;
				}
				else
				{
					DeleteFlag = 0;
				}
			}
			else
			{
				DeleteFlag = 0;
				TopOfAskBook = 0;
				BottomOfAskBook = 0;
			}
		}

	}
}
void LimitOrderBook::LimitOrderExecuteActionOnLimitOrderList(uint32_t ExecutedOrderSerialNumber, uint32_t OrderSize, uint32_t TimeStamp)
{
	//std::map<uint64_t, uint32_t>::iterator findIter = Order_struct.find(OrderNumber);
	//uint32_t DeletedOrderSerialNumber = findIter->second;
	////Order_struct.erase(findIter);

	LimitOrderList[ExecutedOrderSerialNumber].OrderExecuteOrCancelTime = TimeStamp;
	LimitOrderList[ExecutedOrderSerialNumber].OrderStatus = 'E';
	if (OrderSize <= LimitOrderList[ExecutedOrderSerialNumber].OrderVolume)
	{
		LimitOrderList[ExecutedOrderSerialNumber].OrderVolume =
			LimitOrderList[ExecutedOrderSerialNumber].OrderVolume - OrderSize;
	}
	else
	{
		LimitOrderList[ExecutedOrderSerialNumber].OrderVolume = 0;
	}
}





#endif

.
.Here is the cpp file.

#include "stdafx.h"
#include "ConstructLOB02.h"


int main()
{

	//const char* filename = "C:\\Users\\Ahsan\\Documents\\Visual Studio 2015\\Projects\\ConstructLOB02\\20150109_EBAY.txt";
	const char* filename = "C:\\Users\\Ahsan\\Documents\\Visual Studio 2015\\Projects\\ConstructLOB02\\20150106_BMRN.txt";
	std::ifstream ItchFile(filename);
	OrderFormat OrderNew;
	if (!ItchFile.is_open())
	{
		std::cerr << "Error opening file";
		exit(EXIT_FAILURE);
	}
	std::vector<OrderFormat> OrderT;// (1, OrderNew);

	std::string OrderMPID;
	char X;

	uint32_t tnn = 0;

	std::string str;
	while (getline(ItchFile, str)) { // Read the entire line
		std::istringstream ss(str);
		ss >> OrderNew.TimeStamp >> OrderNew.OrderSymbol >> OrderNew.OrderNumber >> OrderNew.OrderAction >> OrderNew.OrderSize >> OrderNew.OrderPrice;
		OrderT.push_back(OrderNew);
	}



	//	while (!ItchFile.eof())
	//	{
	//		ItchFile >> OrderNew.TimeStamp;
	//		ItchFile >> OrderNew.OrderSymbol;
	//		ItchFile >> OrderNew.OrderNumber;
	//		ItchFile >> OrderNew.OrderAction;
	//		ItchFile >> OrderNew.OrderSize;
	//		ItchFile >> OrderNew.OrderPrice;
	//		ItchFile >> OrderMPID;
	//		ItchFile >> X;
	//		OrderT.push_back(OrderNew);
	//	}

	LimitOrderBook LimitOrderBookLive;

	uint32_t OrderPrice0;
	uint32_t OrderSize0;
	uint64_t OrderNumber0;
	uint32_t TimeStamp0;
	//uint32_t OrderSerialNumber = 0;



	uint32_t TimeGranularity_Deltat = 1000;
	uint32_t MarketStart_T = 34200000;
	uint32_t MarketEnd_T = 57600000;
	uint32_t CurrentPeriodStart_T = MarketStart_T;
	uint32_t CurrentPeriodEnd_T = CurrentPeriodStart_T + TimeGranularity_Deltat;
	uint32_t TimeStamp_T;

	




	uint32_t TopOfBidBookBefore = 0;
	uint32_t TopOfBidBookAfter = 0;
	uint32_t TopOfAskBookBefore = 0;
	uint32_t TopOfAskBookAfter = 0;


	std::string ReasonForJump;




	std::ofstream BidAndAskLevelSeries;
	BidAndAskLevelSeries.open("C:\\Users\\Ahsan\\Documents\\Visual Studio 2015\\Projects\\ConstructLOB02\\BidAndAskSeries_20150106_BMRN.csv");

	uint32_t Time_t2 = MarketStart_T;
	//uint32_t Time_t2;
	//Time_t2 = Time_t1 + TimeGranularity_Deltat;



	for (std::vector<OrderFormat>::const_iterator it1 = OrderT.cbegin(); it1 < OrderT.cend() - 1; it1++)//would this -1 cause a problem if the iterator has a null value;
	{
		TimeStamp_T = it1->TimeStamp;

		LimitOrderLevelDataRecord LimitOrderBookLevelData0;
		uint32_t tt;
		if ((TimeStamp_T > Time_t2) && (TimeStamp_T<=MarketEnd_T))
		{
			while (Time_t2< TimeStamp_T)
			{
				tt = Time_t2;

				LimitOrderBookLevelData0 = LimitOrderBookLive.ConstructLimitOrderLevelDataRecord(tt);
				LimitOrderBookLive.LimitOrderLevelDataList.push_back(LimitOrderBookLevelData0);
				BidAndAskLevelSeries << LimitOrderBookLevelData0.DataTime0 << "," << LimitOrderBookLevelData0.BidLevel3 << "," << LimitOrderBookLevelData0.BidLevel3Volume << "," << LimitOrderBookLevelData0.BidLevel2 << "," << LimitOrderBookLevelData0.BidLevel2Volume << "," << LimitOrderBookLevelData0.BidLevel1 << "," << LimitOrderBookLevelData0.BidLevel1Volume << "," << LimitOrderBookLevelData0.AskLevel1 << "," << LimitOrderBookLevelData0.AskLevel1Volume << "," << LimitOrderBookLevelData0.AskLevel2 << "," << LimitOrderBookLevelData0.AskLevel2Volume << "," << LimitOrderBookLevelData0.AskLevel3 << "," << LimitOrderBookLevelData0.AskLevel3Volume << std::endl;

				Time_t2 = Time_t2 + TimeGranularity_Deltat;
			}
		}


		//std::cout << "Size of Limit Order Array " << LimitOrderBookLive.LimitOrderList.size()-1 << " Order Serial Number " << LimitOrderBookLive.OrderSerialNumber;




		std::cout << std::setw(10) << it1->OrderNumber << std::endl;
		std::cout << std::setw(10) << it1->TimeStamp << std::endl;
		OrderPrice0 = (uint32_t)std::round(it1->OrderPrice / 100);
		OrderSize0 = it1->OrderSize;
		OrderNumber0 = it1->OrderNumber;
		TimeStamp0 = it1->TimeStamp;
		if (it1->OrderAction == 'B')
		{

			TopOfBidBookBefore = LimitOrderBookLive.TopOfBidBook;
			LimitOrderBookLive.BidLimitOrderArrivalActionOnLimitOrderBook(OrderPrice0, OrderSize0);
			LimitOrderBookLive.BidLimitOrderArrivalActionOnLimitOrderVector(OrderNumber0, TimeStamp0, OrderPrice0, OrderSize0);
			LimitOrderBookLive.LimitOrderArrivalActionOnOrder_struct(OrderNumber0);
			TopOfBidBookAfter = LimitOrderBookLive.TopOfBidBook;
		}

		if (it1->OrderAction == 'S')
		{

			TopOfAskBookBefore = LimitOrderBookLive.TopOfAskBook;

			LimitOrderBookLive.AskLimitOrderArrivalActionOnLimitOrderBook(OrderPrice0, OrderSize0);
			LimitOrderBookLive.AskLimitOrderArrivalActionOnLimitOrderVector(OrderNumber0, TimeStamp0, OrderPrice0, OrderSize0);
			LimitOrderBookLive.LimitOrderArrivalActionOnOrder_struct(OrderNumber0);
			TopOfAskBookAfter = LimitOrderBookLive.TopOfAskBook;
		}
		if (it1->OrderAction == 'D')
		{


			std::map<uint64_t, uint32_t>::iterator findIter = LimitOrderBookLive.Order_struct.find(OrderNumber0);
			uint32_t DeletedOrderSerialNumber = findIter->second;
			LimitOrderBookLive.Order_struct.erase(findIter);

			uint32_t OrderPrice1 = LimitOrderBookLive.LimitOrderList[DeletedOrderSerialNumber].OrderPrice;
			uint32_t OrderSize1 = LimitOrderBookLive.LimitOrderList[DeletedOrderSerialNumber].OrderVolume;
			char OrderType1 = LimitOrderBookLive.LimitOrderList[DeletedOrderSerialNumber].OrderType;

			if (OrderType1 == 'B')
			{

				TopOfBidBookBefore = LimitOrderBookLive.TopOfBidBook;

				
				LimitOrderBookLive.BidLimitOrderDeleteOrCancelActionOnLimitOrderBook(OrderPrice1, OrderSize1);
				TopOfBidBookAfter = LimitOrderBookLive.TopOfBidBook;

			}
			if (OrderType1 == 'S')
			{
				TopOfAskBookBefore = LimitOrderBookLive.TopOfAskBook;

				LimitOrderBookLive.AskLimitOrderDeleteOrCancelActionOnLimitOrderBook(OrderPrice1, OrderSize1);
				TopOfAskBookAfter = LimitOrderBookLive.TopOfAskBook;
			}
		}
		if (it1->OrderAction == 'C')
		{


			std::map<uint64_t, uint32_t>::iterator findIter = LimitOrderBookLive.Order_struct.find(OrderNumber0);
			uint32_t DeletedOrderSerialNumber = findIter->second;
			//LimitOrderBookLive.Order_struct.erase(findIter);

			uint32_t OrderPrice1 = LimitOrderBookLive.LimitOrderList[DeletedOrderSerialNumber].OrderPrice;
			uint32_t OrderSize1 = OrderSize0;
			char OrderType1 = LimitOrderBookLive.LimitOrderList[DeletedOrderSerialNumber].OrderType;

			if (OrderType1 == 'B')
			{

				TopOfBidBookBefore = LimitOrderBookLive.TopOfBidBook;
				LimitOrderBookLive.BidLimitOrderDeleteOrCancelActionOnLimitOrderBook(OrderPrice1, OrderSize1);
				LimitOrderBookLive.LimitOrderCancelActionOnLimitOrderList(DeletedOrderSerialNumber, OrderSize1, TimeStamp0);
				TopOfBidBookAfter = LimitOrderBookLive.TopOfBidBook;
			}
			if (OrderType1 == 'S')
			{
				TopOfAskBookBefore = LimitOrderBookLive.TopOfAskBook;
				LimitOrderBookLive.AskLimitOrderDeleteOrCancelActionOnLimitOrderBook(OrderPrice1, OrderSize1);
				LimitOrderBookLive.LimitOrderCancelActionOnLimitOrderList(DeletedOrderSerialNumber, OrderSize1, TimeStamp0);
				TopOfAskBookAfter = LimitOrderBookLive.TopOfAskBook;
			}
		}
		if (it1->OrderAction == 'F')
		{
			std::map<uint64_t, uint32_t>::iterator findIter = LimitOrderBookLive.Order_struct.find(OrderNumber0);
			uint32_t LimitOrderSerialNumber = findIter->second;
			LimitOrderBookLive.Order_struct.erase(findIter);

			uint32_t OrderSize1 = LimitOrderBookLive.LimitOrderList[LimitOrderSerialNumber].OrderVolume;
			char OrderType1 = LimitOrderBookLive.LimitOrderList[LimitOrderSerialNumber].OrderType;
			//LimitOrderBookLive.LimitOrderList[DeletedOrderSerialNumber].OrderPrice;
			uint32_t OrderPrice1 = LimitOrderBookLive.LimitOrderList[LimitOrderSerialNumber].OrderPrice;


			if (OrderType1 == 'B')
			{
				TopOfBidBookBefore = LimitOrderBookLive.TopOfBidBook;
				if (OrderPrice1 != TopOfBidBookBefore)
				{
					std::cout << " There is an error. Market Order Does not match with Price; Market Fill " << std::endl;
					std::cout << "  Limit Order Price is " << OrderPrice1 << std::endl;
					std::cout << "  Top Of Bid Book is " << TopOfBidBookBefore << std::endl;
				}
				LimitOrderBookLive.BidLimitOrderFillOrExecuteActionOnLimitOrderBook(OrderSize1);
				TopOfBidBookAfter = LimitOrderBookLive.TopOfBidBook;
			}
			if (OrderType1 == 'S')
			{
				TopOfAskBookBefore = LimitOrderBookLive.TopOfAskBook;
				if (OrderPrice1 != TopOfAskBookBefore)
				{
					std::cout << " There is an error. Market Order Does not match with Price; Market Fill " << std::endl;
					std::cout << "  Limit Order Price is " << OrderPrice1 << std::endl;
					std::cout << "  Top Of Ask Book is " << TopOfAskBookBefore << std::endl;
				}
				LimitOrderBookLive.AskLimitOrderFillOrExecuteActionOnLimitOrderBook(OrderSize1);
				TopOfAskBookAfter = LimitOrderBookLive.TopOfAskBook;
			}
		}
		if (it1->OrderAction == 'E')
		{
			//		std::cout << "I am in Part E" << "Top OF Bid Book   " << LimitOrderBookLive.TopOfBidBook << "Top OF Ask Book   " << LimitOrderBookLive.TopOfAskBook << std::endl;
			//FilledOrderIn.OrderNumber = it1->OrderNumber;
			//FilledOrderIn.OrderSerialNumber = 0;
			//uint64_t FilledOrderSerialNumber = it1->OrderNumber;

			//goodliffe::skip_list<OrderStruct>::iterator findIter = Order_struct.find(FilledOrderIn);
			std::map<uint64_t, uint32_t>::iterator findIter = LimitOrderBookLive.Order_struct.find(OrderNumber0);
			uint32_t LimitOrderSerialNumber = findIter->second;
			//LimitOrderBookLive.Order_struct.erase(findIter);

			uint32_t OrderSize1 = OrderSize0;
			char OrderType1 = LimitOrderBookLive.LimitOrderList[LimitOrderSerialNumber].OrderType;
			uint32_t OrderPrice1 = LimitOrderBookLive.LimitOrderList[LimitOrderSerialNumber].OrderPrice;
			if (OrderType1 == 'B')
			{
				TopOfBidBookBefore = LimitOrderBookLive.TopOfBidBook;
				if (OrderPrice1 != TopOfBidBookBefore)
				{
					std::cout << " There is an error. Market Order Does not match with Price; MarketExecute " << std::endl;
					std::cout << "  Limit Order Price is " << OrderPrice1 << std::endl;
					std::cout << "  Top Of Bid Book is " << TopOfBidBookBefore << std::endl;
				}
				LimitOrderBookLive.BidLimitOrderFillOrExecuteActionOnLimitOrderBook(OrderSize1);
				TopOfBidBookAfter = LimitOrderBookLive.TopOfBidBook;
				LimitOrderBookLive.LimitOrderExecuteActionOnLimitOrderList(LimitOrderSerialNumber, OrderSize1, TimeStamp0);

			}
			if (OrderType1 == 'S')
			{
				TopOfAskBookBefore = LimitOrderBookLive.TopOfAskBook;
				if (OrderPrice1 != TopOfAskBookBefore)
				{
					std::cout << " There is an error. Market Order Does not match with Price; Market Execute " << std::endl;
					std::cout << "  Limit Order Price is " << OrderPrice1 << std::endl;
					std::cout << "  Top Of Ask Book is " << TopOfAskBookBefore << std::endl;
				}
				LimitOrderBookLive.AskLimitOrderFillOrExecuteActionOnLimitOrderBook(OrderSize1);
				TopOfAskBookBefore = LimitOrderBookLive.TopOfAskBook;
				LimitOrderBookLive.LimitOrderExecuteActionOnLimitOrderList(LimitOrderSerialNumber, OrderSize1, TimeStamp0);
			}
		}
	}
	BidAndAskLevelSeries.close();

	return 0;
}

Amin · June 7th, 2021, 1:25 pm

Friends, we have been working for a long time with solution of SDEs and it seems reasonable that we need to at least achieve our benchmark of stochastic volatility option pricing with one SDE specified for asset price and another SDE specified for stochastic volatility and a description of correlation. I recently posted (post:1150) how to calculate variance of dt-integrals and that is an important step in our direction. We will be continuing on that theme for several months and I hope that we would be able to achieve stochastic volatility pricing out to 2-4 years with an arbitrary asset price SDE and an arbitrary stochastic volatility SDE and assuming a constant correlation between them.
I will keep friends posted with a lot of interesting material.

Amin · June 7th, 2021, 1:35 pm

Here is a sample Nasdaq ITCH file that I downloaded several years ago and I am attaching it for friends to be able to run it with my C++ program that I uploaded two posts ago.

20150109_EBAY.zip: (4.44 MiB) Downloaded 85 times

Amin · June 9th, 2021, 8:44 am

Friends, we have been working for a long time with solution of SDEs and it seems reasonable that we need to at least achieve our benchmark of stochastic volatility option pricing with one SDE specified for asset price and another SDE specified for stochastic volatility and a description of correlation. I recently posted (post:1150) how to calculate variance of dt-integrals and that is an important step in our direction. We will be continuing on that theme for several months and I hope that we would be able to achieve stochastic volatility pricing out to 2-4 years with an arbitrary asset price SDE and an arbitrary stochastic volatility SDE and assuming a constant correlation between them.
I will keep friends posted with a lot of interesting material.

Though I have several other ideas that I want to try(and see how they work), I want to start on correlated stochastic volatility option pricing by doing a two dimensional Ito-Taylor expansion with eight orders of expansion in each dimension(and crosses) with two arbitrary SDEs. I want two series expansions of which one will be for stochastic volatility in asset price and the other would be for square-root of stochastic variance in the asset price. I hope that this series solution could work for anywhere between .5 years to 2 years for SDEs that take gamma(volatility power) in stochastic volatility SDEs reasonably larger than .5. Since it is not simple to calculate the coefficients of these SDEs by hand, I will have to automate everything so it can take a bit of time (7-10 days) to do everything and verify its accuracy. I will keep friends posted about progress on my work.

Amin · June 13th, 2021, 6:32 pm

Friends, I had been thinking about organization of my new program but then I wanted to ask the question whether we could take a large step for stochastic volatility SDE using higher order asymptotic expansion. If we cannot take a large step with SV SDE with accuracy, then there is no point in applying it to the asset. And therefore I took my already written 8th order expansion (that I had also previously posted on this forum) and I realized that we can take a one year step size with 8th order expansion in original coordinates but not in bessel coordinates. For a larger stepping (with high order expansions) original coordinates faithfully capture the shape of distribution and far more so than expansions in bessel coordinates.
Now, I am very sure that we can easily take the step size in our 1-D quantizers-like simulation (that I like to call probability mass transfer method) to half a year, if we wish, by increasing the order of simulation but it may not be possible in Bessel coordinates easily. So it should also be relatively simple to take a higher order expansion in two dimensions(volatility and asset) in original coordinates and solve it using quantizer-like method (or probability mass transfer method).
Though tail and body of the mean reverting SDE distribution(even reasonably close to zero with good volatility) is perfectly captured out to one year by an 8th order expansion, I still see some minor problems in left tail and I will try increasing the order of expansion to 12 or 16 and see if I could fix the problems in left tail. In another day or two, I will post comparison graphs of 4th order, 8th order, 12th order and 16th order expansions for typical parameters close to zero and away from zero.

Amin · June 14th, 2021, 4:22 pm

I completed work on 12th order Ito-Taylor expansion of a single SDE and I will soon(in another day) posting my all my high order expansion programs so friends can play with and compare 4th, 8th and 12th order expansions in both original and bessel coordinates.
I will also be (hopefully tomorrow or in another day) posting monte carlo simulation programs for two different kinds of multi-dimensional correlated simulations.
a) A monte carlo system of correlated N SDEs of the kind we have in basket pricing. In this drift dynamics are single dimensional but brownian motions are correlated and are traditionally need cholesky or SVD of the correlation/variance-covariance matrix.
b) A system of two correlated SDEs of the kind we have in stochastic volatility option pricing simulations.
Since in high dimensions things quickly get complicated and number of terms increase very fast, for monte carlo simulations, I will simply present programs with 2nd order Ito-Taylor expansions.
I will first be posting monte carlo simulation programs for above multi-dmensional SDE structures and then we would move again towards analytics for stochastic volatility option pricing.

Amin · June 15th, 2021, 2:51 pm

Friends here is the program for SDE series expansion order 4-12. I have given some instructions at the start how to change order of the program.
.

function [] = NewHermiteExpansionO8()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)
%where mu1=kappa*theta and mu2=-kappa with
%beta1=0 and beta2=1;
%Instructions for changing the Expansion Order of the Program
%1. OrderA has to be given appropriate order value for the program
%This can be anywhere between 2-12%
%2. For Order four use the function (Look for the appropriate block below
%where these functions appear and uncomment the ones with right order and
%comment the other ones.)
%ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); for original coordinates
%%ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma); for Bessel coordinates
%For Order eight use the function
%ItoTaylorCoeffsNewOrder8(alpha,beta1,beta2,gamma); for original
%coordinates
%ItoTaylorCoeffsNewOrder8(alpha1,beta1,beta2,gamma); for Bessel coordinates
% For Order twelve use the function
%ItoTaylorCoeffsNewOrder12(alpha,beta1,beta2,gamma); for original
%coordinates
%ItoTaylorCoeffsNewOrder12(alpha1,beta1,beta2,gamma); for Bessel coordinates

%Both OrderA has to be given the right order value and functions for
%appropriate order have to be chosen as I previously explained.



dt=.5;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=1;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
OrderA=8;  %Keep at appropriate order.
OrderM=2;  %Keep at two.
T=Tt*dt;
TtM=8*2;%Monte carlo number of simulation intervals.
dtM=T/TtM;%Monte carlo time interval size dtM.



dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=50;  % No of normal density subdivisions
NnMidl=25;%One half left from mid of normal density(low)
NnMidh=26;%One half right from the mid of normal density(high)
NnMid=4.0;

x0=.150;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.60;   % volatility power.                                                                                                                                                                                                                                                                           
kappa=1.0;   %mean reversion parameter.
theta=.04;   %mean reversion target
%MonteTp5dtp03125x0p1Thp1Ka1Gap75Sip8

mu1=1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.60;%Volatility value
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;
                
w(1:Nn)=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

for nn=2:Nn-1
    ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);
    %Above calculate probability mass in each probability subdivision.
end

tic

wnStart=1;%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\1;

sigma11(1:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one. 
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderA+1)
    if sigma0~=0
        sigma11(k)=sigma0^(k-1);
    end
    if mu1 ~= 0
        mu11(k)=mu1^(k-1);
    end
    if mu2 ~= 0
        mu22(k)=mu2^(k-1);
    end
    if sigma0~=0
        sigma22(k)=sigma0^(2*(k-1));
    end
end
Ft(1:Tt+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA+1)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                Ft(l1,l2,l3,l4) = dt^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));

            end
        end
    end
end
% %Below call the program that calculates the Ito-Taylor coeffs to 12th order of expansion.
% YCoeff = ItoTaylorCoeffsNewOrder12(alpha,beta1,beta2,gamma);
% Yq=ItoTaylorCoeffsNewOrder12(alpha1,beta1,beta2,gamma);
% %Below call the program that calculates the Ito-Taylor coeffs to 8th order of expansion.
 YCoeff = ItoTaylorCoeffsNewOrder8(alpha,beta1,beta2,gamma);
 Yq=ItoTaylorCoeffsNewOrder8(alpha1,beta1,beta2,gamma);
%Below call the program that calculates the Ito-Taylor coeffs to 4th order of expansion.
%YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);

YqCoeff=Yq/(1-gamma);%These Bessel coordinates are divided by (1-gamma)
 %Calculate the density with monte carlo below.

  HermiteP(1,1:Nn)=1;
  HermiteP(2,1:Nn)=Z(1:Nn);
  HermiteP(3,1:Nn)=Z(1:Nn).^2-1;
  HermiteP(4,1:Nn)=Z(1:Nn).^3-3*Z(1:Nn);
  HermiteP(5,1:Nn)=Z(1:Nn).^4-6*Z(1:Nn).^2+3;
  HermiteP(6,1:Nn)=Z(1:Nn).^5-10*Z(1:Nn).^3+15*Z(1:Nn);
  HermiteP(7,1:Nn)=Z(1:Nn).^6-15*Z(1:Nn).^4+45*Z(1:Nn).^2-15;
  HermiteP(8,1:Nn)=Z(1:Nn).^7-21*Z(1:Nn).^5+105*Z(1:Nn).^3-105*Z(1:Nn);
  HermiteP(9,1:Nn)=Z(1:Nn).^8-28*Z(1:Nn).^6+210*Z(1:Nn).^4-420*Z(1:Nn).^2+105;
  HermiteP(10,1:Nn)=Z(1:Nn).^9-36*Z(1:Nn).^7+378*Z(1:Nn).^5-1260*Z(1:Nn).^3+945*Z(1:Nn);
  HermiteP(11,1:Nn)=Z(1:Nn).^10-45*Z(1:Nn).^8+630*Z(1:Nn).^6-3150*Z(1:Nn).^4+4725*Z(1:Nn).^2-945;
  HermiteP(12,1:Nn)=Z(1:Nn).*HermiteP(11,1:Nn)-10*HermiteP(10,1:Nn);
  HermiteP(13,1:Nn)=Z(1:Nn).*HermiteP(12,1:Nn)-11*HermiteP(11,1:Nn);



     yy0=x0;
     ww0=x0^(1-gamma)/(1-gamma);
     yy(1:Nn)=x0;  %Original process monte carlo.
     ww(1:Nn)=ww0;
     for k = 0 : OrderA
         for m = 0:k
             l4 = k - m + 1;
             for n = 0 : m
                 l3 = m - n + 1;
                 for j = 0:n
                     l2 = n - j + 1;
                     l1 = j + 1;
                     if(l4<=5)   
                         
                      ww(1:Nn)=ww(1:Nn) + YqCoeff(l1,l2,l3,l4) * ...
                     mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                     ((1-gamma)*ww0).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
                     HermiteP(l4,1:Nn) * Ft(l1,l2,l3,l4);    
                 
                     yy(1:Nn)=yy(1:Nn) + YCoeff(l1,l2,l3,l4) * ...
                     mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                     yy0.^Fp(l1,l2,l3,l4) .* HermiteP(l4,1:Nn) * Ft(l1,l2,l3,l4);
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end

 y_w(1:Nn)=0;
 y_w(wnStart:Nn) = ((1-gamma)*ww(1:Nn)).^(1/(1-gamma));
 
 Dfy_w(wnStart:Nn)=0;
 Dfy(wnStart:Nn)=0;
 
for nn=wnStart+1:Nn-1
    Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    
    Dfy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    %Change of variable derivative for densities
end
fy_w(1:Nn)=0;

fy(1:Nn)=0;
for nn = wnStart:Nn-1
    fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
    
    fy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy(nn));
end
toc
%str=input('Analytic Values calculated; Press a key to start monte carlo');



sigma11(1:OrderM+1)=0;
mu11(1:OrderM+1)=0;
mu22(1:OrderM+1)=0;
sigma22(1:OrderM+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one. 
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderM+1)
    if sigma0~=0
        sigma11(k)=sigma0^(k-1);
    end
    if mu1 ~= 0
        mu11(k)=mu1^(k-1);
    end
    if mu2 ~= 0
        mu22(k)=mu2^(k-1);
    end
    if sigma0~=0
        sigma22(k)=sigma0^(2*(k-1));
    end
end
Ft(1:TtM+1,1:(OrderM+1),1:(OrderM+1),1:(OrderM+1),1:(OrderM+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderM+1),1:(OrderM+1),1:(OrderM+1),1:(OrderM+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderM+1),1:(OrderM+1),1:(OrderM+1),1:(OrderM+1))=0;%General x powers for bessel transformed coordinates.
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderM+1)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));

            end
        end
    end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
%YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)
 %Calculate the density with monte carlo below.
 rng(85109634, 'twister')
 paths=250000;
 YY(1:paths)=x0;  %Original process monte carlo.
% YY2dt(1:paths)=0.0;
 %YY4dt(1:paths)=0.0;
 X(1:paths)=x0^(1-gamma)/(1-gamma);
 Random1(1:paths)=0;
 for tt=1:TtM
 %    t1=tt*dtM;
 %    t0=(tt-1)*dtM;
     Random1=randn(size(Random1));
     HermiteP1(1,1:paths)=1;
     HermiteP1(2,1:paths)=Random1(1:paths);
     HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
     HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
     HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
     YY0=YY;
  %   YYPrev=YY;
  %   XX0=X;
     
     for k = 0 : OrderM
         for m = 0:k
             l4 = k - m + 1;
             for n = 0 : m
                 l3 = m - n + 1;
                 for j = 0:n
                     l2 = n - j + 1;
                     l1 = j + 1;
                     if(l4<=5)   
                         
  %                    X(1:paths)=X(1:paths) + YqCoeff(l1,l2,l3,l4) * ...
  %                   mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
  %                   ((1-gamma)*XX0(1:paths)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
  %                   HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);    
                 
                     YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
                     mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                     YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
 end
 YY(YY<0)=0;
 
% sum(YY(:))/paths %origianl coordinates monte carlo average.
% sum(y_w(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
% sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
TrueMean= theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
MCVar=sum((YY(:)-MCMean).^2)/paths
disp('Original process average from our simulation');
ItoHermiteMeanBes=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
ItoHermiteMeanOrig=sum(yy(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
ItoHermiteVarBes=sum((y_w(wnStart+1:Nn-1)-ItoHermiteMeanBes).^2.*ZProb(wnStart+1:Nn-1)) 
ItoHermiteVarOrig=sum((yy(wnStart+1:Nn-1)-ItoHermiteMeanOrig).^2.*ZProb(wnStart+1:Nn-1)) 
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 
 MaxCutOff=30;
NoOfBins=round(300*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );

%plot(y_w1(wnStart+1:Nn-1),py_w1(wnStart+1:Nn-1),'b',y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
plot(yy(wnStart+1:Nn-1),fy(wnStart+1:Nn-1),'r',y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(yy(wnStart+1:Nn-1),fy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');


title(sprintf('x0 = %.4f,theta=%.4f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.3f,MOrig=%.5f,MBes=%.5f,TM=%.5f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanOrig,ItoHermiteMeanBes,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'OriginalCoordinates Density','Bessel Coordinates Density','Monte Carlo Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
end

.
Here are the functions for eighth and twelfth order
.
.

function [IntegralCoeff] = ComputeIntegralCoeffsOrder12()


IntegralCoeff(1:3,1:3,1:3,1:3,1:3,1:3,1:3,1:3,1:3,1:3,1:3,1:3)=0;
%IntegralCoeffdt(1:3,1:3,1:3,1:3)=0;
%IntegralCoeffdz(1:3,1:3,1:3,1:3)=0;

IntegralCoeff0(1:3,1:3,1:3,1:3)=0;

IntegralCoeff0(1,1,1,2)=1;
IntegralCoeff0(1,1,1,3)=1;

IntegralCoeff0(1,1,2,2)=1/2;
IntegralCoeff0(1,1,3,2)=1-1/sqrt(3);
IntegralCoeff0(1,1,2,3)=1/sqrt(3);
IntegralCoeff0(1,1,3,3)=1/2;

IntegralCoeff0(1,2,2,2)=1/6;
IntegralCoeff0(1,3,2,2)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5));
IntegralCoeff0(1,2,3,2)=1/sqrt(3)/2*(1-1/sqrt(5));
IntegralCoeff0(1,3,3,2)=1/2*(1-sqrt(2)/2);
IntegralCoeff0(1,2,2,3)=1/2*1/sqrt(5);
IntegralCoeff0(1,3,2,3)=(1-1/sqrt(3))*1/sqrt(2)*1/2;
IntegralCoeff0(1,2,3,3)=1/sqrt(3)/sqrt(2)*1/(2);
IntegralCoeff0(1,3,3,3)=1/6;

IntegralCoeff0(2,2,2,2)=1/24;
IntegralCoeff0(2,2,3,2)=1/2*1/sqrt(5)*1/3*(1-1/sqrt(7));
IntegralCoeff0(2,3,2,2)=1/sqrt(3)*1/2*(1-1/sqrt(5))*1/3*(1-1/sqrt(7));
IntegralCoeff0(3,2,2,2)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5))*1/3*(1-1/sqrt(7));
IntegralCoeff0(3,3,2,2)=1/2*(1-sqrt(2)/2)*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(2,3,3,2)=1/sqrt(3)*1/sqrt(2)*1/2*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(3,2,3,2)=(1-1/sqrt(3))*1/sqrt(2)*1/2*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(3,3,3,2)=1/6*(1-sqrt(3)/sqrt(5));

IntegralCoeff0(2,2,2,3)=1/6*1/sqrt(7);
IntegralCoeff0(2,2,3,3)=1/2*1/sqrt(5)*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(2,3,2,3)=1/sqrt(3)*1/2*(1-1/sqrt(5))*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(3,2,2,3)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5))*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(3,3,2,3)=1/2*(1-sqrt(2)/2)*1/sqrt(3)*1/sqrt(5);
IntegralCoeff0(2,3,3,3)=1/sqrt(3)*1/sqrt(2)*1/sqrt(4)*1/sqrt(3)*(1/sqrt(5));
IntegralCoeff0(3,2,3,3)=(1-1/sqrt(3))*1/sqrt(2)*1/2*1/sqrt(3)*1/sqrt(5);
IntegralCoeff0(3,3,3,3)=1/24;


%Can also be calculated with algorithm below.
%here IntegralCoeffdt indicates the coefficients of a dt-integral.
%This dt-integral means that you are calculating the Ito-hermite expansion
%of  Integral X(t)^alpha dt with X(t) dynamics given by the SDE 
%here IntegralCoeffdz indicates the coefficients of a dz-integral.
%This dz-integral means that you are calculating the Ito-hermite expansion
%of  Integral X(t)^alpha dz(t) with X(t) dynamics given by the SDE 

%IntegralCoeff is is associated with expansion of X(t)^alpha.
%IntegralCoeff below is not returned and IntegralCoeff0 manually calculated
%above is returned but both are the same.

l0(1:2)=1;

for m12=1:2
    l0(1)=1;
    l0(2)=1;
    
    %IntegralCoeff4(m4,1,1,1)=1;
    %IntegralCoeff4(m4,1,1,1)=1;
    %1 is added to m8 since index 1 stands for zero, 2 for one and three
    %for two.
    IntegralCoeff(1,1,1,1,1,1,1,1,1,1,1,m12+1)=1;
    for m11=1:2
        l0(1)=1;
        l0(2)=1;
        l0(m12)=l0(m12)+1;
        
        if(m11==1)
            IntegralCoeff(1,1,1,1,1,1,1,1,1,1,m12+1,m11+1)=1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
        end
        if(m11==2)
            IntegralCoeff(1,1,1,1,1,1,1,1,1,1,m12+1,m11+1)= 1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
        end
        for m10=1:2
            l0(1)=1;
            l0(2)=1;
            l0(m12)=l0(m12)+1;
            l0(m11)=l0(m11)+1;
            
            if(m10==1)
                IntegralCoeff(1,1,1,1,1,1,1,1,1,m12+1,m11+1,m10+1)=IntegralCoeff(1,1,1,1,1,1,1,1,1,1,m12+1,m11+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
            end
            if(m10==2)
                IntegralCoeff(1,1,1,1,1,1,1,1,1,m12+1,m11+1,m10+1)=IntegralCoeff(1,1,1,1,1,1,1,1,1,1,m12+1,m11+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
            end

            for m9=1:2
                l0(1)=1;
                l0(2)=1;
                l0(m12)=l0(m12)+1;
                l0(m11)=l0(m11)+1;
                l0(m10)=l0(m10)+1;
               
                if(m9==1)
                    IntegralCoeff(1,1,1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1)=IntegralCoeff(1,1,1,1,1,1,1,1,1,m12+1,m11+1,m10+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                end
                if(m9==2)
                    IntegralCoeff(1,1,1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1)=IntegralCoeff(1,1,1,1,1,1,1,1,1,m12+1,m11+1,m10+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                end
                for m8=1:2
                    l0(1)=1;
                    l0(2)=1;
                    l0(m12)=l0(m12)+1;
                    l0(m11)=l0(m11)+1;
                    l0(m10)=l0(m10)+1;
                    l0(m9)=l0(m9)+1;
                    if(m8==1)
                        IntegralCoeff(1,1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1)=IntegralCoeff(1,1,1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                    end
                    if(m8==2)
                        IntegralCoeff(1,1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1)=IntegralCoeff(1,1,1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                    end
                    for m7=1:2
                        l0(1)=1;
                        l0(2)=1;
                        l0(m12)=l0(m12)+1;
                        l0(m11)=l0(m11)+1;
                        l0(m10)=l0(m10)+1;
                        l0(m9)=l0(m9)+1;
                        l0(m8)=l0(m8)+1;
                        if(m7==1)
                            IntegralCoeff(1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1)=IntegralCoeff(1,1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                        end
                        if(m7==2)
                            IntegralCoeff(1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1)=IntegralCoeff(1,1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                        end
                        for m6=1:2
                            l0(1)=1;
                            l0(2)=1;
                            l0(m12)=l0(m12)+1;
                            l0(m11)=l0(m11)+1;
                            l0(m10)=l0(m10)+1;
                            l0(m9)=l0(m9)+1;
                            l0(m8)=l0(m8)+1;
                            l0(m7)=l0(m7)+1;
                            if(m6==1)
                                IntegralCoeff(1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1)=IntegralCoeff(1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                            end
                            if(m6==2)
                                IntegralCoeff(1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1)=IntegralCoeff(1,1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                            end
                            for m5=1:2
                                l0(1)=1;
                                l0(2)=1;
                                l0(m12)=l0(m12)+1;
                                l0(m11)=l0(m11)+1;
                                l0(m10)=l0(m10)+1;
                                l0(m9)=l0(m9)+1;
                                l0(m8)=l0(m8)+1;
                                l0(m7)=l0(m7)+1;
                                l0(m6)=l0(m6)+1;
                                
                                if(m5==1)
                                    IntegralCoeff(1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1)=IntegralCoeff(1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                                end
                                if(m5==2)
                                    IntegralCoeff(1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1)=IntegralCoeff(1,1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                                end
                                for m4=1:2
                                    l0(1)=1;
                                    l0(2)=1;
                                    l0(m12)=l0(m12)+1;
                                    l0(m11)=l0(m11)+1;
                                    l0(m10)=l0(m10)+1;
                                    l0(m9)=l0(m9)+1;
                                    l0(m8)=l0(m8)+1;
                                    l0(m7)=l0(m7)+1;
                                    l0(m6)=l0(m6)+1;
                                    l0(m5)=l0(m5)+1;
                                    if(m4==1)
                                        IntegralCoeff(1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1)=IntegralCoeff(1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                                    end
                                    if(m4==2)
                                        IntegralCoeff(1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1)=IntegralCoeff(1,1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                                    end
                                    for m3=1:2
                                        l0(1)=1;
                                        l0(2)=1;
                                        l0(m12)=l0(m12)+1;
                                        l0(m11)=l0(m11)+1;
                                        l0(m10)=l0(m10)+1;
                                        l0(m9)=l0(m9)+1;
                                        l0(m8)=l0(m8)+1;
                                        l0(m7)=l0(m7)+1;
                                        l0(m6)=l0(m6)+1;
                                        l0(m5)=l0(m5)+1;
                                        l0(m4)=l0(m4)+1;
                                        if(m3==1)
                                            IntegralCoeff(1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1)=IntegralCoeff(1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                                        end
                                        if(m3==2)
                                            IntegralCoeff(1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1)=IntegralCoeff(1,1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                                        end
                                        for m2=1:2
                                            l0(1)=1;
                                            l0(2)=1;
                                            l0(m12)=l0(m12)+1;
                                            l0(m11)=l0(m11)+1;
                                            l0(m10)=l0(m10)+1;
                                            l0(m9)=l0(m9)+1;
                                            l0(m8)=l0(m8)+1;
                                            l0(m7)=l0(m7)+1;
                                            l0(m6)=l0(m6)+1;
                                            l0(m5)=l0(m5)+1;
                                            l0(m4)=l0(m4)+1;
                                            l0(m3)=l0(m3)+1;
                                            if(m2==1)
                                                IntegralCoeff(1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1)=IntegralCoeff(1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                                            end
                                            if(m2==2)
                                                IntegralCoeff(1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1)=IntegralCoeff(1,1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                                            end
                                            for m1=1:2
                                                l0(1)=1;
                                                l0(2)=1;
                                                l0(m12)=l0(m12)+1;
                                                l0(m11)=l0(m11)+1;
                                                l0(m10)=l0(m10)+1;
                                                l0(m9)=l0(m9)+1;
                                                l0(m8)=l0(m8)+1;
                                                l0(m7)=l0(m7)+1;
                                                l0(m6)=l0(m6)+1;
                                                l0(m5)=l0(m5)+1;
                                                l0(m4)=l0(m4)+1;
                                                l0(m3)=l0(m3)+1;
                                                l0(m2)=l0(m2)+1;
                                                if(m1==1)
                                                    IntegralCoeff(m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                                                end
                                                if(m1==2)
                                                    IntegralCoeff(m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(1,m12+1,m11+1,m10+1,m9+1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                                                end
                                            end
                                        end
                                    end

                                end
                                
                            end
                        end
                    end
                end
            end

        end
    end
end
end

.
.

function [Y] = ItoTaylorCoeffsNewOrder12(alpha,beta1,beta2,gamma)


%In the coefficient calculation program which calculates Y(l1,l2,l3,l4), 
%I have used four levels of looping each for relevant expansion order. 
%The first loop takes four values and second loop takes 16 values and 
%third loop takes 64 values and so on. And then each coefficient 
%term can be individually calculated while carefully accounting 
%for path dependence. 
%So for example in a nested loop structure 
%m1= 1:mDim
% m2=1:mDim
% m3=1:mDim
%    l(m1)=l(m1)+1;
%    l(m2)=l(m2)+1;
%    l(m3)=l(m3)+1;

%in the above looping loop takes values from one to four with one 
%indicating the first drift term, two indicating the second drift term 
%and three indicating quadratic variation term and 
%four indicating the volatility term. And with this looping structure 
%we can So in the above looping m1=1 would mean that all terms are 
%descendent of first drift term and m2=4 would mean that all terms are 
%descendent of first drift term on first expansion order and descendent 
%of volatility term on the second order and so we can keep track of path 
%dependence perfectly. 
%And on each level, we individually calculate the descendent terms. While 
%keeping track of path dependence and calculating the coefficients with 
%careful path dependence consideration, we update the appropriate element 
%in our polynomial like expansion coefficient array 

%explaining the part of code
%m1= 1:mDim
% m2=1:mDim
% m3=1:mDim
%    l(m1)=l(m1)+1;
%    l(m2)=l(m2)+1;
%    l(m3)=l(m3)+1;
%Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);

%Here l(1) denotes l1 but written as l(1) so it can be conveniently 
%updated with the loop variable when the loop variable takes value one 
%indicating first drift term . And l(2) could be conveniently updated when 
%the loop variable takes value equal to two indicating second 
%drift term and so on.
%Here is the part of code snippet for that

%for m1=1:mDim
%    l(1)=1;
%    l(2)=1;
%    l(3)=1;
%    l(4)=1;
%    l(m1)=l(m1)+1;
%CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
%    CoeffDX2 = CoeffDX1 - 1;
%    ArrIndex0=m1;
%    ArrIndex=(m1-1)*mDim;
%    Coeff1st=Y1(ArrIndex0)*CoeffDX1;
%    Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%    Y2(1+ArrIndex)=Coeff1st;
%    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
%    Y2(2+ArrIndex)=Coeff1st;
%    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
%    Y2(3+ArrIndex)=Coeff2nd;
%    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,2,n1(m1));
%    Y2(4+ArrIndex)=Coeff1st;
%    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,3,n1(m1));

 
%The first four lines update the array indices according to the parent term. 
%And then CoeffDX1 and CoeffDX2 are calculated according to algebraic exponents on parent terms.    
%ArrIndex0=m1; calculates the array index of the parent term 
%And ArrIndex=(m1-1)*mDim; calculates the array index of the descendent terms
%And coefficient of the drift and volatility descendent terms is 
%calculated by multiplying the coefficient of the parent term by 
%Coeff1st=Y1(ArrIndex0)*CoeffDX1;
%And coefficient of the quadratic variation descendent terms is 
%calculated by multiplying the coefficient of the parent term by 
%Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%And then each of the four descendent terms are updated with Coeff1st 
%if they are drift or volatility descendent terms or Coeff2nd if 
%they are quadratic variation descendent terms.
%Here Y1 indicates the temporary coefficient array with parent terms on 
%first level. And Y2 denotes temporary coefficient array with parent terms 
%on second level and Y3 indicates temporary coefficient array with parent terms on third level.

[IntegralCoeff] = ComputeIntegralCoeffsOrder12();
%[IntegralCoeff1] = ComputeIntegralCoeffsOrder8();
%IntegralCoeff(1,1,1,1,:,:,:,:,:,:,:,:)=IntegralCoeff1(:,:,:,:,:,:,:,:);

%IntegralCoeff2(1,1,1,1,:,:,:,:,:,:,:,:)-IntegralCoeff(1,1,1,1,:,:,:,:,:,:,:,:)
% for m1=1:3
%     for m2=1:3
%         for m3=1:3
%             for m4=1:3
%                 for m5=1:3
%                     for m6=1:3
%                         for m7=1:3
%                             for m8=1:3
%                                 if(IntegralCoeff2(1,1,1,1,m8,m7,m6,m5,m4,m3,m2,m1) ~= IntegralCoeff(1,1,1,1,m8,m7,m6,m5,m4,m3,m2,m1))
%                                     IntegralCoeff2(1,1,1,1,m8,m7,m6,m5,m4,m3,m2,m1)- IntegralCoeff(1,1,1,1,m8,m7,m6,m5,m4,m3,m2,m1)
%                                     IntegralCoeff2(1,1,1,1,m8,m7,m6,m5,m4,m3,m2,m1)
%                                     IntegralCoeff(1,1,1,1,m8,m7,m6,m5,m4,m3,m2,m1)
%                                     m8
%                                     m7
%                                     m6
%                                     m5
%                                     m4
%                                     m3
%                                     m2
%                                     m1
%                                     str=input('Look at difference of integrals');
%                                 end
%                             end
%                         end
%                     end
%                 end
%             end
%         end
%     end
% end


%str=input('Look at difference of integrals');


n1(1)=2;
n1(2)=2;
n1(3)=2;
n1(4)=3;
%n1(1), n1(2), n1(3) are drift and quadratic variation term variables 
%and take a value equal to 2 indicating a dt integral.
%n1(4) is volatility term variable and indicates a dz-integral by taking a
%value of 3. 

mDim=4; % four descendent terms in each expansion
Y(1:13,1:13,1:13,1:13)=0;
Y1(1:mDim)=0;
Y2(1:mDim*mDim)=0;
Y3(1:mDim*mDim*mDim)=0;
Y4(1:mDim*mDim*mDim*mDim)=0;
Y5(1:mDim*mDim*mDim*mDim*mDim)=0;
Y6(1:mDim*mDim*mDim*mDim*mDim*mDim)=0;
Y7(1:mDim*mDim*mDim*mDim*mDim*mDim*mDim)=0;
Y8(1:mDim*mDim*mDim*mDim*mDim*mDim*mDim*mDim)=0;
Y9(1:mDim*mDim*mDim*mDim*mDim*mDim*mDim*mDim*mDim)=0;
Y10(1:mDim*mDim*mDim*mDim*mDim*mDim*mDim*mDim*mDim*mDim)=0;
Y11(1:mDim*mDim*mDim*mDim*mDim*mDim*mDim*mDim*mDim*mDim*mDim)=0;


%First Ito-hermite expansion level starts here. No loops but four
%descendent terms.
l(1)=1;
l(2)=1;
l(3)=1;
l(4)=1;
CoeffDX1 = alpha;
CoeffDX2 = CoeffDX1 - 1;
Coeff1st=CoeffDX1;
Coeff2nd=.5*CoeffDX1*CoeffDX2;

Y1(1)=Coeff1st;
Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,1,1,1,2);
Y1(2)=Coeff1st;
Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,1,1,1,2);
Y1(3)=Coeff2nd;
Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,1,1,1,1,1,1,1,2);
Y1(4)=Coeff1st;
Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,1,1,1,3);

%Second Ito-hermite expansion level starts. It has a loop over four parent
%terms and there are four descendent terms for each parent term.
%The coefficient terms are then lumped in a polynomial-like expansion
%array of coefficients.
for m1=1:mDim
    l(1)=1;
    l(2)=1;
    l(3)=1;
    l(4)=1;
    l(m1)=l(m1)+1;
    CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
    CoeffDX2 = CoeffDX1 - 1;
    ArrIndex0=m1;
    ArrIndex=(m1-1)*mDim;
    Coeff1st=Y1(ArrIndex0)*CoeffDX1;
    Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
    Y2(1+ArrIndex)=Coeff1st;
    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,1,1,2,n1(m1));
    Y2(2+ArrIndex)=Coeff1st;
    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,1,1,2,n1(m1));
    Y2(3+ArrIndex)=Coeff2nd;
    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,1,1,1,1,1,1,2,n1(m1));
    Y2(4+ArrIndex)=Coeff1st;
    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,1,1,3,n1(m1));
    %Third Ito-hermite expansion level starts and it is a nested loop with
    %a total of sixteen parents and each parent takes four descendent
    %terms.
    for m2=1:mDim
        l(1)=1;
        l(2)=1;
        l(3)=1;
        l(4)=1;
        l(m1)=l(m1)+1;
        l(m2)=l(m2)+1;
        CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
        CoeffDX2=CoeffDX1-1;            
        ArrIndex0=(m1-1)*mDim+m2;
        ArrIndex=((m1-1)*mDim+(m2-1))*mDim;
        Coeff1st=Y2(ArrIndex0)*CoeffDX1;
        Coeff2nd=Y2(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
        Y3(1+ArrIndex)=Coeff1st;
        Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,1,2,n1(m2),n1(m1));
        Y3(2+ArrIndex)=Coeff1st;
        Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,1,2,n1(m2),n1(m1));
        Y3(3+ArrIndex)=Coeff2nd;
        Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,1,1,1,1,1,2,n1(m2),n1(m1));
        Y3(4+ArrIndex)=Coeff1st;
        Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,1,3,n1(m2),n1(m1));
        %fourht Ito-hermite expansion level starts and it is a triply-nested loop with
        %a total of sixteen parents and each parent takes four descendent
        %terms. We then lump the terms in a relatively sparse polynomial 
        %like expansion coefficient array that has smaller number of
        %non-zero terms.

        for m3=1:mDim
            l(1)=1;
            l(2)=1;
            l(3)=1;
            l(4)=1;
            l(m1)=l(m1)+1;
            l(m2)=l(m2)+1;
            l(m3)=l(m3)+1;
            CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
            CoeffDX2=CoeffDX1-1;
            ArrIndex0=((m1-1)*mDim+(m2-1))*mDim+m3;
            ArrIndex=(((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim;
            Coeff1st=Y3(ArrIndex0)*CoeffDX1;
            Coeff2nd=Y3(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
            Y4(1+ArrIndex)=Coeff1st;
            Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,2,n1(m3),n1(m2),n1(m1));
            Y4(2+ArrIndex)=Coeff1st;
            Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,2,n1(m3),n1(m2),n1(m1));
            Y4(3+ArrIndex)=Coeff2nd;
            Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,1,1,1,1,2,n1(m3),n1(m2),n1(m1));
            Y4(4+ArrIndex)=Coeff1st;
            Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,1,3,n1(m3),n1(m2),n1(m1));
            for m4=1:mDim
                l(1)=1;
                l(2)=1;
                l(3)=1;
                l(4)=1;
                l(m1)=l(m1)+1;
                l(m2)=l(m2)+1;
                l(m3)=l(m3)+1;
                l(m4)=l(m4)+1;
                CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                CoeffDX2=CoeffDX1-1;
                
                ArrIndex0=(((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+m4;
                ArrIndex=((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim;
                
                Coeff1st=Y4(ArrIndex0)*CoeffDX1;
                Coeff2nd=Y4(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                Y5(1+ArrIndex)=Coeff1st;
                Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,2,n1(m4),n1(m3),n1(m2),n1(m1));
                Y5(2+ArrIndex)=Coeff1st;
                Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,2,n1(m4),n1(m3),n1(m2),n1(m1));
                Y5(3+ArrIndex)=Coeff2nd;
                Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,1,1,1,2,n1(m4),n1(m3),n1(m2),n1(m1));
                Y5(4+ArrIndex)=Coeff1st;
                Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,3,n1(m4),n1(m3),n1(m2),n1(m1));
                for m5=1:mDim
                    l(1)=1;
                    l(2)=1;
                    l(3)=1;
                    l(4)=1;
                    l(m1)=l(m1)+1;
                    l(m2)=l(m2)+1;
                    l(m3)=l(m3)+1;
                    l(m4)=l(m4)+1;
                    l(m5)=l(m5)+1;
                    CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                    CoeffDX2=CoeffDX1-1;
                
                    ArrIndex0=((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+m5;
                    ArrIndex=(((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim;
                
                    Coeff1st=Y5(ArrIndex0)*CoeffDX1;
                    Coeff2nd=Y5(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                    Y6(1+ArrIndex)=Coeff1st;
                    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,2,n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                    Y6(2+ArrIndex)=Coeff1st;
                    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,2,n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                    Y6(3+ArrIndex)=Coeff2nd;
                    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,1,1,2,n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                    Y6(4+ArrIndex)=Coeff1st;
                    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,1,1,3,n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                    for m6=1:mDim
                        l(1)=1;
                        l(2)=1;
                        l(3)=1;
                        l(4)=1;
                        l(m1)=l(m1)+1;
                        l(m2)=l(m2)+1;
                        l(m3)=l(m3)+1;
                        l(m4)=l(m4)+1;
                        l(m5)=l(m5)+1;
                        l(m6)=l(m6)+1;
                        CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                            - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                        CoeffDX2=CoeffDX1-1;
                
                        ArrIndex0=(((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+m6;
                        ArrIndex=((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim;
                
                        Coeff1st=Y6(ArrIndex0)*CoeffDX1;
                        Coeff2nd=Y6(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                        Y7(1+ArrIndex)=Coeff1st;
                        Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,2,n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                        Y7(2+ArrIndex)=Coeff1st;
                        Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,2,n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                        Y7(3+ArrIndex)=Coeff2nd;
                        Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,1,2,n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                        Y7(4+ArrIndex)=Coeff1st;
                        Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,1,3,n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                        for m7=1:mDim
                            l(1)=1;
                            l(2)=1;
                            l(3)=1;
                            l(4)=1;
                            l(m1)=l(m1)+1;
                            l(m2)=l(m2)+1;
                            l(m3)=l(m3)+1;
                            l(m4)=l(m4)+1;
                            l(m5)=l(m5)+1;
                            l(m6)=l(m6)+1;
                            l(m7)=l(m7)+1;
                            CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                                - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                            CoeffDX2=CoeffDX1-1;
                
                            ArrIndex0=((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+m7;
                            ArrIndex=(((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+(m7-1))*mDim;
                
                            Coeff1st=Y7(ArrIndex0)*CoeffDX1;
                            Coeff2nd=Y7(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                            Y8(1+ArrIndex)=Coeff1st;
                            Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,2,n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                            Y8(2+ArrIndex)=Coeff1st;
                            Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,2,n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                            Y8(3+ArrIndex)=Coeff2nd;
                            Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,2,n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                            Y8(4+ArrIndex)=Coeff1st;
                            Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,3,n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                        
                            for m8=1:mDim
                                l(1)=1;
                                l(2)=1;
                                l(3)=1;
                                l(4)=1;
                                l(m1)=l(m1)+1;
                                l(m2)=l(m2)+1;
                                l(m3)=l(m3)+1;
                                l(m4)=l(m4)+1;
                                l(m5)=l(m5)+1;
                                l(m6)=l(m6)+1;
                                l(m7)=l(m7)+1;
                                l(m8)=l(m8)+1;
                                CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                                    - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                                CoeffDX2=CoeffDX1-1;
                
                                ArrIndex0=(((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+(m7-1))*mDim+m8;
                                ArrIndex=((((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+(m7-1))*mDim+(m8-1))*mDim;
                
                                Coeff1st=Y8(ArrIndex0)*CoeffDX1;
                                Coeff2nd=Y8(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                                Y9(1+ArrIndex)=Coeff1st;
                                Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2,n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                Y9(2+ArrIndex)=Coeff1st;
                                Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2,n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                Y9(3+ArrIndex)=Coeff2nd;
                                Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,2,n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                Y9(4+ArrIndex)=Coeff1st;
                                Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,3,n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                for m9=1:mDim
                                    l(1)=1;
                                    l(2)=1;
                                    l(3)=1;
                                    l(4)=1;
                                    l(m1)=l(m1)+1;
                                    l(m2)=l(m2)+1;
                                    l(m3)=l(m3)+1;
                                    l(m4)=l(m4)+1;
                                    l(m5)=l(m5)+1;
                                    l(m6)=l(m6)+1;
                                    l(m7)=l(m7)+1;
                                    l(m8)=l(m8)+1;
                                    l(m9)=l(m9)+1;
                                    CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                                    CoeffDX2=CoeffDX1-1;
                
                                     ArrIndex0=((((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+(m7-1))*mDim+(m8-1))*mDim+m9;
                                    ArrIndex=(((((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+(m7-1))*mDim+(m8-1))*mDim+(m9-1))*mDim;
                
                                    Coeff1st=Y9(ArrIndex0)*CoeffDX1;
                                    Coeff2nd=Y9(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                                    Y10(1+ArrIndex)=Coeff1st;
                                    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                    Y10(2+ArrIndex)=Coeff1st;
                                    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                    Y10(3+ArrIndex)=Coeff2nd;
                                    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,2,n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                    Y10(4+ArrIndex)=Coeff1st;
                                    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,3,n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                    for m10=1:mDim
                                        l(1)=1;
                                        l(2)=1;
                                        l(3)=1;
                                        l(4)=1;
                                        l(m1)=l(m1)+1;
                                        l(m2)=l(m2)+1;
                                        l(m3)=l(m3)+1;
                                        l(m4)=l(m4)+1;
                                        l(m5)=l(m5)+1;
                                        l(m6)=l(m6)+1;
                                        l(m7)=l(m7)+1;
                                        l(m8)=l(m8)+1;
                                        l(m9)=l(m9)+1;
                                        l(m10)=l(m10)+1;
                                        CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                                            - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                                        CoeffDX2=CoeffDX1-1;
                
                                        ArrIndex0=(((((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+(m7-1))*mDim+(m8-1))*mDim+(m9-1))*mDim+m10;
                                        ArrIndex=((((((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+(m7-1))*mDim+(m8-1))*mDim+(m9-1))*mDim+(m10-1))*mDim;
                
                                        Coeff1st=Y10(ArrIndex0)*CoeffDX1;
                                        Coeff2nd=Y10(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                                        Y11(1+ArrIndex)=Coeff1st;
                                        Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,2,n1(m10),n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                        Y11(2+ArrIndex)=Coeff1st;
                                        Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,2,n1(m10),n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                        Y11(3+ArrIndex)=Coeff2nd;
                                        Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,2,n1(m10),n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                        Y11(4+ArrIndex)=Coeff1st;
                                        Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,3,n1(m10),n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                        for m11=1:mDim
                                            l(1)=1;
                                            l(2)=1;
                                            l(3)=1;
                                            l(4)=1;
                                            l(m1)=l(m1)+1;
                                            l(m2)=l(m2)+1;
                                            l(m3)=l(m3)+1;
                                            l(m4)=l(m4)+1;
                                            l(m5)=l(m5)+1;
                                            l(m6)=l(m6)+1;
                                            l(m7)=l(m7)+1;
                                            l(m8)=l(m8)+1;
                                            l(m9)=l(m9)+1;
                                            l(m10)=l(m10)+1;
                                            l(m11)=l(m11)+1;
                                            CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                                                - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                                            CoeffDX2=CoeffDX1-1;
                
                                            ArrIndex0=((((((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+(m7-1))*mDim+(m8-1))*mDim+(m9-1))*mDim+(m10-1))*mDim+m11;
                                            ArrIndex=(((((((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+(m7-1))*mDim+(m8-1))*mDim+(m9-1))*mDim+(m10-1))*mDim+(m11-1))*mDim;
                
                                            Coeff1st=Y11(ArrIndex0)*CoeffDX1;
                                            Coeff2nd=Y11(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                                            %Y12(1+ArrIndex)=Coeff1st;
                                            Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m11),n1(m10),n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                            %Y12(2+ArrIndex)=Coeff1st;
                                            Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m11),n1(m10),n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                            %Y12(3+ArrIndex)=Coeff2nd;
                                            Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(2,n1(m11),n1(m10),n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                            %Y12(4+ArrIndex)=Coeff1st;
                                            Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(3,n1(m11),n1(m10),n1(m9),n1(m8),n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end

         end
    end
end
       
end

.
.

Amin · June 15th, 2021, 3:03 pm

I could not post the functions required for eighth order expansion since they increased the size of maximum characters in a single post. I am posting them here.
.
.

function [IntegralCoeff] = ComputeIntegralCoeffsOrder8()


IntegralCoeff(1:3,1:3,1:3,1:3,1:3,1:3,1:3,1:3)=0;
%IntegralCoeffdt(1:3,1:3,1:3,1:3)=0;
%IntegralCoeffdz(1:3,1:3,1:3,1:3)=0;

IntegralCoeff0(1:3,1:3,1:3,1:3)=0;

IntegralCoeff0(1,1,1,2)=1;
IntegralCoeff0(1,1,1,3)=1;

IntegralCoeff0(1,1,2,2)=1/2;
IntegralCoeff0(1,1,3,2)=1-1/sqrt(3);
IntegralCoeff0(1,1,2,3)=1/sqrt(3);
IntegralCoeff0(1,1,3,3)=1/2;

IntegralCoeff0(1,2,2,2)=1/6;
IntegralCoeff0(1,3,2,2)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5));
IntegralCoeff0(1,2,3,2)=1/sqrt(3)/2*(1-1/sqrt(5));
IntegralCoeff0(1,3,3,2)=1/2*(1-sqrt(2)/2);
IntegralCoeff0(1,2,2,3)=1/2*1/sqrt(5);
IntegralCoeff0(1,3,2,3)=(1-1/sqrt(3))*1/sqrt(2)*1/2;
IntegralCoeff0(1,2,3,3)=1/sqrt(3)/sqrt(2)*1/(2);
IntegralCoeff0(1,3,3,3)=1/6;

IntegralCoeff0(2,2,2,2)=1/24;
IntegralCoeff0(2,2,3,2)=1/2*1/sqrt(5)*1/3*(1-1/sqrt(7));
IntegralCoeff0(2,3,2,2)=1/sqrt(3)*1/2*(1-1/sqrt(5))*1/3*(1-1/sqrt(7));
IntegralCoeff0(3,2,2,2)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5))*1/3*(1-1/sqrt(7));
IntegralCoeff0(3,3,2,2)=1/2*(1-sqrt(2)/2)*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(2,3,3,2)=1/sqrt(3)*1/sqrt(2)*1/2*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(3,2,3,2)=(1-1/sqrt(3))*1/sqrt(2)*1/2*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(3,3,3,2)=1/6*(1-sqrt(3)/sqrt(5));

IntegralCoeff0(2,2,2,3)=1/6*1/sqrt(7);
IntegralCoeff0(2,2,3,3)=1/2*1/sqrt(5)*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(2,3,2,3)=1/sqrt(3)*1/2*(1-1/sqrt(5))*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(3,2,2,3)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5))*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(3,3,2,3)=1/2*(1-sqrt(2)/2)*1/sqrt(3)*1/sqrt(5);
IntegralCoeff0(2,3,3,3)=1/sqrt(3)*1/sqrt(2)*1/sqrt(4)*1/sqrt(3)*(1/sqrt(5));
IntegralCoeff0(3,2,3,3)=(1-1/sqrt(3))*1/sqrt(2)*1/2*1/sqrt(3)*1/sqrt(5);
IntegralCoeff0(3,3,3,3)=1/24;


%Can also be calculated with algorithm below.
%here IntegralCoeffdt indicates the coefficients of a dt-integral.
%This dt-integral means that you are calculating the Ito-hermite expansion
%of  Integral X(t)^alpha dt with X(t) dynamics given by the SDE 
%here IntegralCoeffdz indicates the coefficients of a dz-integral.
%This dz-integral means that you are calculating the Ito-hermite expansion
%of  Integral X(t)^alpha dz(t) with X(t) dynamics given by the SDE 

%IntegralCoeff is is associated with expansion of X(t)^alpha.
%IntegralCoeff below is not returned and IntegralCoeff0 manually calculated
%above is returned but both are the same.

l0(1:2)=1;

for m8=1:2
    l0(1)=1;
    l0(2)=1;
    
    %IntegralCoeff4(m4,1,1,1)=1;
    %IntegralCoeff4(m4,1,1,1)=1;
    %1 is added to m8 since index 1 stands for zero, 2 for one and three
    %for two.
    IntegralCoeff(1,1,1,1,1,1,1,m8+1)=1;
    for m7=1:2
        l0(1)=1;
        l0(2)=1;
        l0(m8)=l0(m8)+1;
        
        if(m7==1)
            IntegralCoeff(1,1,1,1,1,1,m8+1,m7+1)=1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
        end
        if(m7==2)
            IntegralCoeff(1,1,1,1,1,1,m8+1,m7+1)= 1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
        end
        for m6=1:2
            l0(1)=1;
            l0(2)=1;
            l0(m8)=l0(m8)+1;
            l0(m7)=l0(m7)+1;
            
            if(m6==1)
                IntegralCoeff(1,1,1,1,1,m8+1,m7+1,m6+1)=IntegralCoeff(1,1,1,1,1,1,m8+1,m7+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
            end
            if(m6==2)
                IntegralCoeff(1,1,1,1,1,m8+1,m7+1,m6+1)=IntegralCoeff(1,1,1,1,1,1,m8+1,m7+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
            end

            for m5=1:2
                l0(1)=1;
                l0(2)=1;
                l0(m8)=l0(m8)+1;
                l0(m7)=l0(m7)+1;
                l0(m6)=l0(m6)+1;
                if(m5==1)
                    IntegralCoeff(1,1,1,1,m8+1,m7+1,m6+1,m5+1)=IntegralCoeff(1,1,1,1,1,m8+1,m7+1,m6+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                end
                if(m5==2)
                    IntegralCoeff(1,1,1,1,m8+1,m7+1,m6+1,m5+1)=IntegralCoeff(1,1,1,1,1,m8+1,m7+1,m6+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                end
                for m4=1:2
                    l0(1)=1;
                    l0(2)=1;
                    l0(m8)=l0(m8)+1;
                    l0(m7)=l0(m7)+1;
                    l0(m6)=l0(m6)+1;
                    l0(m5)=l0(m5)+1;
                    if(m4==1)
                        IntegralCoeff(1,1,1,m8+1,m7+1,m6+1,m5+1,m4+1)=IntegralCoeff(1,1,1,1,m8+1,m7+1,m6+1,m5+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                    end
                    if(m4==2)
                        IntegralCoeff(1,1,1,m8+1,m7+1,m6+1,m5+1,m4+1)=IntegralCoeff(1,1,1,1,m8+1,m7+1,m6+1,m5+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                    end
                    for m3=1:2
                        l0(1)=1;
                        l0(2)=1;
                        l0(m8)=l0(m8)+1;
                        l0(m7)=l0(m7)+1;
                        l0(m6)=l0(m6)+1;
                        l0(m5)=l0(m5)+1;
                        l0(m4)=l0(m4)+1;
                        if(m3==1)
                            IntegralCoeff(1,1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1)=IntegralCoeff(1,1,1,m8+1,m7+1,m6+1,m5+1,m4+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                        end
                        if(m3==2)
                            IntegralCoeff(1,1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1)=IntegralCoeff(1,1,1,m8+1,m7+1,m6+1,m5+1,m4+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                        end
                        for m2=1:2
                            l0(1)=1;
                            l0(2)=1;
                            l0(m8)=l0(m8)+1;
                            l0(m7)=l0(m7)+1;
                            l0(m6)=l0(m6)+1;
                            l0(m5)=l0(m5)+1;
                            l0(m4)=l0(m4)+1;
                            l0(m3)=l0(m3)+1;
                            if(m2==1)
                                IntegralCoeff(1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1)=IntegralCoeff(1,1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                            end
                            if(m2==2)
                                IntegralCoeff(1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1)=IntegralCoeff(1,1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                            end
                            for m1=1:2
                                l0(1)=1;
                                l0(2)=1;
                                l0(m8)=l0(m8)+1;
                                l0(m7)=l0(m7)+1;
                                l0(m6)=l0(m6)+1;
                                l0(m5)=l0(m5)+1;
                                l0(m4)=l0(m4)+1;
                                l0(m3)=l0(m3)+1;
                                l0(m2)=l0(m2)+1;
                                if(m1==1)
                                    IntegralCoeff(m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                                end
                                if(m1==2)
                                    IntegralCoeff(m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(1,m8+1,m7+1,m6+1,m5+1,m4+1,m3+1,m2+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                                end
                            end
                        end
                    end
                end
            end

        end
    end
end
end

.
,
.

function [Y] = ItoTaylorCoeffsNewOrder8(alpha,beta1,beta2,gamma)


%In the coefficient calculation program which calculates Y(l1,l2,l3,l4), 
%I have used four levels of looping each for relevant expansion order. 
%The first loop takes four values and second loop takes 16 values and 
%third loop takes 64 values and so on. And then each coefficient 
%term can be individually calculated while carefully accounting 
%for path dependence. 
%So for example in a nested loop structure 
%m1= 1:mDim
% 	m2=1:mDim
%		m3=1:mDim
%    l(m1)=l(m1)+1;
%    l(m2)=l(m2)+1;
%    l(m3)=l(m3)+1;

%in the above looping loop takes values from one to four with one 
%indicating the first drift term, two indicating the second drift term 
%and three indicating quadratic variation term and 
%four indicating the volatility term. And with this looping structure 
%we can So in the above looping m1=1 would mean that all terms are 
%descendent of first drift term and m2=4 would mean that all terms are 
%descendent of first drift term on first expansion order and descendent 
%of volatility term on the second order and so we can keep track of path 
%dependence perfectly. 
%And on each level, we individually calculate the descendent terms. While 
%keeping track of path dependence and calculating the coefficients with 
%careful path dependence consideration, we update the appropriate element 
%in our polynomial like expansion coefficient array 

%explaining the part of code
%m1= 1:mDim
%	m2=1:mDim
%		m3=1:mDim
%    l(m1)=l(m1)+1;
%    l(m2)=l(m2)+1;
%    l(m3)=l(m3)+1;
%Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);

%Here l(1) denotes l1 but written as l(1) so it can be conveniently 
%updated with the loop variable when the loop variable takes value one 
%indicating first drift term . And l(2) could be conveniently updated when 
%the loop variable takes value equal to two indicating second 
%drift term and so on.
%Here is the part of code snippet for that

%for m1=1:mDim
%    l(1)=1;
%    l(2)=1;
%    l(3)=1;
%    l(4)=1;
%    l(m1)=l(m1)+1;
%CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
%    CoeffDX2 = CoeffDX1 - 1;
%    ArrIndex0=m1;
%    ArrIndex=(m1-1)*mDim;
%    Coeff1st=Y1(ArrIndex0)*CoeffDX1;
%    Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%    Y2(1+ArrIndex)=Coeff1st;
%    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
%    Y2(2+ArrIndex)=Coeff1st;
%    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
%    Y2(3+ArrIndex)=Coeff2nd;
%    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,2,n1(m1));
%    Y2(4+ArrIndex)=Coeff1st;
%    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,3,n1(m1));

 
%The first four lines update the array indices according to the parent term. 
%And then CoeffDX1 and CoeffDX2 are calculated according to algebraic exponents on parent terms.    
%ArrIndex0=m1; calculates the array index of the parent term 
%And ArrIndex=(m1-1)*mDim; calculates the array index of the descendent terms
%And coefficient of the drift and volatility descendent terms is 
%calculated by multiplying the coefficient of the parent term by 
%Coeff1st=Y1(ArrIndex0)*CoeffDX1;
%And coefficient of the quadratic variation descendent terms is 
%calculated by multiplying the coefficient of the parent term by 
%Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%And then each of the four descendent terms are updated with Coeff1st 
%if they are drift or volatility descendent terms or Coeff2nd if 
%they are quadratic variation descendent terms.
%Here Y1 indicates the temporary coefficient array with parent terms on 
%first level. And Y2 denotes temporary coefficient array with parent terms 
%on second level and Y3 indicates temporary coefficient array with parent terms on third level.

[IntegralCoeff] = ComputeIntegralCoeffsOrder8();

n1(1)=2;
n1(2)=2;
n1(3)=2;
n1(4)=3;
%n1(1), n1(2), n1(3) are drift and quadratic variation term variables 
%and take a value equal to 2 indicating a dt integral.
%n1(4) is volatility term variable and indicates a dz-integral by taking a
%value of 3. 

mDim=4; % four descendent terms in each expansion
Y(1:9,1:9,1:9,1:9)=0;
Y1(1:mDim)=0;
Y2(1:mDim*mDim)=0;
Y3(1:mDim*mDim*mDim)=0;
Y4(1:mDim*mDim*mDim*mDim)=0;
Y5(1:mDim*mDim*mDim*mDim*mDim)=0;
Y6(1:mDim*mDim*mDim*mDim*mDim*mDim)=0;
Y7(1:mDim*mDim*mDim*mDim*mDim*mDim*mDim)=0;


%First Ito-hermite expansion level starts here. No loops but four
%descendent terms.
l(1)=1;
l(2)=1;
l(3)=1;
l(4)=1;
CoeffDX1 = alpha;
CoeffDX2 = CoeffDX1 - 1;
Coeff1st=CoeffDX1;
Coeff2nd=.5*CoeffDX1*CoeffDX2;

Y1(1)=Coeff1st;
Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,2);
Y1(2)=Coeff1st;
Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,2);
Y1(3)=Coeff2nd;
Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,1,1,1,2);
Y1(4)=Coeff1st;
Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,1,1,1,3);

%Second Ito-hermite expansion level starts. It has a loop over four parent
%terms and there are four descendent terms for each parent term.
%The coefficient terms are then lumped in a polynomial-like expansion
%array of coefficients.
for m1=1:mDim
    l(1)=1;
    l(2)=1;
    l(3)=1;
    l(4)=1;
    l(m1)=l(m1)+1;
    CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
    CoeffDX2 = CoeffDX1 - 1;
    ArrIndex0=m1;
    ArrIndex=(m1-1)*mDim;
    Coeff1st=Y1(ArrIndex0)*CoeffDX1;
    Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
    Y2(1+ArrIndex)=Coeff1st;
    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,2,n1(m1));
    Y2(2+ArrIndex)=Coeff1st;
    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,1,2,n1(m1));
    Y2(3+ArrIndex)=Coeff2nd;
    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,1,1,2,n1(m1));
    Y2(4+ArrIndex)=Coeff1st;
    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,1,1,3,n1(m1));
    %Third Ito-hermite expansion level starts and it is a nested loop with
    %a total of sixteen parents and each parent takes four descendent
    %terms.
    for m2=1:mDim
        l(1)=1;
        l(2)=1;
        l(3)=1;
        l(4)=1;
        l(m1)=l(m1)+1;
        l(m2)=l(m2)+1;
        CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
        CoeffDX2=CoeffDX1-1;            
        ArrIndex0=(m1-1)*mDim+m2;
        ArrIndex=((m1-1)*mDim+(m2-1))*mDim;
        Coeff1st=Y2(ArrIndex0)*CoeffDX1;
        Coeff2nd=Y2(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
        Y3(1+ArrIndex)=Coeff1st;
        Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,2,n1(m2),n1(m1));
        Y3(2+ArrIndex)=Coeff1st;
        Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,1,2,n1(m2),n1(m1));
        Y3(3+ArrIndex)=Coeff2nd;
        Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,1,2,n1(m2),n1(m1));
        Y3(4+ArrIndex)=Coeff1st;
        Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,1,3,n1(m2),n1(m1));
        %fourht Ito-hermite expansion level starts and it is a triply-nested loop with
        %a total of sixteen parents and each parent takes four descendent
        %terms. We then lump the terms in a relatively sparse polynomial 
        %like expansion coefficient array that has smaller number of
        %non-zero terms.

        for m3=1:mDim
            l(1)=1;
            l(2)=1;
            l(3)=1;
            l(4)=1;
            l(m1)=l(m1)+1;
            l(m2)=l(m2)+1;
            l(m3)=l(m3)+1;
            CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
            CoeffDX2=CoeffDX1-1;
            ArrIndex0=((m1-1)*mDim+(m2-1))*mDim+m3;
            ArrIndex=(((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim;
            Coeff1st=Y3(ArrIndex0)*CoeffDX1;
            Coeff2nd=Y3(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
            Y4(1+ArrIndex)=Coeff1st;
            Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,2,n1(m3),n1(m2),n1(m1));
            Y4(2+ArrIndex)=Coeff1st;
            Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,1,2,n1(m3),n1(m2),n1(m1));
            Y4(3+ArrIndex)=Coeff2nd;
            Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,1,2,n1(m3),n1(m2),n1(m1));
            Y4(4+ArrIndex)=Coeff1st;
            Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,1,3,n1(m3),n1(m2),n1(m1));
            for m4=1:mDim
                l(1)=1;
                l(2)=1;
                l(3)=1;
                l(4)=1;
                l(m1)=l(m1)+1;
                l(m2)=l(m2)+1;
                l(m3)=l(m3)+1;
                l(m4)=l(m4)+1;
                CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                CoeffDX2=CoeffDX1-1;
                
                ArrIndex0=(((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+m4;
                ArrIndex=((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim;
                
                Coeff1st=Y4(ArrIndex0)*CoeffDX1;
                Coeff2nd=Y4(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                Y5(1+ArrIndex)=Coeff1st;
                Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2,n1(m4),n1(m3),n1(m2),n1(m1));
                Y5(2+ArrIndex)=Coeff1st;
                Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2,n1(m4),n1(m3),n1(m2),n1(m1));
                Y5(3+ArrIndex)=Coeff2nd;
                Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,2,n1(m4),n1(m3),n1(m2),n1(m1));
                Y5(4+ArrIndex)=Coeff1st;
                Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,3,n1(m4),n1(m3),n1(m2),n1(m1));
                for m5=1:mDim
                    l(1)=1;
                    l(2)=1;
                    l(3)=1;
                    l(4)=1;
                    l(m1)=l(m1)+1;
                    l(m2)=l(m2)+1;
                    l(m3)=l(m3)+1;
                    l(m4)=l(m4)+1;
                    l(m5)=l(m5)+1;
                    CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                    CoeffDX2=CoeffDX1-1;
                
                    ArrIndex0=((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+m5;
                    ArrIndex=(((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim;
                
                    Coeff1st=Y5(ArrIndex0)*CoeffDX1;
                    Coeff2nd=Y5(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                    Y6(1+ArrIndex)=Coeff1st;
                    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                    Y6(2+ArrIndex)=Coeff1st;
                    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                    Y6(3+ArrIndex)=Coeff2nd;
                    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,2,n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                    Y6(4+ArrIndex)=Coeff1st;
                    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,3,n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                    for m6=1:mDim
                        l(1)=1;
                        l(2)=1;
                        l(3)=1;
                        l(4)=1;
                        l(m1)=l(m1)+1;
                        l(m2)=l(m2)+1;
                        l(m3)=l(m3)+1;
                        l(m4)=l(m4)+1;
                        l(m5)=l(m5)+1;
                        l(m6)=l(m6)+1;
                        CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                            - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                        CoeffDX2=CoeffDX1-1;
                
                        ArrIndex0=(((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+m6;
                        ArrIndex=((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim;
                
                        Coeff1st=Y6(ArrIndex0)*CoeffDX1;
                        Coeff2nd=Y6(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                        Y7(1+ArrIndex)=Coeff1st;
                        Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,2,n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                        Y7(2+ArrIndex)=Coeff1st;
                        Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,2,n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                        Y7(3+ArrIndex)=Coeff2nd;
                        Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,2,n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                        Y7(4+ArrIndex)=Coeff1st;
                        Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,3,n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                        for m7=1:mDim
                            l(1)=1;
                            l(2)=1;
                            l(3)=1;
                            l(4)=1;
                            l(m1)=l(m1)+1;
                            l(m2)=l(m2)+1;
                            l(m3)=l(m3)+1;
                            l(m4)=l(m4)+1;
                            l(m5)=l(m5)+1;
                            l(m6)=l(m6)+1;
                            l(m7)=l(m7)+1;
                            CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                                - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
                            CoeffDX2=CoeffDX1-1;
                
                            ArrIndex0=((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+m7;
                            ArrIndex=(((((((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim+(m4-1))*mDim+(m5-1))*mDim+(m6-1))*mDim+(m7-1))*mDim;
                
                            Coeff1st=Y7(ArrIndex0)*CoeffDX1;
                            Coeff2nd=Y7(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
                            %Y8(1+ArrIndex)=Coeff1st;
                            Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                            %Y8(2+ArrIndex)=Coeff1st;
                            Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                            %Y8(3+ArrIndex)=Coeff2nd;
                            Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(2,n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                            %Y8(4+ArrIndex)=Coeff1st;
                            Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(3,n1(m7),n1(m6),n1(m5),n1(m4),n1(m3),n1(m2),n1(m1));
                        end
                    end
                end
            end

         end
    end
end
        
        
        

end

.
.
Here is the output and graph you will see when you run the first main program.

Elapsed time is 0.195409 seconds.
TrueMean =
0.106718372568390
MCMean =
0.106777299997107 - 0.000000070142151i
MCVar =
0.009680496098567 + 0.000000008389328i
Original process average from our simulation
ItoHermiteMeanBes =
0.108129928675028
ItoHermiteMeanOrig =
0.106806958066734
ItoHermiteVarBes =
0.010086160986542
ItoHermiteVarOrig =
0.009977942913806
true Mean only applicble to standard SV mean reverting type models otherwise disregard
IndexMax =
398
red line is density of SDE from Ito-Hermite method, green is monte carlo.

Amin · June 15th, 2021, 3:24 pm

Explaining again, the above program makes a single step high order expansion of a single dimensional SDE and compares its density with short stepping monte carlo simulation density of the same SDE.

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