onefactorcopula.hpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2008 Roland Lichters
 
 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_one_factor_copula_hpp
#define quantlib_one_factor_copula_hpp
 
#include <ql/experimental/credit/distribution.hpp>
#include <ql/patterns/lazyobject.hpp>
#include <ql/quote.hpp>
#include <utility>
 
namespace QuantLib {
 
 
    class OneFactorCopula : public LazyObject {
      public:
        OneFactorCopula(Handle<Quote> correlation,
                        Real maximum = 5.0,
                        Size integrationSteps = 50,
                        Real minimum = -5.0)
        : correlation_(std::move(correlation)), max_(maximum), steps_(integrationSteps),
          min_(minimum) {
            QL_REQUIRE(correlation_->value() >= -1
                       && correlation_->value() <= 1,
                       "correlation out of range [-1, +1]");
            registerWith(correlation_);
        }
 
 
        virtual Real density(Real m) const = 0;
 
        virtual Real cumulativeZ(Real z) const = 0;
 
        virtual Real cumulativeY(Real y) const;
 
        virtual Real inverseCumulativeY(Real p) const;
 
        Real correlation() const;
 
 
        Real conditionalProbability(Real prob,
                                    Real m) const;
 
 
        std::vector<Real> conditionalProbability(const std::vector<Real>& prob,
                                                 Real m) const;
 
        Real integral(Real p) const {
            QL_REQUIRE(p >= 0 && p <= 1, "probability p=" << p
                       << " out of range [0,1]");
            calculate();
 
            Real avg = 0;
            for (Size k = 0; k < steps(); k++) {
                Real pp = conditionalProbability(p, m(k));
                avg += pp * densitydm(k);
            }
            return avg;
        }
 
        template <class F>
        Real integral(const F& f, std::vector<Real>& probabilities) const {
            calculate();
 
            Real avg = 0.0;
            for (Size i = 0; i < steps_; i++) {
                std::vector<Real> conditional
                    = conditionalProbability(probabilities, m(i));
                Real prob = f(conditional);
                avg += prob * densitydm(i);
            }
            return avg;
        }
 
        template <class F>
        Distribution integral(const F& f,
                              const std::vector<Real>& nominals,
                              const std::vector<Real>& probabilities) const {
            calculate();
 
            Distribution dist(f.buckets(), 0.0, f.maximum());
            for (Size i = 0; i < steps(); i++) {
                std::vector<Real> conditional
                    = conditionalProbability(probabilities, m(i));
                Distribution d = f(nominals, conditional);
                for (Size j = 0; j < dist.size(); j++)
                    dist.addDensity(j, d.density(j) * densitydm(i));
            }
            return dist;
        }
 
        int checkMoments(Real tolerance) const;
 
      protected:
        Handle<Quote> correlation_;
        mutable Real max_;
        mutable Size steps_;
        mutable Real min_;
 
        // Tabulated numerical solution of the cumulated distribution of Y
        mutable std::vector<Real> y_;
        mutable std::vector<Real> cumulativeY_;
 
        //private:
        // utilities for simple Euler integrations over the density of M
        Size steps() const;
 
        // i not used yet, might allow varying grid size
        // for the copula integration in the future
        Real dm(Size i) const;
 
        Real m(Size i) const;
        Real densitydm(Size i) const;
    };
 
    inline Real OneFactorCopula::correlation() const {
        calculate();
        return correlation_->value();
    }
 
    inline Size OneFactorCopula::steps() const {
        return steps_;
    }
 
    inline Real OneFactorCopula::dm(Size) const {
        return (max_ - min_)/ steps_;
    }
 
    inline Real OneFactorCopula::m(Size i) const {
        QL_REQUIRE(i < steps_, "index out of range");
        return min_ + dm(i) * i + dm(i) / 2;
    }
 
    inline Real OneFactorCopula::densitydm(Size i) const {
        QL_REQUIRE(i < steps_, "index out of range");
        return density(m(i)) * dm(i);
    }
 
}
 
#endif