Pol

void AssetBasket::DoSimulation(){	Utilities U;	// reset the ratio matrix and set up for importance sampling	if (useImportanceSampling) {		for (int iP=1; iP<=numberOfPaths; iP++) {			likelihoodRatios(iP,1)=1.;		}		dynamicISState = initialISState;	}	// generate the cholesky factorisation of the "covariance matrix"	// ready for the multifactor pricing	// generation the correlated basis	for (int i=1; i<=numberOfUnderliers; i++) {		GenericProcess* b = assetBasket[i];		double v = b->getVol();		correlationMatrix(i,i)=v*v;	}	basisMatrix = correlationMatrix.Cholesky();		#pragma omp parallel for	for (int lPath=1; lPath<=numberOfPaths; lPath++) {		// set initial path values		for (int lFactor=1; lFactor<=numberOfUnderliers; lFactor++) {			assetPaths[lFactor]->operator ()(lPath,1)=assetBasket[lFactor]->getStartValue();				}						for (int lTime=2; lTime<=numberOfTimesteps+1; lTime++) {			// draw normal deviates			Matrix normalDeviates = dz();									// importance sampling goes here			Matrix isDeviates = normalDeviates;			if (useImportanceSampling) {				Matrix& ISD = getISDynamicState();				Matrix isMeans = ISD.Column(1);				Matrix isVariances = ISD.Column(2);				Matrix isSTDs = isVariances; 				// convert to standard deviations in order to				// multiply with the dz vector later				for(int i=1; i<=isVariances.Rows(); i++) {					isSTDs(i,1) = sqrt(isSTDs(i,1));				}				isSTDs = isSTDs.ColumnToDiagonal();				isVariances=isVariances.ColumnToDiagonal();								// apply the importance sampling and calculate				// the next likelihood ratio				// apply the custom importance sampling				for (int i=1; i<=isMeans.Rows(); i++) {					isMeans(i,1)+=ISD(i,3)*assetPaths[i]->operator ()(lPath,lTime-1) / ISD(1,4);				}				// (optimised for diagonal isSTDs)				for (int i=1; i<=isSTDs.Rows(); i++) {					isDeviates(i,1)=(isSTDs(i,i)*isDeviates(i,1))+isMeans(i,1);				}								likelihoodRatios(lPath,lTime)=likelihoodRatios(lPath,lTime-1)*U.mutliLRatioDiagonal(isDeviates,isMeans,isVariances);			} 			// advance the paths			for (int lFactor=1; lFactor<=numberOfUnderliers; lFactor++) {					double dW=0;				for (int i=1; i<=lFactor; i++) {					dW +=  basisMatrix(lFactor,i) * isDeviates(i,1);				}				double pV = assetPaths[lFactor]->operator ()(lPath,lTime-1);				double nV = assetBasket[lFactor]->Advance(pV,dW,dt);				assetPaths[lFactor]->operator ()(lPath,lTime) = nV;			}		}	}}