upgrade platform for multi asset

AADExpr.h

@@ -1107,5 +1107,29 @@ class Number : public Expression<Number>
         *this = *this / e;
         return *this;
     }
+
+    Number& operator+=(const double& e)
+    {
+        *this = *this + e;
+        return *this;
+    }
+
+    Number& operator*=(const double& e)
+    {
+        *this = *this * e;
+        return *this;
+    }
+
+    Number& operator-=(const double& e)
+    {
+        *this = *this - e;
+        return *this;
+    }
+
+    Number& operator/=(const double& e)
+    {
+        *this = *this / e;
+        return *this;
+    }
 };
 

mcBase.h

@@ -47,9 +47,6 @@ struct SampleDef
     //  need numeraire?
     bool            numeraire = true;
 
-    vector<Time>    forwardMats;
-    vector<Time>    discountMats;
-
     struct RateDef
     {
         Time    start;
@@ -60,7 +57,11 @@ struct SampleDef
             start(s), end(e), curve(c) {};
     };
 
+    vector<Time>    discountMats;
     vector<RateDef> liborDefs;
+
+    //  multi-asset: forwardMats[a] = maturities for asset a
+    vector<vector<Time>>        forwardMats;
 };
 
 //  Sample = simulated value
@@ -69,25 +70,30 @@ template <class T>
 struct Sample
 {
     T           numeraire;
-    vector<T>   forwards;
     vector<T>   discounts;
     vector<T>   libors;
 
+    //  multi-asset: forwardMats[a][t] = forward for asset a, maturity t
+    vector<vector<T>>   forwards;
+
     //  Allocate given SampleDef
     void allocate(const SampleDef& data)
     {
-        forwards.resize(data.forwardMats.size());
         discounts.resize(data.discountMats.size());
         libors.resize(data.liborDefs.size());
+
+        forwards.resize(data.forwardMats.size());
+        for (size_t a = 0; a < forwards.size(); ++a) forwards[a].resize(data.forwardMats[a].size());
     }
 
     //  Initialize defaults
     void initialize()
     {
         numeraire = T(1.0);
-		fill(forwards.begin(), forwards.end(), T(100.0));
         fill(discounts.begin(), discounts.end(), T(1.0));
         fill(libors.begin(), libors.end(), T(0.0));
+
+		for (auto& forward: forwards) fill(forward.begin(), forward.end(), T(100.0));
     }
 };
 
@@ -123,6 +129,10 @@ class Product
     virtual const vector<Time>& timeline() const = 0;
     virtual const vector<SampleDef>& defline() const = 0;
 
+    //  Number and names of underlying assets, default = 1 and "spot"
+    virtual size_t numAssets() const { return 1; }
+    virtual vector<string> assetNames() const { return { "spot" }; }
+
     //  Labels of all payoffs in the product
     virtual const vector<string>& payoffLabels() const = 0;
 
@@ -147,6 +157,10 @@ class Model
 {
 public:
 
+    //  Number and names of underlying assets, default = 1 and "spot"
+    virtual size_t numAssets() const { return 1; }
+    virtual vector<string> assetNames() const { return { "spot" }; }
+
     //  Initialize with product timeline
     virtual void allocate(
         const vector<Time>&         prdTimeline, 
@@ -230,6 +244,17 @@ class RNG
 //  Template algorithms
 //  ===================
 
+//  Check compatibility of model and product
+//  At the moment, only check that assets are the samein both cases
+//  May be easily extended in the future
+template <class T>
+inline bool checkCompatiblity(
+    const Product<T>&      prd,
+    const Model<T>&        mdl)
+{
+    return prd.assetNames() == mdl.assetNames();
+}
+
 //  Serial valuation, chapter 6
 
 //	MC simulator: free function that conducts simulations 
@@ -241,6 +266,8 @@ inline vector<vector<double>> mcSimul(
     const RNG&                  rng,			            
     const size_t                nPath)                      
 {
+    if (!checkCompatiblity(prd, mdl)) throw runtime_error("Model and product are not compatible");
+
     //  Work with copies of the model and RNG
     //      which are modified when we set up the simulation
     //  Copies are OK at high level
@@ -286,6 +313,8 @@ inline vector<vector<double>> mcParallelSimul(
     const RNG&                  rng,
     const size_t                nPath)
 {
+    if (!checkCompatiblity(prd, mdl)) throw runtime_error("Model and product are not compatible");
+
     auto cMdl = mdl.clone();
 
     const size_t nPay = prd.payoffLabels().size();
@@ -400,6 +429,8 @@ mcSimulAAD(
     const size_t            nPath,
     const F&                aggFun = defaultAggregator)
 {
+    if (!checkCompatiblity(prd, mdl)) throw runtime_error("Model and product are not compatible");
+
     //  Work with copies of the model and RNG
     //      which are modified when we set up the simulation
     //  Copies are OK at high level
@@ -535,6 +566,8 @@ mcParallelSimulAAD(
     const size_t            nPath,
     const F&                aggFun = defaultAggregator)
 {
+    if (!checkCompatiblity(prd, mdl)) throw runtime_error("Model and product are not compatible");
+
     const size_t nPay = prd.payoffLabels().size();
     const size_t nParam = mdl.numParams();
 
@@ -742,6 +775,8 @@ mcSimulAADMulti(
 	const RNG&              rng,
 	const size_t            nPath)
 {
+    if (!checkCompatiblity(prd, mdl)) throw runtime_error("Model and product are not compatible");
+
 	auto cMdl = mdl.clone();
 	auto cRng = rng.clone();
 
@@ -823,6 +858,8 @@ mcParallelSimulAADMulti(
 	const RNG& rng,
 	const size_t            nPath)
 {
+    if (!checkCompatiblity(prd, mdl)) throw runtime_error("Model and product are not compatible");
+
 	const size_t nPay = prd.payoffLabels().size();
 	const size_t nParam = mdl.numParams();
 

mcMdlBS.h

@@ -258,11 +258,11 @@ class BlackScholes : public Model<T>
 			}
 
 			//  Forward factors
-			const size_t pFF = defline[i].forwardMats.size();
+			const size_t pFF = defline[i].forwardMats.front().size();
 			for (size_t j = 0; j < pFF; ++j)
 			{
 				myForwardFactors[i][j] =
-					exp(mu * (defline[i].forwardMats[j] - productTimeline[i]));
+					exp(mu * (defline[i].forwardMats.front()[j] - productTimeline[i]));
 			}
 
 			//  Libors
@@ -299,7 +299,7 @@ class BlackScholes : public Model<T>
         }
 
         transform(myForwardFactors[idx].begin(), myForwardFactors[idx].end(), 
-            scen.forwards.begin(), 
+            scen.forwards.front().begin(), 
             [&spot](const T& ff)
             {
                 return spot * ff;

mcMdlDupire.h

@@ -227,7 +227,7 @@ class Dupire : public Model<T>
     //  Helper function, fills a sample given the spot
     inline static void fillScen(const T& spot, Sample<T>& scen)
     {
-        fill(scen.forwards.begin(), scen.forwards.end(), spot);
+        fill(scen.forwards.front().begin(), scen.forwards.front().end(), spot);
     }
 
 public:

mcMdlMultiDisplaced.h

@@ -8,6 +8,7 @@ check
 check
 
 - upgrade engine to multiple, named, assets, update all models and products, check adequation of model and product in engine
+check
 
 - reg tests
 
@@ -58,6 +59,13 @@ class MultiDisplaced : public Model<T>
 public:
 	//	
 
+    enum Dynamics
+    {
+        Normal = 0,
+        Surnormal = 1,
+        Subnormal = 2
+    };
+
     //  Model parameters
 
 	//	Number of assets and reference prices
@@ -90,6 +98,7 @@ class MultiDisplaced : public Model<T>
 	//	Derived
 	vector<T>			myAlphas;
 	vector<T>			myBetas;
+    vector<Dynamics>    myDynamics;
 
 	//	Finally, correl
 	//	Lower diagonal matrix
@@ -101,9 +110,6 @@ class MultiDisplaced : public Model<T>
 	//	Cholesky coefs, also lower diagonal
 	matrix<T>			myChol;
 
-	//	Temporary
-	vector<T>			myCorrelatedGaussians;
-
     //  Similuation timeline = today + div dates + event dates
     vector<Time>		myTimeline;
     //  true (1) if the time step is an event date
@@ -272,7 +278,7 @@ class MultiDisplaced : public Model<T>
 
     //  Read access to parameters
 
-	const size_t numAssets() const 
+	size_t numAssets() const override
 	{
 		return myNumAssets;
 	}
@@ -343,12 +349,11 @@ class MultiDisplaced : public Model<T>
     {
 		myAlphas.resize(myNumAssets);
 		myBetas.resize(myNumAssets);
+        myDynamics.resize(myNumAssets);
 
 		myUsedCorrel.resize(myNumAssets, myNumAssets);		
 		myChol.resize(myNumAssets, myNumAssets);
 
-		myCorrelatedGaussians.resize(myNumAssets);
-
         //  Simulation timeline = today + div dates + product timeline
 
 		set<Time> times;
@@ -423,11 +428,23 @@ class MultiDisplaced : public Model<T>
 		for (size_t asset=0; asset<myNumAssets; ++asset)
 		{
 			myBetas[asset] = myAtms[asset] + 2 * mySkews[asset];
-			if (fabs(myBetas[asset]) < 1.0e-08)
+			if (fabs(myBetas[asset]) < 1.0e-05)
 			{
-				myBetas[asset] = 1.0e-08;
+                myDynamics[asset] = Normal;
+			    myAlphas[asset] = - 2 * mySpots[asset] * mySkews[asset];	
+				myBetas[asset] = 0.0;
 			}
-			myAlphas[asset] = - 2 * mySpots[asset] / fabs(myBetas[asset]) * mySkews[asset];	
+            else if (myBetas[asset] > 0)
+            {
+                myDynamics[asset] = Surnormal;
+			    myAlphas[asset] = - 2 * mySpots[asset] / myBetas[asset] * mySkews[asset];	
+            }
+            else    //  neg beta
+            {
+                myDynamics[asset] = Subnormal;
+                myBetas[asset] *= -1.0;   //  pos now
+			    myAlphas[asset] = - 2 * mySpots[asset] / myBetas[asset] * mySkews[asset];	            
+            }
 		}
 
 		//	Apply lambda factor to correl
@@ -478,14 +495,6 @@ class MultiDisplaced : public Model<T>
 					exp(-myRate * (defline[i].discountMats[j] - productTimeline[i]));
 			}
 
-			//  Forward factors
-			const size_t pFF = defline[i].forwardMats.size();
-			for (size_t j = 0; j < pFF; ++j)
-			{
-				myForwardFactors[i][j] =
-					exp(mu * (defline[i].forwardMats[j] - productTimeline[i]));
-			}
-
 			//  Libors
 			const size_t pL = defline[i].liborDefs.size();
 			for (size_t j = 0; j < pL; ++j)
@@ -494,6 +503,18 @@ class MultiDisplaced : public Model<T>
 					= defline[i].liborDefs[j].end - defline[i].liborDefs[j].start;
 				myLibors[i][j] = (exp(myRate*dt) - 1.0) / dt;
 			}
+
+			//  Forward factors
+
+            /*
+            const size_t pFF = defline[i].forwardMats.size();
+			for (size_t j = 0; j < pFF; ++j)
+			{
+				myForwardFactors[i][j] =
+					exp(mu * (defline[i].forwardMats[j] - productTimeline[i]));
+			}
+            */
+
 		}   //  loop on event dates
 	}