#include <orea/engine/parsensitivityanalysis.hpp>

Collaboration diagram for ParSensitivityConverter:

Public Member Functions
	ParSensitivityConverter (const ParSensitivityAnalysis::ParContainer &parSensitivities, const std::map< ore::analytics::RiskFactorKey, std::pair< QuantLib::Real, QuantLib::Real > > &shiftSizes)

const std::set< ore::analytics::RiskFactorKey > &	rawKeys ()
	Inspectors. More...

const std::set< ore::analytics::RiskFactorKey > &	parKeys ()

boost::numeric::ublas::vector< Real >	convertSensitivity (const boost::numeric::ublas::vector< Real > &zeroSensitivities)
	Takes an array of zero sensitivities and returns an array of par sensitivities. More...

void	writeConversionMatrix (ore::data::Report &reportOut) const
	Write the inverse of the transposed Jacobian to the `reportOut`. More...

ParSensitivityAnalysis::ParContainer	inverseJacobian () const

Private Attributes
std::set< ore::analytics::RiskFactorKey >	rawKeys_

std::set< ore::analytics::RiskFactorKey >	parKeys_

QuantLib::SparseMatrix	jacobi_transp_inv_

boost::numeric::ublas::vector< QuantLib::Real >	zeroShifts_
	Vector of absolute zero shift sizes. More...

boost::numeric::ublas::vector< QuantLib::Real >	parShifts_
	Vector of absolute par shift sizes. More...

Detailed Description

ParSensitivityConverter class.

1) Build Jacobi matrix containing sensitivities of par rates (first index) w.r.t. zero shifts (second index) 2) Convert zero rate and optionlet vol sensitivities into par rate/vol sensitivities

Let: p_ij denote curve i's par instrument j (curve may be a discount or an index curve) z_kl denote curve k's zero shift l (curve may be a discount or an index curve) i,k run from 0 to I j,l run from 0 to J

Organisation of the matrix: z_00 z_01 z_02 ... z_0J ... z_10 z_11 ... z_1J ... z_I0 z_I1 ... z_IJ p_00 p_01 p_02 ... p_0J p_10 p_11 ... p_1J ... p_I0 p_I1 ... p_IJ

The curve numbering starts with discount curves (by ccy), followed by index curves (by index name). The number J of par instruments respectively zero shifts can differ between discount and index curves. The number of zero shifts matches the number of par instruments. The Jacobi matrix is therefore quadratic by construction.

Definition at line 174 of file parsensitivityanalysis.hpp.

Constructor & Destructor Documentation

◆ ParSensitivityConverter()

ParSensitivityConverter	(	const ParSensitivityAnalysis::ParContainer &	parSensitivities,
		const std::map< ore::analytics::RiskFactorKey, std::pair< QuantLib::Real, QuantLib::Real > > &	shiftSizes
	)

Constructor where parSensitivities is the par rate sensitivities w.r.t. zero shifts shiftSizes gives the absolute zero and par shift sizes for each risk factor key.

Definition at line 572 of file parsensitivityanalysis.cpp.

                                                                                                         {
 
    // Populate the set of par keys (rows of Jacobi) and raw zero keys (columns of Jacobi)
    for (auto parEntry : parSensitivities) {
        parKeys_.insert(parEntry.first.first);
        rawKeys_.insert(parEntry.first.second);
    }
 
    Size n_par = parKeys_.size();
    Size n_raw = rawKeys_.size();
    SparseMatrix jacobi_transp(n_raw, n_par); // transposed Jacobi
    LOG("Transposed Jacobi matrix dimension " << n_raw << " x " << n_par);
    if (parKeys_ != rawKeys_) {
        std::set<RiskFactorKey> parMinusRaw, rawMinusPar;
        std::set_difference(parKeys_.begin(), parKeys_.end(), rawKeys_.begin(), rawKeys_.end(),
                            std::inserter(parMinusRaw, parMinusRaw.begin()));
        std::set_difference(rawKeys_.begin(), rawKeys_.end(), parKeys_.begin(), parKeys_.end(),
                            std::inserter(rawMinusPar, rawMinusPar.begin()));
        for (auto const& p : parMinusRaw) {
            ALOG("par key '" << p << "' not in raw key set");
        }
        for (auto const& p : rawMinusPar) {
            ALOG("raw key '" << p << "' not in par key set");
        }
        QL_FAIL("Zero and par keys should be equal for par conversion, see log for differences");
    }
    QL_REQUIRE(n_raw > 0, "Transposed Jacobi matrix has size 0");
 
    LOG("Populating the vector of zero and par shift sizes");
    zeroShifts_.resize(n_raw);
    parShifts_.resize(n_par);
    Size i = 0;
    for (const auto& key : rawKeys_) {
        QL_REQUIRE(shiftSizes.count(key) == 1, "ParSensitivityConverter is missing shift sizes for key " << key);
        // avoid division by zero below: if we have a zero shift, the corresponding zero sensi will be zero, too,
        // but if we divide this sensi by zero, we will get nan instead of zero
        zeroShifts_[i] = std::max(shiftSizes.at(key).first, 1E-10);
        parShifts_[i] = shiftSizes.at(key).second;
        i++;
    }
 
    LOG("Populating Transposed Jacobi matrix");
    for (auto const& p : parSensitivities) {
        Size parIdx = std::distance(parKeys_.begin(), parKeys_.find(p.first.first));
        Size rawIdx = std::distance(rawKeys_.begin(), rawKeys_.find(p.first.second));
        QL_REQUIRE(parIdx < n_par, "internal error: parKey " << p.first.first << " not found in parKeys_");
        QL_REQUIRE(rawIdx < n_raw, "internal error: rawKey " << p.first.second << " not found in parKeys_");
        jacobi_transp(rawIdx, parIdx) = p.second;
        TLOG("Matrix entry [" << rawIdx << ", " << parIdx << "] ~ [raw:" << p.first.second << ", par:" << p.first.first
                              << "]: " << p.second);
    }
    LOG("Finished populating transposed Jacobi matrix, non-zero entries = "
        << parSensitivities.size() << " ("
        << 100.0 * static_cast<Real>(parSensitivities.size()) / static_cast<Real>(n_par * n_raw) << "%)");
 
    LOG("Populating block indices");
    vector<Size> blockIndices;
    pair<RiskFactorKey::KeyType, string> previousGroup(RiskFactorKey::KeyType::None, "");
    Size blockIndex = 0;
    for (auto r : rawKeys_) {
        // Update block indices with the index where the risk factor group changes
        // e.g. when move from (DiscountCurve, EUR) to (DiscountCurve, USD)
        pair<RiskFactorKey::KeyType, string> thisGroup(r.keytype, r.name);
        if (blockIndex > 0 && previousGroup != thisGroup) {
            blockIndices.push_back(blockIndex);
            TLOG("Adding block index " << blockIndex);
        }
        ++blockIndex;
 
        // Update the risk factor group
        previousGroup = thisGroup;
    }
    blockIndices.push_back(blockIndex);
    TLOG("Adding block index " << blockIndex);
    LOG("Finished Populating block indices.");
 
    LOG("Invert Transposed Jacobi matrix");
    bool success = true;
    try {
        jacobi_transp_inv_ = blockMatrixInverse(jacobi_transp, blockIndices);
    } catch (const std::exception& e) {
        // something went wrong during the matrix inversion, so we run an extended analysis on the original matrix
        // to see whether there are zero or linearly dependent rows / columns
        StructuredAnalyticsErrorMessage("Par sensitivity conversion", "Transposed Jacobi matrix inversion failed",
                                        e.what())
            .log();
        LOG("Running extended matrix diagnostics (looking for zero or linearly dependent rows / columns...)");
        constexpr Size nOp = 1000; // number of operations for close_enough comparisons below
        LOG("Checking for zero rows...");
        for (auto it1 = jacobi_transp.begin1(); it1 != jacobi_transp.end1(); ++it1) {
            Real tmp = 0.0;
            for (auto it2 = it1.begin(); it2 != it1.end(); ++it2) {
                tmp += (*it2) * (*it2);
            }
            if (close_enough(tmp, 0.0, n_par * nOp)) {
                WLOG("row " << it1.index1() << " is zero");
            }
        }
        LOG("Checking for zero columns...");
        for (auto it2 = jacobi_transp.begin2(); it2 != jacobi_transp.end2(); ++it2) {
            Real tmp = 0.0;
            for (auto it1 = it2.begin(); it1 != it2.end(); ++it1) {
                tmp += (*it1) * (*it1);
            }
            if (close_enough(tmp, 0.0, n_par * nOp)) {
                WLOG("column " << it2.index2() << " is zero");
            }
        }
        LOG("Checking for linearly dependent rows...");
        for (auto it1 = jacobi_transp.begin1(); it1 != jacobi_transp.end1(); ++it1) {
            for (auto it1b = jacobi_transp.begin1(); it1b != jacobi_transp.end1(); ++it1b) {
                if (it1b.index1() <= it1.index1())
                    continue;
                Real ratio = Null<Real>();
                if (it1.begin() != it1.end() && it1b.begin() != it1b.end()) {
                    bool dependent = true;
                    auto it2b = it1b.begin();
                    for (auto it2 = it1.begin(); it2 != it1.end() && dependent; ++it2) {
                        if (close_enough(*it2, 0.0, nOp))
                            continue;
                        bool foundMatchingIndex = false;
                        while (it2b != it1b.end() && it2b.index2() <= it2.index2() && dependent) {
                            if (it2b.index2() < it2.index2()) {
                                if (!close_enough(*it2b, 0.0, nOp))
                                    dependent = false;
                            } else {
                                foundMatchingIndex = true;
                                if (close_enough(*it2b, 0.0, nOp)) {
                                    dependent = false;
                                } else if (ratio == Null<Real>()) {
                                    ratio = *it2b / *it2;
                                } else if (!close_enough(*it2b / *it2, ratio, nOp)) {
                                    dependent = false;
                                }
                            }
                            ++it2b;
                        }
                        if (!foundMatchingIndex)
                            dependent = false;
                    }
                    while (it2b != it1b.end() && dependent) {
                        if (!close_enough(*it2b, 0.0))
                            dependent = false;
                        ++it2b;
                    }
                    if (dependent) {
                        WLOG("rows " << it1.index1() << " and " << it1b.index1() << " are linearly dependent.");
                    }
                }
            }
        }
        LOG("Checking for linearly dependent columns...");
        for (auto it1 = jacobi_transp.begin2(); it1 != jacobi_transp.end2(); ++it1) {
            for (auto it1b = jacobi_transp.begin2(); it1b != jacobi_transp.end2(); ++it1b) {
                if (it1b.index2() <= it1.index2())
                    continue;
                Real ratio = Null<Real>();
                if (it1.begin() != it1.end() && it1b.begin() != it1b.end()) {
                    bool dependent = true;
                    auto it2b = it1b.begin();
                    for (auto it2 = it1.begin(); it2 != it1.end() && dependent; ++it2) {
                        if (close_enough(*it2, 0.0, nOp))
                            continue;
                        bool foundMatchingIndex = false;
                        while (it2b != it1b.end() && it2b.index1() <= it2.index1() && dependent) {
                            if (it2b.index1() < it2.index1()) {
                                if (!close_enough(*it2b, 0.0, nOp))
                                    dependent = false;
                            } else {
                                foundMatchingIndex = true;
                                if (close_enough(*it2b, 0.0, nOp)) {
                                    dependent = false;
                                } else if (ratio == Null<Real>()) {
                                    ratio = *it2b / *it2;
                                } else if (!close_enough(*it2b / *it2, ratio, nOp)) {
                                    dependent = false;
                                }
                            }
                            ++it2b;
                        }
                        if (!foundMatchingIndex)
                            dependent = false;
                    }
                    while (it2b != it1b.end() && dependent) {
                        if (!close_enough(*it2b, 0.0))
                            dependent = false;
                        ++it2b;
                    }
                    if (dependent) {
                        WLOG("columns " << it1.index2() << " and " << it1b.index2() << " are linearly dependent.");
                    }
                }
            }
        }
        LOG("Extended matrix diagnostics done. Exiting application.");
        success = false;
    }
    QL_REQUIRE(success, "Jacobi matrix inversion failed, see log file for more details.");
    Real conditionNumber = modifiedMaxNorm(jacobi_transp) * modifiedMaxNorm(jacobi_transp_inv_);
    LOG("Inverse Jacobi done, condition number of Jacobi matrix is " << conditionNumber);
    DLOG("Diagonal entries of Jacobi and inverse Jacobi:");
    DLOG("row/col              Jacobi             Inverse");
    for (Size j = 0; j < jacobi_transp.size1(); ++j) {
        DLOG(right << setw(7) << j << setw(20) << jacobi_transp(j, j) << setw(20) << jacobi_transp_inv_(j, j));
    }
}

Here is the call graph for this function:

Member Function Documentation

◆ rawKeys()

const std::set< ore::analytics::RiskFactorKey > & rawKeys ( )

Inspectors.

Return the set of raw, i.e. zero, risk factor keys The ordering in this set defines the order of the columns in the Jacobi matrix

Definition at line 187 of file parsensitivityanalysis.hpp.

187{ return rawKeys_; }

◆ parKeys()

const std::set< ore::analytics::RiskFactorKey > & parKeys ( )

Return the set of par risk factor keys The ordering in this set defines the order of the rows in the Jacobi matrix

Definition at line 190 of file parsensitivityanalysis.hpp.

190{ return parKeys_; }

◆ convertSensitivity()

boost::numeric::ublas::vector< Real > convertSensitivity ( const boost::numeric::ublas::vector< Real > & zeroSensitivities )

Takes an array of zero sensitivities and returns an array of par sensitivities.

Parameters

zeroSensitivities array of zero sensitivities ordered according to rawKeys()

Returns: array of par sensitivities ordered according to parKeys()

Definition at line 781 of file parsensitivityanalysis.cpp.

                                                                                                  {
 
    DLOG("Start sensitivity conversion");
 
    Size dim = zeroSensitivities.size();
    QL_REQUIRE(jacobi_transp_inv_.size1() == dim, "Size mismatch between Transoposed Jacobi inverse matrix ["
                                                      << jacobi_transp_inv_.size1() << " x "
                                                      << jacobi_transp_inv_.size2() << "] and zero sensitivity array ["
                                                      << dim << "]");
 
    // Vector storing approximation for \frac{\partial V}{\partial z_i} for each zero factor z_i
    boost::numeric::ublas::vector<Real> zeroDerivs(dim);
    zeroDerivs = element_div(zeroSensitivities, zeroShifts_);
 
    // Vector initially storing approximation for \frac{\partial V}{\partial c_i} for each par factor c_i
    boost::numeric::ublas::vector<Real> parSensitivities(dim);
    boost::numeric::ublas::axpy_prod(jacobi_transp_inv_, zeroDerivs, parSensitivities, true);
 
    // Update parSensitivities vector to hold the first order approximation of the NPV change due to the configured
    // shift in each of the par factors c_i
    parSensitivities = element_prod(parSensitivities, parShifts_);
 
    DLOG("Sensitivity conversion done");
 
    return parSensitivities;
}

◆ writeConversionMatrix()

void writeConversionMatrix ( ore::data::Report & reportOut ) const

Write the inverse of the transposed Jacobian to the reportOut.

Definition at line 808 of file parsensitivityanalysis.cpp.

                                                                        {
 
    // Report headers
    report.addColumn("RawFactor(z)", string());
    report.addColumn("ParFactor(c)", string());
    report.addColumn("dz/dc", double(), 12);
 
    // Write report contents i.e. entries where sparse matrix is non-zero
    Size parIdx = 0;
    for (const auto& parKey : parKeys_) {
        Size rawIdx = 0;
        for (const auto& rawKey : rawKeys_) {
            if (!close(jacobi_transp_inv_(parIdx, rawIdx), 0.0)) {
                report.next();
                report.add(to_string(rawKey));
                report.add(to_string(parKey));
                report.add(jacobi_transp_inv_(parIdx, rawIdx));
            }
            rawIdx++;
        }
        parIdx++;
    }
 
    // Close report
    report.end();
}

Here is the call graph for this function:

◆ inverseJacobian()

ParSensitivityAnalysis::ParContainer inverseJacobian ( ) const

Definition at line 206 of file parsensitivityanalysis.hpp.

                                                               { ParSensitivityAnalysis::ParContainer results;
        Size parIdx = 0;
        for (const auto& parKey : parKeys_) {
            Size rawIdx = 0;
            for (const auto& rawKey : rawKeys_) {
                results[{rawKey, parKey}] = jacobi_transp_inv_(parIdx, rawIdx);
                rawIdx++;
            }
            parIdx++;
        }
        return results;
    }