1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2008 Roland Lichters
5 Copyright (C) 2014 Jose Aparicio
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#ifndef quantlib_latent_model_hpp
22#define quantlib_latent_model_hpp
23
24#include <ql/experimental/math/multidimquadrature.hpp>
25#include <ql/experimental/math/multidimintegrator.hpp>
26#include <ql/math/integrals/trapezoidintegral.hpp>
27#include <ql/math/randomnumbers/randomsequencegenerator.hpp>
28#include <ql/experimental/math/gaussiancopulapolicy.hpp>
29#include <ql/experimental/math/tcopulapolicy.hpp>
30#include <ql/math/randomnumbers/boxmullergaussianrng.hpp>
31#include <ql/experimental/math/polarstudenttrng.hpp>
32#include <ql/handle.hpp>
33#include <ql/quote.hpp>
34#include <vector>
35
36/*! \file latentmodel.hpp
37 \brief Generic multifactor latent variable model.
38*/
39
40namespace QuantLib {
41
42 namespace detail {
43 // havent figured out how to do this in-place
44 struct multiplyV {
45 /*! \deprecated Use `auto` or `decltype` instead.
46 Deprecated in version 1.29.
47 */
48 QL_DEPRECATED
49 typedef std::vector<Real> result_type;
50
51 std::vector<Real> operator()(Real d, std::vector<Real> v)
52 {
53 std::transform(first: v.begin(), last: v.end(), result: v.begin(),
54 unary_op: [=](Real x) -> Real { return x * d; });
55 return v;
56 }
57 };
58 }
59
60 //! \name Latent model direct integration facility.
61 //@{
62 /* Things trying to achieve here:
63 1.- Unify the two branches of integrators in the library, they do not
64 hang from a common base class and here a common ptr for the
65 factory is needed.
66 2.- Have a common signature for the integration call.
67 3.- Factory construction so integrable latent models can choose the
68 integration algorithm separately.
69 */
70 class LMIntegration {
71 public:
72 // Interface with actual integrators:
73 // integral of a scalar function
74 virtual Real integrate(const ext::function<Real (
75 const std::vector<Real>& arg)>& f) const = 0;
76 // integral of a vector function
77 /* I had to use a different name, since the compiler does not
78 recognise the overload; MSVC sees the argument as
79 ext::function<Signature> in both cases....
80 I could do the as with the quadratures and have this as a template
81 function and spez for the vector case but I prefer to understand
82 why the overload fails....
83 FIX ME
84 */
85 virtual std::vector<Real> integrateV(
86 const ext::function<std::vector<Real> (
87 const std::vector<Real>& arg)>& f) const {
88 QL_FAIL("No vector integration provided");
89 }
90 virtual ~LMIntegration() = default;
91 };
92
93 //CRTP-ish for joining the integrations, class above to have the factory
94 template <class I_T>
95 class IntegrationBase :
96 public I_T, public LMIntegration {// diamond on 'integrate'
97 // this class template always to be fully specialized:
98 private:
99 IntegrationBase() = default;
100 };
101 //@}
102
103 // gcc reports value collision with heston engine (?!) thats why the name
104 namespace LatentModelIntegrationType {
105 typedef
106 enum LatentModelIntegrationType {
107 #ifndef QL_PATCH_SOLARIS
108 GaussianQuadrature,
109 #endif
110 Trapezoid
111 // etc....
112 } LatentModelIntegrationType;
113 }
114
115 #ifndef QL_PATCH_SOLARIS
116
117 /* class template specializations. I havent use CRTP type cast directly
118 because the signature of the integrators is different, grid integration
119 needs the domain. */
120 template<> class IntegrationBase<GaussianQuadMultidimIntegrator> :
121 public GaussianQuadMultidimIntegrator, public LMIntegration {
122 public:
123 IntegrationBase(Size dimension, Size order)
124 : GaussianQuadMultidimIntegrator(dimension, order) {}
125 Real integrate(const ext::function<Real(const std::vector<Real>& arg)>& f) const override {
126 return GaussianQuadMultidimIntegrator::integrate<Real>(f);
127 }
128 std::vector<Real> integrateV(
129 const ext::function<std::vector<Real>(const std::vector<Real>& arg)>& f)
130 const override {
131 return GaussianQuadMultidimIntegrator::integrate<std::vector<Real>>(f);
132 }
133 ~IntegrationBase() override = default;
134 };
135
136 #endif
137
138 template<> class IntegrationBase<MultidimIntegral> :
139 public MultidimIntegral, public LMIntegration {
140 public:
141 IntegrationBase(
142 const std::vector<ext::shared_ptr<Integrator> >& integrators,
143 Real a, Real b)
144 : MultidimIntegral(integrators),
145 a_(integrators.size(),a), b_(integrators.size(),b) {}
146 Real integrate(const ext::function<Real(const std::vector<Real>& arg)>& f) const override {
147 return MultidimIntegral::operator()(f, a: a_, b: b_);
148 }
149 // vector version here....
150 ~IntegrationBase() override = default;
151 const std::vector<Real> a_, b_;
152 };
153
154 // Intended to replace OneFactorCopula
155
156 /*!
157 \brief Generic multifactor latent variable model.\par
158 In this model set up one considers latent (random) variables
159 \f$ Y_i \f$ described by:
160 \f[
161 \begin{array}{ccccc}
162 Y_1 & = & \sum_k M_k a_{1,k} & + \sqrt{1-\sum_k a_{1,k}^2} Z_1 &
163 \sim \Phi_{Y_1}\nonumber \\
164 ... & = & ... & ... & \nonumber \\
165 Y_i & = & \sum_k M_k a_{i,k} & + \sqrt{1-\sum_k a_{i,k}^2} Z_i &
166 \sim \Phi_{Y_i}\nonumber \\
167 ... & = & ... & ... & \nonumber \\
168 Y_N & = & \sum_k M_k a_{N,k} & + \sqrt{1-\sum_k a_{N,k}^2} Z_N &
169 \sim \Phi_{Y_N}
170 \end{array}
171 \f]
172 where the systemic \f$ M_k \f$ and idiosyncratic \f$ Z_i \f$ (this last
173 one known as error term in some contexts) random variables have
174 independent zero-mean unit-variance distributions. A restriction of the
175 model implemented here is that the N idiosyncratic variables all follow
176 the same probability law \f$ \Phi_Z(z)\f$ (but they are still
177 independent random variables) Also the model is normalized
178 so that: \f$-1\leq a_{i,k} \leq 1\f$ (technically the \f$Y_i\f$ are
179 convex linear combinations). The correlation between \f$Y_i\f$ and
180 \f$Y_j\f$ is then \f$\sum_k a_{i,k} a_{j,k}\f$.
181 \f$\Phi_{Y_i}\f$ denotes the cumulative distribution function of
182 \f$Y_i\f$ which in general differs for each latent variable.\par
183 In its single factor set up this model is usually employed in derivative
184 pricing and it is best to use it through integration of the desired
185 statistical properties of the model; in its multifactorial version (with
186 typically around a dozen factors) it is used in the context of portfolio
187 risk metrics; because of the number of variables it is best to opt for a
188 simulation to compute model properties/magnitudes.
189 For this reason this class template provides a random factor sample
190 interface and an integration interface that will be instantiated by
191 derived concrete models as needed. The class is neutral on the
192 integration and random generation algorithms\par
193 The latent variables are typically treated as unobservable magnitudes
194 and they serve to model one or several magnitudes related to them
195 through some function
196 \f[
197 \begin{array}{ccc}
198 F_i(Y_i) & = &
199 F_i(\sum_k M_k a_{i,k} + \sqrt{1-\sum_k a_{i,k}^2} Z_i )\nonumber \\
200 & = & F_i(M_1,..., M_k, ..., M_K, Z_i)
201 \end{array}
202 \f]
203 The transfer function can have a more generic form:
204 \f$F_i(Y_1,....,Y_N)\f$ but here the model is restricted to a one to
205 one relation between the latent variables and the modelled ones. Also
206 it is assumed that \f$F_i(y_i; \tau)\f$ is monotonic in \f$y_i\f$; it
207 can then be inverted and the relation of the cumulative probability of
208 \f$F_i\f$ and \f$Y_i\f$ is simple:
209 \f[
210 \int_{\infty}^b \phi_{F_i} df =
211 \int_{\infty}^{F_i^{-1}(b)} \phi_{Y_i} dy
212 \f]
213 If \f$t\f$ is some value of the functional or modelled variable,
214 \f$y\f$ is mapped to \f$t\f$ such that percentiles match, i.e.
215 \f$F_Y(y)=Q_i(t)\f$ or \f$y=F_Y^{-1}(Q_i(t))\f$.
216 The class provides an integration facility of arbitrary functions
217 dependent on the model states. It also provides random number generation
218 interfaces for usage of the model in monte carlo simulations.\par
219 Now let \f$\Phi_Z(z)\f$ be the cumulated distribution function of (all
220 equal as mentioned) \f$Z_i\f$. For a given realization of \f$M_k\f$,
221 this determines the distribution of \f$y\f$:
222 \f[
223 Prob \,(Y_i < y|M_k) = \Phi_Z \left( \frac{y-\sum_k a_{i,k}\,M_k}
224 {\sqrt{1-\sum_k a_{i,k}^2}}\right)
225 \qquad
226 \mbox{or}
227 \qquad
228 Prob \,(t_i < t|M) = \Phi_Z \left( \frac
229 {F_{Y_{i}}^{-1}(Q_i(t))-\sum_k a_{i,k}\,M_k}
230 {\sqrt{1-\sum_k a_{i,k}^2}}
231 \right)
232 \f]
233 The distribution functions of \f$ M_k, Z_i \f$ are specified in
234 specific copula template classes. The distribution function
235 of \f$ Y_i \f$ is then given by the convolution
236 \f[
237 F_{Y_{i}}(y) = Prob\,(Y_i<y) =
238 \int_{-\infty}^\infty\,\cdots\,\int_{-\infty}^{\infty}\:
239 D_Z(z)\,\prod_k D_{M_{k}}(m_k) \quad
240 \Theta \left(y - \sum_k a_{i,k}m_k -
241 \sqrt{1-\sum_k a_{i,k}^2}\,z\right)\,d\bar{m}\,dz,
242 \qquad
243 \Theta (x) = \left\{
244 \begin{array}{ll}
245 1 & x \geq 0 \\
246 0 & x < 0
247 \end{array}\right.
248 \f]
249 where \f$ D_Z(z) \f$ and \f$ D_M(m) \f$ are the probability
250 densities of \f$ Z\f$ and \f$ M, \f$ respectively.\par
251 This convolution can also be written
252 \f[
253 F_{Y_{i}}(y) = Prob \,(Y_i < y) =
254 \int_{-\infty}^\infty\,\cdots\,\int_{-\infty}^{\infty}
255 D_{M_{k}}(m_k)\,dm_k\:
256 \int_{-\infty}^{g(y,\vec{a},\vec{m})} D_Z(z)\,dz, \qquad
257 g(y,\vec{a},\vec{m}) = \frac{y - \sum_k a_{i,k}m_k}
258 {\sqrt{1-\sum_k a_{i,k}^2}}, \qquad \sum_k a_{i,k}^2 < 1
259 \f]
260 In general, \f$ F_{Y_{i}}(y) \f$ needs to be computed numerically.\par
261 The policy class template separates the copula function (the
262 distributions involved) and the functionality (i.e. what the latent
263 model represents: a default probability, a recovery...). Since the
264 copula methods for the
265 probabilities are to be called repeatedly from an integration or a MC
266 simulation, virtual tables are avoided and template parameter mechnics
267 is preferred.\par
268 There is nothing at this level enforncing the requirement
269 on the factor distributions to be of zero mean and unit variance. Thats
270 the user responsibility and the model fails to behave correctly if it
271 is not the case.\par
272 Derived classes should implement a modelled magnitude (default time,
273 etc) and will provide probability distributions and conditional values.
274 They could also provide functionality for the parameter inversion
275 problem, the (e.g.) time at which the modeled variable first takes a
276 given value. This problem has solution/sense depending on the transfer
277 function \f$F_i(Y_i)\f$ characteristics.
278
279 To make direct integration and simulation time efficient virtual
280 functions have been avoided in accessing methods in the copula policy
281 and in the sampling of the random factors
282 */
283 template <class copulaPolicyImpl>
284 class LatentModel
285 : public virtual Observer , public virtual Observable
286 {//observer if factors as quotes
287 public:
288 void update() override;
289 //! \name Copula interface.
290 //@{
291 typedef copulaPolicyImpl copulaType;
292 /*! Cumulative probability of the \f$ Y_i \f$ modelled latent random
293 variable to take a given value.
294 */
295 Probability cumulativeY(Real val, Size iVariable) const {
296 return copula_.cumulativeY(val, iVariable);
297 }
298 //! Cumulative distribution of Z, the idiosyncratic/error factors.
299 Probability cumulativeZ(Real z) const {
300 return copula_.cumulativeZ(z);
301 }
302 //! Density function of M, the market/systemic factors.
303 Probability density(const std::vector<Real>& m) const {
304 #if defined(QL_EXTRA_SAFETY_CHECKS)
305 QL_REQUIRE(m.size() == nFactors_,
306 "Factor size must match that of model.");
307 #endif
308 return copula_.density(m);
309 }
310 //! Inverse cumulative distribution of the systemic factor iFactor.
311 Real inverseCumulativeDensity(Probability p, Size iFactor) const {
312 return copula_.inverseCumulativeDensity(p, iFactor);
313 }
314 /*! Inverse cumulative value of the i-th random latent variable with a
315 given probability. */
316 Real inverseCumulativeY(Probability p, Size iVariable) const {
317 return copula_.inverseCumulativeY(p, iVariable);
318 }
319 /*! Inverse cumulative value of the idiosyncratic variable with a given
320 probability. */
321 Real inverseCumulativeZ(Probability p) const {
322 return copula_.inverseCumulativeZ(p);
323 }
324 /*! All factor cumulative inversion. Used in integrations and sampling.
325 Inverts all the cumulative random factors probabilities in the
326 model. These are all the systemic factors plus all the idiosyncratic
327 ones, so the size of the inversion is the number of systemic factors
328 plus the number of latent modelled variables*/
329 std::vector<Real> allFactorCumulInverter(const std::vector<Real>& probs) const {
330 return copula_.allFactorCumulInverter(probs);
331 }
332 //@}
333
334 /*! The value of the latent variable Y_i conditional to
335 (given) a set of values of the factors.
336
337 The passed allFactors vector contains values for all the
338 independent factors in the model (systemic and
339 idiosyncratic, in that order). A full sample is required,
340 i.e. all the idiosyncratic values are expected to be
341 present even if only the relevant one is used.
342 */
343 Real latentVarValue(const std::vector<Real>& allFactors,
344 Size iVar) const
345 {
346 return std::inner_product(first1: factorWeights_[iVar].begin(),
347 // systemic term:
348 last1: factorWeights_[iVar].end(), first2: allFactors.begin(),
349 // idiosyncratic term:
350 init: Real(allFactors[numFactors()+iVar] * idiosyncFctrs_[iVar]));
351 }
352 // \to do write variants of the above, although is the most common case
353
354 const copulaType& copula() const {
355 return copula_;
356 }
357
358
359 // protected:
360 //! \name Latent model random factor number generator facility.
361 //@{
362 /*! Allows generation or random samples of the latent variable.
363
364 Generates samples of all the factors in the latent model according
365 to the given copula as random sequence. The default implementation
366 given uses the inversion in the copula policy (which must be
367 present).
368 USNG is expected to be a uniform sequence generator in the default
369 implementation.
370 */
371 /*
372 Several (very different) usages make the spez non trivial
373 The final goal is to obtain a sequence generator of the factor
374 samples, several routes are possible depending on the algorithms:
375
376 1.- URNG -> Sequence Gen -> CopulaInversion
377 e.g.: CopulaInversion(RandomSequenceGenerator<MersenneTwisterRNG>)
378 2.- PseudoRSG ------------> CopulaInversion
379 e.g.: CopulaInversion(SobolRSG)
380 3.- URNG -> SpecificMapping -> Sequence Gen (bypasses the copula
381 for performance)
382 e.g.: RandomSequenceGenerator<BoxMullerGaussianRng<
383 MersenneTwisterRNG> >
384
385 Notice that the order the three algorithms involved (uniform gen,
386 sequence construction, distribution mapping) is not always the same.
387 (in fact there could be some other ways to generate but these are
388 the ones in the library now.)
389 Difficulties arise when wanting to use situation 3.- whith a generic
390 RNG, leaving it unspecified
391
392 Derived classes might specialize (on the copula
393 type) to another type of generator if a more efficient algorithm
394 that the distribution inversion is available; rewritig then the
395 nextSequence method for a particular copula implementation.
396 Some combinations of generators might make no sense, while it
397 could be possible to block template classes corresponding to those
398 cases its not done (yet?) (e.g. a BoxMuller under a TCopula.)
399 Dimensionality coherence (between the generator and the copula)
400 should have been checked by the client code.
401 In multithread usage the sequence generator is expect to be already
402 in position.
403 To sample the latent variable itself users should call
404 LatentModel::latentVarValue with these samples.
405 */
406 // Cant use InverseCumulativeRsg since the inverse there has to return a
407 // real number and here a vector is needed, the function inverted here
408 // is multivalued.
409 template <class USNG,
410 // dummy template parameter to allow for 'full' specialization of
411 // inner class without specialization of the outer.
412 bool = true>
413 class FactorSampler {
414 public:
415 typedef Sample<std::vector<Real> > sample_type;
416 explicit FactorSampler(const copulaType& copula,
417 BigNatural seed = 0)
418 : sequenceGen_(copula.numFactors(), seed), // base case construction
419 x_(std::vector<Real>(copula.numFactors()), 1.0),
420 copula_(copula) { }
421 /*! Returns a sample of the factor set \f$ M_k\,Z_i\f$.
422 This method has the vocation of being specialized at particular
423 types of the copula with a more efficient inversion to generate the
424 random variables modelled (e.g. Box-Muller for a gaussian).
425 Here a default implementation is provided based directly on the
426 inversion of the cumulative distribution from the copula.
427 Care has to be taken in potential specializations that the generator
428 algorithm is compatible with an eventual concurrence of the
429 simulations.
430 */
431 const sample_type& nextSequence() const {
432 typename USNG::sample_type sample =
433 sequenceGen_.nextSequence();
434 x_.value = copula_.allFactorCumulInverter(sample.value);
435 return x_;
436 }
437 private:
438 USNG sequenceGen_;// copy, we might be mutithreaded
439 mutable sample_type x_;
440 // no copies
441 const copulaType& copula_;
442 };
443 //@}
444 protected:
445 /* \todo Move integrator traits like number of quadrature points,
446 integration domain dimensions, etc to the copula through a static
447 member function. Since they depend on the nature of the probability
448 density distribution thats where they belong.
449 This is why theres one factory per copula policy template parameter
450 (even if this is not used...yet)
451 */
452 class IntegrationFactory {
453 public:
454 static ext::shared_ptr<LMIntegration> createLMIntegration(
455 Size dimension,
456 LatentModelIntegrationType::LatentModelIntegrationType type =
457 #ifndef QL_PATCH_SOLARIS
458 LatentModelIntegrationType::GaussianQuadrature)
459 #else
460 LatentModelIntegrationType::Trapezoid)
461 #endif
462 {
463 switch(type) {
464 #ifndef QL_PATCH_SOLARIS
465 case LatentModelIntegrationType::GaussianQuadrature:
466 return
467 ext::make_shared<
468 IntegrationBase<GaussianQuadMultidimIntegrator> >(
469 args&: dimension, args: 25);
470 #endif
471 case LatentModelIntegrationType::Trapezoid:
472 {
473 std::vector<ext::shared_ptr<Integrator> > integrals;
474 for(Size i=0; i<dimension; i++)
475 integrals.push_back(
476 x: ext::make_shared<TrapezoidIntegral<Default> >(
477 args: 1.e-4, args: 20));
478 /* This integration domain is tailored for the T
479 distribution; it is too wide for normals or Ts of high
480 order.
481 \todo This needs to be solved by having the copula to
482 provide the integration traits for any integration
483 algorithm since it is the copula that knows the relevant
484 domain for its density distributions. Also to be able to
485 block integrations which will fail; like a quadrature
486 here in some cases.
487 */
488 return
489 ext::make_shared<IntegrationBase<MultidimIntegral> >
490 (args&: integrals, args: -35., args: 35.);
491 }
492 default:
493 QL_FAIL("Unknown latent model integration type.");
494 }
495 }
496 private:
497 IntegrationFactory() = default;
498 };
499 //@}
500
501
502 public:
503 // model size, number of latent variables modelled
504 Size size() const {return nVariables_;}
505 //! Number of systemic factors.
506 Size numFactors() const {return nFactors_;}
507 //! Number of total free random factors; systemic and idiosyncratic.
508 Size numTotalFactors() const { return nVariables_ + nFactors_; }
509
510 /*! Constructs a LM with an arbitrary number of latent variables
511 and factors given by the dimensions of the passed matrix.
512 @param factorsWeights Ordering is factorWeights_[iVar][iFactor]
513 @param ini Initialization variables. Trait type from the copula
514 policy to allow for static policies (this solution needs to be
515 revised, possibly drop the static policy and create a policy
516 member in LatentModel)
517 */
518 explicit LatentModel(
519 const std::vector<std::vector<Real> >& factorsWeights,
520 const typename copulaType::initTraits& ini =
521 copulaType::initTraits());
522 /*! Constructs a LM with an arbitrary number of latent variables
523 depending only on one random factor but contributing to each latent
524 variable through different weights.
525 @param factorsWeight Ordering is factorWeights_[iVariable]
526 @param ini Initialization variables. Trait type from the copula
527 policy to allow for static policies (this solution needs to be
528 revised, possibly drop the static policy and create a policy
529 member in LatentModel)
530 */
531 explicit LatentModel(const std::vector<Real>& factorsWeight,
532 const typename copulaType::initTraits& ini =
533 copulaType::initTraits());
534 /*! Constructs a LM with an arbitrary number of latent variables
535 depending only on one random factor with the same weight for all
536 latent variables.
537
538 correlSqr is the weight, same for all.
539
540 ini is a trait type from the copula policy, to allow for
541 static policies (this solution needs to be revised,
542 possibly drop the static policy and create a policy member
543 in LatentModel)
544 */
545 explicit LatentModel(Real correlSqr,
546 Size nVariables,
547 const typename copulaType::initTraits& ini = copulaType::initTraits());
548 /*! Constructs a LM with an arbitrary number of latent variables
549 depending only on one random factor with the same weight for all
550 latent variables. The weight is observed and this constructor is
551 intended to be used when the model relates to a market value.
552
553 singleFactorCorrel is the weight/mkt-factor, same for all.
554
555 ini is a trait type from the copula policy, to allow for
556 static policies (this solution needs to be revised,
557 possibly drop the static policy and create a policy member
558 in LatentModel)
559 */
560 explicit LatentModel(const Handle<Quote>& singleFactorCorrel,
561 Size nVariables,
562 const typename copulaType::initTraits& ini =
563 copulaType::initTraits());
564
565 //! Provides values of the factors \f$ a_{i,k} \f$
566 const std::vector<std::vector<Real> >& factorWeights() const {
567 return factorWeights_;
568 }
569 //! Provides values of the normalized idiosyncratic factors \f$ Z_i \f$
570 const std::vector<Real>& idiosyncFctrs() const {return idiosyncFctrs_;}
571
572 //! Latent variable correlations:
573 Real latentVariableCorrel(Size iVar1, Size iVar2) const {
574 // true for any normalized combination
575 Real init = (iVar1 == iVar2 ?
576 idiosyncFctrs_[iVar1] * idiosyncFctrs_[iVar1] : Real(0.));
577 return std::inner_product(first1: factorWeights_[iVar1].begin(),
578 last1: factorWeights_[iVar1].end(), first2: factorWeights_[iVar2].begin(),
579 init: init);
580 }
581 //! \name Integration facility interface
582 //@{
583 /*! Integrates an arbitrary scalar function over the density domain(i.e.
584 computes its expected value).
585 */
586 Real integratedExpectedValue(
587 const ext::function<Real(const std::vector<Real>& v1)>& f) const {
588 // function composition: composes the integrand with the density
589 // through a product.
590 return integration()->integrate(
591 [&](const std::vector<Real>& x){ return copula_.density(x) * f(x); });
592 }
593 /*! Integrates an arbitrary vector function over the density domain(i.e.
594 computes its expected value).
595 */
596 std::vector<Real> integratedExpectedValueV(
597 // const ext::function<std::vector<Real>(
598 const ext::function<std::vector<Real>(
599 const std::vector<Real>& v1)>& f ) const {
600 detail::multiplyV M;
601 return integration()->integrateV(//see note in LMIntegrators base class
602 [&](const std::vector<Real>& x){ return M(copula_.density(x), f(x)); });
603 }
604 protected:
605 // Integrable models must provide their integrator.
606 // Arguable, not having the integration in the LM class saves that
607 // memory but have an entry in the VT...
608 virtual const ext::shared_ptr<LMIntegration>& integration() const {
609 QL_FAIL("Integration non implemented in Latent model.");
610 }
611 //@}
612
613 // Ordering is: factorWeights_[iVariable][iFactor]
614 mutable std::vector<std::vector<Real> > factorWeights_;
615 /* This is a duplicated value from the data above chosen for memory
616 reasons.
617 I have opted for this one value redundant memory rather than have the
618 memory load of the observable in all factors. Typically Latent models
619 are used in two very different ways: with many factors and not linked
620 to a market observable (typical matrix size above is of tens of
621 thousands entries) or with just one observable value and the matrix is
622 just a scalar. Otherwise, to remove the redundancy, the matrix
623 factorWeights_ should be one of Quotes Handles.
624 Yet it is not entirely true that quotes might be used only in pricing,
625 think sensitivity analysis....
626 \todo Reconsider this, see how expensive truly is.
627 */
628 mutable Handle<Quote> cachedMktFactor_;
629
630 // updated only by correlation observability and constructors.
631 // \sqrt{1-\sum_k \beta_{i,k}^2} the addition being along the factors.
632 // It has therefore the size of the basket. Cached for perfomance
633 mutable std::vector<Real> idiosyncFctrs_;
634 //! Number of systemic factors.
635 mutable Size nFactors_;//matches idiosyncFctrs_[0].size();i=0 or any
636 //! Number of latent model variables, idiosyncratic terms or model dim
637 mutable Size nVariables_;// matches idiosyncFctrs_.size()
638
639 mutable copulaType copula_;
640 };
641
642
643
644
645 // Defines ----------------------------------------------------------------
646
647#ifndef __DOXYGEN__
648
649 template <class Impl>
650 LatentModel<Impl>::LatentModel(
651 const std::vector<std::vector<Real> >& factorWeights,
652 const typename Impl::initTraits& ini)
653 : factorWeights_(factorWeights),
654 nFactors_(factorWeights[0].size()),
655 nVariables_(factorWeights.size()), copula_(factorWeights, ini)
656 {
657 for(Size i=0; i<factorWeights.size(); i++) {
658 idiosyncFctrs_.push_back(x: std::sqrt(x: 1.-
659 std::inner_product(first1: factorWeights[i].begin(),
660 last1: factorWeights[i].end(),
661 first2: factorWeights[i].begin(), init: Real(0.))));
662 // while at it, check sizes are coherent:
663 QL_REQUIRE(factorWeights[i].size() == nFactors_,
664 "Name " << i << " provides a different number of factors");
665 }
666 }
667
668 template <class Impl>
669 LatentModel<Impl>::LatentModel(
670 const std::vector<Real>& factorWeights,
671 const typename Impl::initTraits& ini)
672 : nFactors_(1),
673 nVariables_(factorWeights.size())
674 {
675 for (Real factorWeight : factorWeights)
676 factorWeights_.emplace_back(args: 1, args&: factorWeight);
677 for (Real factorWeight : factorWeights)
678 idiosyncFctrs_.push_back(x: std::sqrt(x: 1. - factorWeight * factorWeight));
679 //convert row to column vector....
680 copula_ = copulaType(factorWeights_, ini);
681 }
682
683 template <class Impl>
684 LatentModel<Impl>::LatentModel(
685 const Real correlSqr,
686 Size nVariables,
687 const typename Impl::initTraits& ini)
688 : factorWeights_(nVariables, std::vector<Real>(1, correlSqr)),
689 idiosyncFctrs_(nVariables,
690 std::sqrt(x: 1.-correlSqr*correlSqr)),
691 nFactors_(1),
692 nVariables_(nVariables),
693 copula_(factorWeights_, ini)
694 { }
695
696 template <class Impl>
697 LatentModel<Impl>::LatentModel(
698 const Handle<Quote>& singleFactorCorrel,
699 Size nVariables,
700 const typename Impl::initTraits& ini)
701 : factorWeights_(nVariables, std::vector<Real>(1,
702 std::sqrt(x: singleFactorCorrel->value()))),
703 cachedMktFactor_(singleFactorCorrel),
704 idiosyncFctrs_(nVariables,
705 std::sqrt(x: 1.-singleFactorCorrel->value())),
706 nFactors_(1),
707 nVariables_(nVariables),
708 copula_(factorWeights_, ini)
709 {
710 registerWith(h: cachedMktFactor_);
711 }
712
713#endif
714
715 template <class Impl>
716 void LatentModel<Impl>::update() {
717 /* only registration with the single market correl quote. If we get
718 register with something else remember that the quote stores correlation
719 and the model need factor values; which for one factor models are the
720 square root of the correlation.
721 */
722 factorWeights_ = std::vector<std::vector<Real> >(nVariables_,
723 std::vector<Real>(1, std::sqrt(x: cachedMktFactor_->value())));
724 idiosyncFctrs_ = std::vector<Real>(nVariables_,
725 std::sqrt(x: 1.-cachedMktFactor_->value()));
726 copula_ = copulaType(factorWeights_, copula_.getInitTraits());
727 notifyObservers();
728 }
729
730#ifndef __DOXYGEN__
731
732 //----Template partial specializations of the random FactorSampler--------
733 /*
734 Notice that while the default template needs a sequence generator the
735 specializations need a number generator. This is forced at the time the
736 concrete policy class is used in the template parameter, if it has been
737 specialized it needs the sample type typedef to match at compilation.
738
739 Notice here the outer class template is specialized only, leaving the inner
740 generator still a class template. Apparently old versions of gcc (3.x) bug
741 on this one not recognizing the specialization.
742 */
743 /*! \brief Specialization for direct Gaussian Box-Muller generation.\par
744 The implementation of Box-Muller in the library is the rejection variant so
745 do not use it within a multithreaded simulation.
746 */
747 template<class TC> template<class URNG, bool dummy>
748 class LatentModel<TC>
749 ::FactorSampler <RandomSequenceGenerator<BoxMullerGaussianRng<URNG> > ,
750 dummy> {
751 typedef URNG urng_type;
752 public:
753 //Size below must be == to the numb of factors idiosy + systemi
754 typedef Sample<std::vector<Real> > sample_type;
755 explicit FactorSampler(const GaussianCopulaPolicy& copula,
756 BigNatural seed = 0)
757 : boxMullRng_(copula.numFactors(),
758 BoxMullerGaussianRng<urng_type>(urng_type(seed))){ }
759 const sample_type& nextSequence() const {
760 return boxMullRng_.nextSequence();
761 }
762 private:
763 RandomSequenceGenerator<BoxMullerGaussianRng<urng_type> > boxMullRng_;
764 };
765
766 /*! \brief Specialization for direct T samples generation.\par
767 The PolarT is a rejection algorithm so do not use it within a
768 multithreaded simulation.
769 The RandomSequenceGenerator class does not admit heterogeneous
770 distribution samples so theres a trick here since the template parameter is
771 not what it is used internally.
772 */
773 template<class TC> template<class URNG, bool dummy>//uniform number expected
774 class LatentModel<TC>
775 ::FactorSampler<RandomSequenceGenerator<PolarStudentTRng<URNG> > ,
776 dummy> {
777 typedef URNG urng_type;
778 public:
779 typedef Sample<std::vector<Real> > sample_type;
780 explicit FactorSampler(const TCopulaPolicy& copula, BigNatural seed = 0)
781 : sequence_(std::vector<Real> (copula.numFactors()), 1.0),
782 urng_(seed) {
783 // 1 == urng.dimension() is enforced by the sample type
784 const std::vector<Real>& varF = copula.varianceFactors();
785 for (Real i : varF) // ...use back inserter lambda
786 trng_.push_back(PolarStudentTRng<urng_type>(2. / (1. - i * i), urng_));
787 }
788 const sample_type& nextSequence() const {
789 Size i=0;
790 for(; i<trng_.size(); i++)//systemic samples plus one idiosyncratic
791 sequence_.value[i] = trng_[i].next().value;
792 for(; i<sequence_.value.size(); i++)//rest of idiosyncratic samples
793 sequence_.value[i] = trng_.back().next().value;
794 return sequence_;
795 }
796 private:
797 mutable sample_type sequence_;
798 urng_type urng_;
799 mutable std::vector<PolarStudentTRng<urng_type> > trng_;
800 };
801
802#endif
803
804}
805
806
807#endif
808

source code of quantlib/ql/experimental/math/latentmodel.hpp