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