Logo Search packages:      
Sourcecode: quantlib version File versions


/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

 Copyright (C) 2006 Marco Bianchetti
 Copyright (C) 2006 Silvia Frasson
 Copyright (C) 2006 Mario Pucci

 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

 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/MarketModels/driftcalculator.hpp>
#include <ql/MarketModels/duffsdeviceinnerproduct.hpp>

namespace QuantLib {

00027     DriftCalculator::DriftCalculator(const Matrix& pseudo,
                                     const std::vector<Rate>& displacements,
                                     const std::vector<Time>& taus,
                                     Size numeraire,
                                     Size alive)
    : dim_(taus.size()), factors_(pseudo.columns()),
      isFullFactor_(factors_==dim_ ? true : false),
      numeraire_(numeraire), alive_(alive),
      displacements_(displacements), oneOverTaus_(taus.size()),
      pseudo_(pseudo), tmp_(taus.size(), 0.0),
      e_(pseudo_.columns(), pseudo_.rows(), 0.0),
      downs_(taus.size()), ups_(taus.size()) {

        // Check requirements
        QL_REQUIRE(dim_>0, "Dim out of range");
        QL_REQUIRE(displacements.size() == dim_,
            "Displacements out of range");
            "pseudo.rows() not consistent with dim");
        QL_REQUIRE(pseudo.columns()>0 && pseudo.columns()<=dim_,
            "pseudo.rows() not consistent with pseudo.columns()");
        QL_REQUIRE(alive<dim_, "Alive out of bounds");
        QL_REQUIRE(numeraire_<=dim_, "Numeraire larger than dim");
        QL_REQUIRE(numeraire_>=alive, "Numeraire smaller than alive");

        // Precompute 1/taus
        Size i;
        for (i=0; i<taus.size(); ++i)
            oneOverTaus_[i] = 1.0/taus[i];

        // Compute covariance matrix from pseudoroot
        const Disposable<Matrix> pT = transpose(pseudo_);
        C_ = pseudo_*pT;

        // Compute lower and upper extrema for (non reduced) drift calculation
        for (i=alive_; i<dim_; ++i) {
            downs_[i] = std::min(i+1, numeraire_);
            ups_[i]   = std::max(i+1, numeraire_);

00068     void DriftCalculator::compute(const std::vector<Rate>& forwards,
                                  std::vector<Real>& drifts) const {
        #if defined(QL_EXTRA_SAFETY_CHECKS)
            QL_REQUIRE(forwards.size()==dim_, "forwards.size() <> dim");
            QL_REQUIRE(drifts.size()==dim_, "drifts.size() <> dim");

        if (isFullFactor_)
            computePlain(forwards, drifts);
            computeReduced(forwards, drifts);

00081     void DriftCalculator::computePlain(const std::vector<Rate>& forwards,
                                       std::vector<Real>& drifts) const {

        // Compute drifts without factor reduction,
        // using directly the covariance matrix.

        // Precompute forwards factor
        Size i;
        for(i=alive_; i<dim_; ++i)
            tmp_[i] = (forwards[i]+displacements_[i]) /
        // Compute drifts
        for (i=alive_; i<dim_; ++i) {
            drifts[i] = std::inner_product(tmp_.begin()+downs_[i],
                                           C_.row_begin(i)+downs_[i], 0.0);
            if (numeraire_>i+1)
                drifts[i] = -drifts[i];

00102     void DriftCalculator::computeReduced(const std::vector<Rate>& forwards,
                                         std::vector<Real>& drifts) const {

        // Compute drifts with factor reduction,
        // using the pseudo square root of the covariance matrix.

        // Precompute forwards factor
        Size i;
        for (i=alive_; i<dim_; ++i)
            tmp_[i] = (forwards[i]+displacements_[i]) /

        // Enforce initialization
        for (Size r=0; r<factors_; ++r)
            e_[r][std::max(0,static_cast<Integer>(numeraire_)-1)] = 0.0;

        // Now compute drifts: take the numeraire P_N (numeraire_=N)
        // as the reference point, divide the summation into 3 steps,
        // et impera:

        // 1st step: the drift corresponding to the numeraire P_N is zero.
        // (if N=0 no drift is null, if N=dim_ the last drift is null).
        if (numeraire_>0) drifts[numeraire_-1] = 0.0;

        // 2nd step: then, move backward from N-2 (included) back to
        // alive (included) (if N=0 jumps to 3rd step, if N=dim_ the
        // e_[r][N-1] are correctly initialized):
        for (Integer j=static_cast<Integer>(numeraire_)-2;
             j>=static_cast<Integer>(alive_); --j) {
            for (Size r=0; r<factors_; ++r) {
                e_[r][j] = e_[r][j+1] + tmp_[j+1] * pseudo_[j+1][r];
            drifts[j] = - std::inner_product(e_.column_begin(j),

        // 3rd step: now, move forward from N (included) up to n (excluded)
        // (if N=0 this is the only relevant computation):
        for (i=numeraire_; i<dim_; ++i) {
            for (Size r=0; r<factors_; ++r) {
                if (i==0) {
                    e_[r][i] = tmp_[i] * pseudo_[i][r];
                } else {
                    e_[r][i] = e_[r][i-1] + tmp_[i] * pseudo_[i][r];
            drifts[i] = std::inner_product(e_.column_begin(i),


Generated by  Doxygen 1.6.0   Back to index