lossdistribution.cpp

/* -*- 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.
*/
 
#include <ql/experimental/credit/lossdistribution.hpp>
#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
 
using namespace std;
 
namespace QuantLib {
 
    //--------------------------------------------------------------------------
    Real LossDist::binomialProbabilityOfNEvents(int n, vector<Real>& p) {
    //--------------------------------------------------------------------------
        BinomialDistribution binomial (p[0], p.size());
        return binomial(n);
    }
 
    //--------------------------------------------------------------------------
    Real LossDist::binomialProbabilityOfAtLeastNEvents(int n, vector<Real>& p) {
    //--------------------------------------------------------------------------
        CumulativeBinomialDistribution binomial(p[0], p.size());
        return 1.0 - binomial(n-1);
        /*
        Real defp = 0;
        for (Size i = n; i <= p.size(); i++)
            defp += binomialProbabilityOfNEvents (i, p);
 
        return defp;
        */
    }
 
    //--------------------------------------------------------------------------
    vector<Real> LossDist::probabilityOfNEvents(vector<Real>& p) {
    //--------------------------------------------------------------------------
        Size n = p.size();
        vector<Real> probability(n+1, 0.0);
        vector<Real> prev;
        probability[0] = 1.0;
        for (Size j = 0; j < n; j++) {
            prev = probability;
            probability[0] = prev[0] * (1.0 - p[j]);
            for (Size i = 1; i <= j; i++)
                probability[i] = prev[i-1] * p[j] + prev[i] * (1.0 - p[j]);
            probability[j+1] = prev[j] * p[j];
        }
 
        return probability;
    }
 
    //--------------------------------------------------------------------------
    Real LossDist::probabilityOfNEvents(int k, vector<Real>& p) {
    //--------------------------------------------------------------------------
        return probabilityOfNEvents(p)[k];
 
//      vector<Real> w (p.size(), 0);
//      vector<Real> u (k+1, 0);
//      vector<Real> v (k+1, 0);
 
//      Real pZero = 1.0;
//      for (Size i = 0; i < w.size(); i++) {
//          pZero *= (1.0 - p[i]);
//          w[i] = p[i] / (1.0 - p[i]);
//      }
 
//      if (k == 0) return pZero;
 
//      int kk = k;
//      Real prodw = 1.0;
 
//      Cumulated probability of up to n events:
//      Cut off when the cumulated probability reaches 1,
//      i.e. set all following probabilities of exactly n events to zero.
//      Real sum = 1.0;
 
//      u[0] = 1.0;
//      for (int i = 1; i <= kk; i++) {
//          v[i] = 0;
//          for (Size j = 0; j < w.size(); j++)
//              v[i] += pow (w[j], i);
//          u[i] = 0;
//          for (int j = 1; j <= i; j++)
//              u[i] +=  pow (-1.0, j+1) * v[j] * u[i-j];
//          u[i] /= i;
 
//          cut off
//          if (sum * pZero >= 1.0 || u[i] < 0 || u[i] * pZero >= 1.0)
//              u[i] = 0;
 
//          sum += u[i];
//      }
 
//      return pZero * prodw * u[kk];
    }
 
    //--------------------------------------------------------------------------
    Real LossDist::probabilityOfAtLeastNEvents (int k, vector<Real>& p) {
    //--------------------------------------------------------------------------
        vector<Real> probability = probabilityOfNEvents(p);
        Real sum = 1.0;
        for (int j = 0; j < k; j++)
            sum -= probability[j];
        return sum;
        /*
        Real sum = 0;
        for (Size i = k; i <= p.size(); i++)
            sum += probabilityOfNEvents (i, p);
        return sum;
        */
    }
 
    //--------------------------------------------------------------------------
    Real ProbabilityOfNEvents::operator()(vector<Real> p) const {
    //--------------------------------------------------------------------------
        return LossDist::probabilityOfNEvents (n_, p);
    }
 
    //--------------------------------------------------------------------------
    Real ProbabilityOfAtLeastNEvents::operator()(vector<Real> p) const {
    //--------------------------------------------------------------------------
        return LossDist::probabilityOfAtLeastNEvents (n_, p);
    }
 
    //--------------------------------------------------------------------------
    Real BinomialProbabilityOfAtLeastNEvents::operator()(vector<Real> p) const {
        //--------------------------------------------------------------------------
        return LossDist::binomialProbabilityOfAtLeastNEvents(n_, p);
    }
 
    //--------------------------------------------------------------------------
    Distribution LossDistBinomial::operator()(Size n, Real volume,
                                              Real probability) const {
    //--------------------------------------------------------------------------
        n_ = n;
        probability_.clear();
        probability_.resize(n_+1, 0.0);
        Distribution dist (nBuckets_, 0.0, maximum_);
        BinomialDistribution binomial (probability, n);
        for (Size i = 0; i <= n; i++) {
            if (volume_ * i <= maximum_) {
                probability_[i] = binomial(i);
                Size bucket = dist.locate(volume * i);
                dist.addDensity (bucket, probability_[i] / dist.dx(bucket));
                dist.addAverage (bucket, volume * i);
            }
        }
 
        excessProbability_.clear();
        excessProbability_.resize(n_+1, 0.0);
        excessProbability_[n_] = probability_[n_];
        for (int k = n_-1; k >= 0; k--)
            excessProbability_[k] = excessProbability_[k+1] + probability_[k];
 
        dist.normalize();
 
        return dist;
    }
 
    //--------------------------------------------------------------------------
    Distribution LossDistBinomial::operator()(const vector<Real>& nominals,
                                    const vector<Real>& probabilities) const {
    //--------------------------------------------------------------------------
        return operator()(nominals.size(), nominals[0], probabilities[0]);
    }
 
    //--------------------------------------------------------------------------
    Distribution LossDistHomogeneous::operator()(Real volume,
                                                 const vector<Real>& p) const {
    //--------------------------------------------------------------------------
        volume_ = volume;
        n_ = p.size();
        probability_.clear();
        probability_.resize(n_+1, 0.0);
        vector<Real> prev;
        probability_[0] = 1.0;
        for (Size k = 0; k < n_; k++) {
            prev = probability_;
            probability_[0] = prev[0] * (1.0 - p[k]);
            for (Size i = 1; i <= k; i++)
                probability_[i] = prev[i-1] * p[k] + prev[i] * (1.0 - p[k]);
            probability_[k+1] = prev[k] * p[k];
        }
 
        excessProbability_.clear();
        excessProbability_.resize(n_+1, 0.0);
        excessProbability_[n_] = probability_[n_];
        for (int k = n_ - 1; k >= 0; k--)
            excessProbability_[k] = excessProbability_[k+1] + probability_[k];
 
        Distribution dist (nBuckets_, 0.0, maximum_);
        for (Size i = 0; i <= n_; i++) {
            if (volume * i <= maximum_) {
                Size bucket = dist.locate(volume * i);
                dist.addDensity (bucket, probability_[i] / dist.dx(bucket));
                dist.addAverage (bucket, volume*i);
            }
        }
 
        dist.normalize();
 
        return dist;
    }
 
    //--------------------------------------------------------------------------
    Distribution LossDistHomogeneous::operator()(const vector<Real>& nominals,
                                    const vector<Real>& probabilities) const {
    //--------------------------------------------------------------------------
        return operator()(nominals[0], probabilities);
    }
 
    //--------------------------------------------------------------------------
    Distribution LossDistBucketing::operator()(const vector<Real>& nominals,
                                    const vector<Real>& probabilities) const {
    //--------------------------------------------------------------------------
        QL_REQUIRE (nominals.size() == probabilities.size(), "sizes differ: "
                    << nominals.size() << " vs " << probabilities.size());
 
        vector<Real> p (nBuckets_, 0.0);
        vector<Real> a (nBuckets_, 0.0);
        vector<Real> ap (nBuckets_, 0.0);
 
        p[0] = 1.0;
        a[0] = 0.0;
        Real dx = maximum_ / nBuckets_;
        for (Size k = 1; k < nBuckets_; k++)
            a[k] = dx * k + dx/2;
 
        for (Size i = 0; i < nominals.size(); i++) {
            Real L = nominals[i];
            Real P = probabilities[i];
            for (int k = a.size()-1; k >= 0; k--) {
                if (p[k] > 0) {
                    int u = locateTargetBucket (a[k] + L, k);
                    QL_REQUIRE (u >= 0, "u=" << u << " at i=" << i << " k=" << k);
                    QL_REQUIRE (u >= k, "u=" << u << "<k=" << k << " at i=" << i);
 
                    Real dp = p[k] * P;
                    if (u == k)
                        a[k] += P * L;
                    else {
                        // no update of a[u] and p[u] if u is beyond grid end
                        if (u < int(nBuckets_)) {
                            // a[u] remains unchanged, if dp = 0
                            if (dp > 0.0) {
                                // on Windows, p[u]/dp could cause a NaN for
                                // some very small values of p[k].
                                // Writing the above as (p[u]/p[k])/P prevents
                                // the NaN. What can I say?
                                Real f = 1.0 / (1.0 + (p[u]/p[k]) / P);
                                a[u] = (1.0 - f) * a[u] + f * (a[k] + L);
                            }
                            /* formulation of Hull-White:
                               if (p[u] + dp > 0)
                                  a[u] = (p[u] * a[u] + dp * (a[k] + L))
                                         / (p[u] + dp);
                            */
                            p[u] += dp;
                        }
                        p[k] -= dp;
                    }
                }
                QL_REQUIRE(a[k] + epsilon_ >= dx * k && a[k] < dx * (k+1),
                           "a out of range at k=" << k << ", contract " << i);
            }
        }
 
        Distribution dist (nBuckets_, 0.0, maximum_);
        for (Size i = 0; i < nBuckets_; i++) {
            dist.addDensity (i, p[i] / dx);
            dist.addAverage (i, a[i]);
        }
 
        return dist;
    }
 
    //--------------------------------------------------------------------------
    int LossDistBucketing::locateTargetBucket (Real loss, Size i0) const {
    //--------------------------------------------------------------------------
        QL_REQUIRE (loss >= 0, "loss " << loss << " must be >= 0");
        Real dx = maximum_ / nBuckets_;
        for (Size i = i0; i < nBuckets_; i++)
            if (dx * i > loss + epsilon_) return i - 1;
        return nBuckets_;
    }
 
    //--------------------------------------------------------------------------
    Distribution LossDistMonteCarlo::operator()(const vector<Real>& nominals,
                                   const vector<Real>& probabilities) const {
    //--------------------------------------------------------------------------
        Distribution dist (nBuckets_, 0.0, maximum_);
        // KnuthUniformRng rng(seed_);
        // LecuyerUniformRng rng(seed_);
        MersenneTwisterUniformRng rng(seed_);
        for (Size i = 0; i < simulations_; i++) {
            Real e = 0;
            for (Size j = 0; j < nominals.size(); j++) {
                Real r = rng.next().value;
                if (r <= probabilities[j])
                    e += nominals[j];
            }
            dist.add (e + epsilon_);
        }
 
        dist.normalize();
 
        return dist;
    }
 
}