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