latentmodel.hpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2008 Roland Lichters
 Copyright (C) 2014 Jose Aparicio
 
 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/
 
 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.
 
 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/
 
#ifndef quantlib_latent_model_hpp
#define quantlib_latent_model_hpp
 
#include <ql/experimental/math/multidimquadrature.hpp>
#include <ql/experimental/math/multidimintegrator.hpp>
#include <ql/math/integrals/trapezoidintegral.hpp>
#include <ql/math/randomnumbers/randomsequencegenerator.hpp>
#include <ql/experimental/math/gaussiancopulapolicy.hpp>
#include <ql/experimental/math/tcopulapolicy.hpp>
#include <ql/math/randomnumbers/boxmullergaussianrng.hpp>
#include <ql/experimental/math/polarstudenttrng.hpp>
#include <ql/handle.hpp>
#include <ql/quote.hpp>
#include <vector>
 
namespace QuantLib {
 
    namespace detail {
        // havent figured out how to do this in-place
        struct multiplyV {
            QL_DEPRECATED
            typedef std::vector<Real> result_type;
 
            std::vector<Real> operator()(Real d, std::vector<Real> v) 
            {
                std::transform(v.begin(), v.end(), v.begin(), 
                               [=](Real x) -> Real { return x * d; });
                return v;
            }
        };
    }
 
 
    /* Things trying to achieve here:
    1.- Unify the two branches of integrators in the library, they do not 
      hang from a common base class and here a common ptr for the 
      factory is needed.
    2.- Have a common signature for the integration call.
    3.- Factory construction so integrable latent models can choose the 
      integration algorithm separately.
    */
    class LMIntegration {
    public:
        // Interface with actual integrators:
        // integral of a scalar function
        virtual Real integrate(const ext::function<Real (
            const std::vector<Real>& arg)>& f) const = 0;
        // integral of a vector function
        /* I had to use a different name, since the compiler does not
        recognise the overload; MSVC sees the argument as 
        ext::function<Signature> in both cases....   
        I could do the as with the quadratures and have this as a template 
        function and spez for the vector case but I prefer to understand
        why the overload fails....
                    FIX ME
        */
        virtual std::vector<Real> integrateV(
            const ext::function<std::vector<Real>  (
            const std::vector<Real>& arg)>& f) const {
            QL_FAIL("No vector integration provided");
        }
        virtual ~LMIntegration() = default;
    };
 
    //CRTP-ish for joining the integrations, class above to have the factory
    template <class I_T>
    class IntegrationBase : 
        public I_T, public LMIntegration {// diamond on 'integrate'
     // this class template always to be fully specialized:
     private:
       IntegrationBase() = default;
    };
    
    // gcc reports value collision with heston engine (?!) thats why the name
    namespace LatentModelIntegrationType {
        typedef 
        enum LatentModelIntegrationType {
            #ifndef QL_PATCH_SOLARIS
            GaussianQuadrature,
            #endif
            Trapezoid
            // etc....
        } LatentModelIntegrationType;
    }
 
    #ifndef QL_PATCH_SOLARIS
 
    /* class template specializations. I havent use CRTP type cast directly
    because the signature of the integrators is different, grid integration
    needs the domain. */
    template<> class IntegrationBase<GaussianQuadMultidimIntegrator> : 
    public GaussianQuadMultidimIntegrator, public LMIntegration {
    public:
        IntegrationBase(Size dimension, Size order) 
        : GaussianQuadMultidimIntegrator(dimension, order) {}
        Real integrate(const ext::function<Real(const std::vector<Real>& arg)>& f) const override {
            return GaussianQuadMultidimIntegrator::integrate<Real>(f);
        }
        std::vector<Real> integrateV(
            const ext::function<std::vector<Real>(const std::vector<Real>& arg)>& f)
            const override {
            return GaussianQuadMultidimIntegrator::integrate<std::vector<Real>>(f);
        }
        ~IntegrationBase() override = default;
    };
 
    #endif
 
    template<> class IntegrationBase<MultidimIntegral> : 
        public MultidimIntegral, public LMIntegration {
    public:
        IntegrationBase(
            const std::vector<ext::shared_ptr<Integrator> >& integrators, 
            Real a, Real b) 
        : MultidimIntegral(integrators), 
          a_(integrators.size(),a), b_(integrators.size(),b) {}
        Real integrate(const ext::function<Real(const std::vector<Real>& arg)>& f) const override {
            return MultidimIntegral::operator()(f, a_, b_);
        }
        // vector version here....
        ~IntegrationBase() override = default;
        const std::vector<Real> a_, b_;
    };
 
    // Intended to replace OneFactorCopula
 
    template <class copulaPolicyImpl>
    class LatentModel 
        : public virtual Observer , public virtual Observable 
    {//observer if factors as quotes
    public:
      void update() override;
 
      typedef copulaPolicyImpl copulaType;
      Probability cumulativeY(Real val, Size iVariable) const {
          return copula_.cumulativeY(val, iVariable);
        }
        Probability cumulativeZ(Real z) const {
            return copula_.cumulativeZ(z);
        }
        Probability density(const std::vector<Real>& m) const {
            #if defined(QL_EXTRA_SAFETY_CHECKS)
                QL_REQUIRE(m.size() == nFactors_, 
                    "Factor size must match that of model.");
            #endif
            return copula_.density(m);
        }
        Real inverseCumulativeDensity(Probability p, Size iFactor) const {
            return copula_.inverseCumulativeDensity(p, iFactor);
        }
        Real inverseCumulativeY(Probability p, Size iVariable) const {
            return copula_.inverseCumulativeY(p, iVariable);
        }
        Real inverseCumulativeZ(Probability p) const {
            return copula_.inverseCumulativeZ(p);
        }
        std::vector<Real> allFactorCumulInverter(const std::vector<Real>& probs) const {
            return copula_.allFactorCumulInverter(probs);
        }
 
        Real latentVarValue(const std::vector<Real>& allFactors, 
                            Size iVar) const 
        {
            return std::inner_product(factorWeights_[iVar].begin(), 
                // systemic term:
                factorWeights_[iVar].end(), allFactors.begin(),
                // idiosyncratic term:
                Real(allFactors[numFactors()+iVar] * idiosyncFctrs_[iVar]));
        }
        // \to do write variants of the above, although is the most common case
 
        const copulaType& copula() const {
            return copula_;
        }
 
 
    //  protected:
 
 
        /*
            Several (very different) usages make the spez non trivial
            The final goal is to obtain a sequence generator of the factor 
            samples, several routes are possible depending on the algorithms:
            
            1.- URNG -> Sequence Gen -> CopulaInversion  
              e.g.: CopulaInversion(RandomSequenceGenerator<MersenneTwisterRNG>)
            2.- PseudoRSG ------------> CopulaInversion
              e.g.: CopulaInversion(SobolRSG)
            3.- URNG -> SpecificMapping -> Sequence Gen  (bypasses the copula 
                for performance)
              e.g.: RandomSequenceGenerator<BoxMullerGaussianRng<
                MersenneTwisterRNG> > 
            
            Notice that the order the three algorithms involved (uniform gen, 
            sequence construction, distribution mapping) is not always the same.
            (in fact there could be some other ways to generate but these are 
            the ones in the library now.)
            Difficulties arise when wanting to use situation 3.- whith a generic
            RNG, leaving it unspecified
            
            Derived classes might specialize (on the copula
            type) to another type of generator if a more efficient algorithm 
            that the distribution inversion is available; rewritig then the 
            nextSequence method for a particular copula implementation.
            Some combinations of generators might make no sense, while it 
            could be possible to block template classes corresponding to those
            cases its not done (yet?) (e.g. a BoxMuller under a TCopula.)
            Dimensionality coherence (between the generator and the copula) 
            should have been checked by the client code.
            In multithread usage the sequence generator is expect to be already
            in position.
            To sample the latent variable itself users should call 
            LatentModel::latentVarValue with these samples.
        */
        // Cant use InverseCumulativeRsg since the inverse there has to return a
        //   real number and here a vector is needed, the function inverted here
        //   is multivalued.
        template <class USNG, 
            // dummy template parameter to allow for 'full' specialization of 
            // inner class without specialization of the outer.
            bool = true>
        class FactorSampler {
        public:
            typedef Sample<std::vector<Real> > sample_type;
            explicit FactorSampler(const copulaType& copula, 
                BigNatural seed = 0) 
            : sequenceGen_(copula.numFactors(), seed), // base case construction
              x_(std::vector<Real>(copula.numFactors()), 1.0),
              copula_(copula) { }
            const sample_type& nextSequence() const {
                typename USNG::sample_type sample =
                    sequenceGen_.nextSequence();
                x_.value = copula_.allFactorCumulInverter(sample.value);
                return x_;
            }
        private:
            USNG sequenceGen_;// copy, we might be mutithreaded
            mutable sample_type x_;
            // no copies
            const copulaType& copula_;
        };
    protected:
        /* \todo Move integrator traits like number of quadrature points, 
        integration domain dimensions, etc to the copula through a static 
        member function. Since they depend on the nature of the probability 
        density distribution thats where they belong.
        This is why theres one factory per copula policy template parameter 
        (even if this is not used...yet)
        */
        class IntegrationFactory {
        public:
            static ext::shared_ptr<LMIntegration> createLMIntegration(
                Size dimension, 
                LatentModelIntegrationType::LatentModelIntegrationType type = 
                    #ifndef QL_PATCH_SOLARIS
                    LatentModelIntegrationType::GaussianQuadrature)
                    #else
                    LatentModelIntegrationType::Trapezoid)
                    #endif
            {
                switch(type) {
                    #ifndef QL_PATCH_SOLARIS
                    case LatentModelIntegrationType::GaussianQuadrature:
                        return 
                            ext::make_shared<
                            IntegrationBase<GaussianQuadMultidimIntegrator> >(
                                dimension, 25);
                    #endif
                    case LatentModelIntegrationType::Trapezoid:
                        {
                        std::vector<ext::shared_ptr<Integrator> > integrals;
                        for(Size i=0; i<dimension; i++)
                            integrals.push_back(
                            ext::make_shared<TrapezoidIntegral<Default> >(
                                1.e-4, 20));
                        /* This integration domain is tailored for the T 
                        distribution; it is too wide for normals or Ts of high
                        order. 
                        \todo This needs to be solved by having the copula to 
                        provide the integration traits for any integration 
                        algorithm since it is the copula that knows the relevant
                        domain for its density distributions. Also to be able to
                        block integrations which will fail; like a quadrature  
                        here in some cases.
                        */
                        return 
                          ext::make_shared<IntegrationBase<MultidimIntegral> >
                               (integrals, -35., 35.);
                        }
                    default:
                        QL_FAIL("Unknown latent model integration type.");
                }
            }
        private:
          IntegrationFactory() = default;
        };
 
 
    public:
        // model size, number of latent variables modelled
        Size size() const {return nVariables_;}
        Size numFactors() const {return nFactors_;}
        Size numTotalFactors() const { return nVariables_ + nFactors_; }
 
        explicit LatentModel(
            const std::vector<std::vector<Real> >& factorsWeights, 
            const typename copulaType::initTraits& ini = 
                copulaType::initTraits());
        explicit LatentModel(const std::vector<Real>& factorsWeight,
            const typename copulaType::initTraits& ini = 
                copulaType::initTraits());
        explicit LatentModel(Real correlSqr,
                             Size nVariables,
                             const typename copulaType::initTraits& ini = copulaType::initTraits());
        explicit LatentModel(const Handle<Quote>& singleFactorCorrel,
            Size nVariables,
            const typename copulaType::initTraits& ini = 
                copulaType::initTraits());
 
        const std::vector<std::vector<Real> >& factorWeights() const {
            return factorWeights_;
        }
        const std::vector<Real>& idiosyncFctrs() const {return idiosyncFctrs_;}
 
        Real latentVariableCorrel(Size iVar1, Size iVar2) const {
            // true for any normalized combination
            Real init = (iVar1 == iVar2 ? 
                idiosyncFctrs_[iVar1] * idiosyncFctrs_[iVar1] : Real(0.));
            return std::inner_product(factorWeights_[iVar1].begin(), 
                factorWeights_[iVar1].end(), factorWeights_[iVar2].begin(), 
                    init);
        }
 
 
        Real integratedExpectedValue(
            const ext::function<Real(const std::vector<Real>& v1)>& f) const {
            // function composition: composes the integrand with the density 
            //   through a product.
            return integration()->integrate(
                [&](const std::vector<Real>& x){ return copula_.density(x) * f(x); });
        }
        std::vector<Real> integratedExpectedValueV(
            // const ext::function<std::vector<Real>(
            const ext::function<std::vector<Real>(
                const std::vector<Real>& v1)>& f ) const {
            detail::multiplyV M;
            return integration()->integrateV(//see note in LMIntegrators base class
                [&](const std::vector<Real>& x){ return M(copula_.density(x), f(x)); });
        }
    protected:
        // Integrable models must provide their integrator.
        // Arguable, not having the integration in the LM class saves that 
        //   memory but have an entry in the VT... 
        virtual const ext::shared_ptr<LMIntegration>& integration() const {
            QL_FAIL("Integration non implemented in Latent model.");
        }
 
        // Ordering is: factorWeights_[iVariable][iFactor]
        mutable std::vector<std::vector<Real> > factorWeights_;
        /* This is a duplicated value from the data above chosen for memory 
        reasons.
        I have opted for this one value redundant memory rather than have the 
        memory load of the observable in all factors. Typically Latent models 
        are used in two very different ways: with many factors and not linked 
        to a market observable (typical matrix size above is of tens of 
        thousands entries) or with just one observable value and the matrix is 
        just a scalar. Otherwise, to remove the redundancy, the matrix 
        factorWeights_ should be one of Quotes Handles.
        Yet it is not entirely true that quotes might be used only in pricing, 
        think sensitivity analysis....
        \todo Reconsider this, see how expensive truly is.
        */
        mutable Handle<Quote> cachedMktFactor_;
 
        // updated only by correlation observability and constructors.
        // \sqrt{1-\sum_k \beta_{i,k}^2} the addition being along the factors. 
        // It has therefore the size of the basket. Cached for perfomance
        mutable std::vector<Real> idiosyncFctrs_;
        mutable Size nFactors_;//matches idiosyncFctrs_[0].size();i=0 or any
        mutable Size nVariables_;// matches idiosyncFctrs_.size() 
 
        mutable copulaType copula_;
    };
 
 
 
 
    // Defines ----------------------------------------------------------------
 
#ifndef __DOXYGEN__
 
    template <class Impl>
    LatentModel<Impl>::LatentModel(
        const std::vector<std::vector<Real> >& factorWeights,
        const typename Impl::initTraits& ini)
    : factorWeights_(factorWeights),
      nFactors_(factorWeights[0].size()), 
      nVariables_(factorWeights.size()), copula_(factorWeights, ini)
    {
        for(Size i=0; i<factorWeights.size(); i++) {
            idiosyncFctrs_.push_back(std::sqrt(1.-
                    std::inner_product(factorWeights[i].begin(), 
                factorWeights[i].end(), 
                factorWeights[i].begin(), Real(0.))));
            // while at it, check sizes are coherent:
            QL_REQUIRE(factorWeights[i].size() == nFactors_, 
                "Name " << i << " provides a different number of factors");
        }
    }
 
    template <class Impl>
    LatentModel<Impl>::LatentModel(
        const std::vector<Real>& factorWeights,
        const typename Impl::initTraits& ini)
    : nFactors_(1),
      nVariables_(factorWeights.size())
    {
        for (Real factorWeight : factorWeights)
            factorWeights_.emplace_back(1, factorWeight);
        for (Real factorWeight : factorWeights)
            idiosyncFctrs_.push_back(std::sqrt(1. - factorWeight * factorWeight));
        //convert row to column vector....
        copula_ = copulaType(factorWeights_, ini);
    }
 
    template <class Impl>
    LatentModel<Impl>::LatentModel(
        const Real correlSqr,
        Size nVariables,
        const typename Impl::initTraits& ini)
    : factorWeights_(nVariables, std::vector<Real>(1, correlSqr)),
      idiosyncFctrs_(nVariables, 
        std::sqrt(1.-correlSqr*correlSqr)),
      nFactors_(1), 
      nVariables_(nVariables),
      copula_(factorWeights_, ini)
    { }
 
    template <class Impl>
    LatentModel<Impl>::LatentModel(
        const Handle<Quote>& singleFactorCorrel,
        Size nVariables,
        const typename Impl::initTraits& ini)
    : factorWeights_(nVariables, std::vector<Real>(1, 
        std::sqrt(singleFactorCorrel->value()))),
      cachedMktFactor_(singleFactorCorrel),
      idiosyncFctrs_(nVariables, 
        std::sqrt(1.-singleFactorCorrel->value())),
      nFactors_(1), 
      nVariables_(nVariables),
      copula_(factorWeights_, ini)
    {
        registerWith(cachedMktFactor_);
    }
 
#endif
 
    template <class Impl>
    void LatentModel<Impl>::update() {
        /* only registration with the single market correl quote. If we get 
        register with something else remember that the quote stores correlation
        and the model need factor values; which for one factor models are the
        square root of the correlation.
        */
        factorWeights_ = std::vector<std::vector<Real> >(nVariables_, 
            std::vector<Real>(1, std::sqrt(cachedMktFactor_->value())));
        idiosyncFctrs_ = std::vector<Real>(nVariables_, 
            std::sqrt(1.-cachedMktFactor_->value()));
        copula_ = copulaType(factorWeights_, copula_.getInitTraits());
        notifyObservers();
    }
 
#ifndef __DOXYGEN__
 
    //----Template partial specializations of the random FactorSampler--------
    /*
    Notice that while the default template needs a sequence generator the 
    specializations need a number generator. This is forced at the time the 
    concrete policy class is used in the template parameter, if it has been 
    specialized it needs the sample type typedef to match at compilation. 
    
    Notice here the outer class template is specialized only, leaving the inner
    generator still a class template. Apparently old versions of gcc (3.x) bug 
    on this one not recognizing the specialization.
    */
    template<class TC> template<class URNG, bool dummy>
    class LatentModel<TC>
        ::FactorSampler <RandomSequenceGenerator<BoxMullerGaussianRng<URNG> > ,
            dummy> {
        typedef URNG urng_type;
    public:
        //Size below must be == to the numb of factors idiosy + systemi
        typedef Sample<std::vector<Real> > sample_type;
        explicit FactorSampler(const GaussianCopulaPolicy& copula,
                               BigNatural seed = 0) 
        : boxMullRng_(copula.numFactors(), 
            BoxMullerGaussianRng<urng_type>(urng_type(seed))){ }
        const sample_type& nextSequence() const {
                return boxMullRng_.nextSequence();
        }
    private:
        RandomSequenceGenerator<BoxMullerGaussianRng<urng_type> > boxMullRng_;
    };
 
    template<class TC> template<class URNG, bool dummy>//uniform number expected
    class LatentModel<TC>
        ::FactorSampler<RandomSequenceGenerator<PolarStudentTRng<URNG> > , 
            dummy> {
        typedef URNG urng_type;
    public:
        typedef Sample<std::vector<Real> > sample_type;
        explicit FactorSampler(const TCopulaPolicy& copula, BigNatural seed = 0)
        : sequence_(std::vector<Real> (copula.numFactors()), 1.0),
          urng_(seed) {
            // 1 == urng.dimension() is enforced by the sample type
            const std::vector<Real>& varF = copula.varianceFactors();
            for (Real i : varF) // ...use back inserter lambda
                trng_.push_back(PolarStudentTRng<urng_type>(2. / (1. - i * i), urng_));
        }
        const sample_type& nextSequence() const {
            Size i=0;
            for(; i<trng_.size(); i++)//systemic samples plus one idiosyncratic
                sequence_.value[i] = trng_[i].next().value;
            for(; i<sequence_.value.size(); i++)//rest of idiosyncratic samples
                sequence_.value[i] = trng_.back().next().value;
            return sequence_;
        }
    private:
        mutable sample_type sequence_;
        urng_type urng_;
        mutable std::vector<PolarStudentTRng<urng_type> > trng_;
    };
 
#endif
 
}                    
 
 
#endif