Merge branch 'multiAssets'

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

choldc.h

@@ -0,0 +1,50 @@
+#pragma once
+
+//	Choleski decomposition, inspired by Numerical Recipes
+
+#include "matrix.h"
+
+template <class T>
+void choldc(const matrix<T>& in, matrix<T>& out)
+{
+	int n = in.rows();
+	T sum;
+
+	fill(out.begin(), out.end(), T(0.0));
+
+	for(int i=0; i<n; ++i)
+	{
+		auto* ai = in[i];
+		auto* pi = out[i];
+		for(int j=0; j<=i; ++j)
+		{
+			auto* aj = in[j];
+			auto* pj = out[j];
+			sum = ai[j];
+			for(int k=0; k<j; ++k)
+			{
+				sum -= pi[k] * pj[k];
+			}
+			if(i == j)
+			{      
+                if(sum < - 1.0e-15)
+                {
+                    throw runtime_error("choldc : matrix not positive definite");
+				}
+                if(sum < 1.0e-15) sum = 0.0;
+				pi[i] = sqrt(sum);
+			}
+			else
+			{
+				if(fabs(pj[j]) < 1.0e-15)
+                {
+                    pi[j] = 0.0;
+                }
+                else
+                {
+                    pi[j] = sum/pj[j];
+                }
+			}
+		}
+	}
+}

main.h

@@ -21,6 +21,7 @@ As long as this comment is preserved at the top of the file
 #include "mcBase.h"
 #include "mcMdl.h"
 #include "mcPrd.h"
+#include "mcPrdMulti.h"
 #include "mrg32k3a.h"
 #include "sobol.h"
 #include <numeric>

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));
     }
 };
 
@@ -116,13 +122,19 @@ inline void initializePath(Scenario<T>& path)
 template <class T>
 class Product
 {
+    inline static const vector<string> defaultAssetNames = { "spot" };
+
 public:
 
     //  Access to the product timeline
     //      along with the sample definitions (defline)
     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 const size_t numAssets() const { return 1; }
+    virtual const vector<string>& assetNames() const { return defaultAssetNames; }
+
     //  Labels of all payoffs in the product
     virtual const vector<string>& payoffLabels() const = 0;
 
@@ -145,8 +157,14 @@ class Product
 template <class T>
 class Model
 {
+    inline static const vector<string> defaultAssetNames = { "spot" };
+
 public:
 
+    //  Number and names of underlying assets, default = 1 and "spot"
+    virtual const size_t numAssets() const { return 1; }
+    virtual const vector<string>& assetNames() const { return defaultAssetNames; }
+
     //  Initialize with product timeline
     virtual void allocate(
         const vector<Time>&         prdTimeline, 
@@ -230,6 +248,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 +270,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 +317,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 +433,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 +570,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 +779,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 +862,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();
 

mcMdl.h

@@ -2,3 +2,4 @@
 
 #include "mcMdlBS.h"
 #include "mcMdlDupire.h"
+#include "mcMdlMultiDisplaced.h"

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: