QuantExt::CommodityAveragePriceOptionMomementMatching Namespace Reference

Classes
struct	MomentMatchingResults

Functions
MomentMatchingResults	matchFirstTwoMomentsTurnbullWakeman (const ext::shared_ptr< CommodityIndexedAverageCashFlow > &flow, const ext::shared_ptr< QuantLib::BlackVolTermStructure > &vol, const std::function< double(const QuantLib::Date &expiry1, const QuantLib::Date &expiry2)> &rho, QuantLib::Real strike)

Function Documentation

matchFirstTwoMomentsTurnbullWakeman()

MomentMatchingResults matchFirstTwoMomentsTurnbullWakeman	(	const ext::shared_ptr< CommodityIndexedAverageCashFlow > &	flow,
		const ext::shared_ptr< QuantLib::BlackVolTermStructure > &	vol,
		const std::function< double(const QuantLib::Date &expiry1, const QuantLib::Date &expiry2)> &	rho,
		QuantLib::Real	strike
	)

Definition at line 44 of file commodityapoengine.cpp.

                         {
    Date today = Settings::instance().evaluationDate();
    MomentMatchingResults res;
 
    res.tn = 0.0;
    res.accruals = 0.0;
    double EA = 0;
    std::vector<Date> futureExpiries;
    std::map<Date, Real> futureVols;
    std::vector<double> spotVariances;
    size_t n = flow->indices().size();
    double atmUnderlyingCcy = 0;
 
    for (const auto& [pricingDate, index] : flow->indices()) {
        Date fixingDate = index->fixingCalendar().adjust(pricingDate, Preceding);
        Real fxRate = (flow->fxIndex()) ? flow->fxIndex()->fixing(fixingDate) : 1.0;
        res.indexNames.push_back(index->name());
        res.pricingDates.push_back(fixingDate);
        res.indexExpiries.push_back(index->expiryDate());
        res.fixings.push_back(index->fixing(fixingDate) * fxRate);
        if (pricingDate <= today) {
            res.accruals += res.fixings.back();
        } else {
            atmUnderlyingCcy = index->fixing(fixingDate);
            res.forwards.push_back(res.fixings.back());
            res.times.push_back(vol->timeFromReference(pricingDate));
            // use ATM vol if no strike is given
            double K = strike == Null<Real>() ? atmUnderlyingCcy : strike;
            if (flow->useFuturePrice()) {
                Date expiry = index->expiryDate();
                futureExpiries.push_back(expiry);
                if (futureVols.count(expiry) == 0) {
                    futureVols[expiry] = vol->blackVol(expiry, K);
                }
            } else {
                spotVariances.push_back(
                    vol->blackVariance(res.times.back(), K));
                res.spotVols.push_back(std::sqrt(spotVariances.back() / res.times.back()));
            }
            EA += res.forwards.back();
        }
    }
    res.accruals /= static_cast<double>(n);
    EA /= static_cast<double>(n);
 
    res.forward = EA;
 
    double EA2 = 0.0;
 
    if (flow->useFuturePrice()) {
        // References future prices
        for (Size i = 0; i < res.forwards.size(); ++i) {
            Date e_i = futureExpiries[i];
            Volatility v_i = futureVols.at(e_i);
            res.futureVols.push_back(v_i);
            EA2 += res.forwards[i] * res.forwards[i] * exp(v_i * v_i * res.times[i]);
            for (Size j = 0; j < i; ++j) {
                Date e_j = futureExpiries[j];
                Volatility v_j = futureVols.at(e_j);
                EA2 += 2 * res.forwards[i] * res.forwards[j] * exp(rho(e_i, e_j) * v_i * v_j * res.times[j]);
            }
        }
    } else {
        // References spot prices
        for (Size i = 0; i < res.forwards.size(); ++i) {
            EA2 += res.forwards[i] * res.forwards[i] * exp(spotVariances[i]);
            for (Size j = 0; j < i; ++j) {
                EA2 += 2 * res.forwards[i] * res.forwards[j] * exp(spotVariances[j]);
            }
        }
    }
    EA2 /= std::pow(static_cast<double>(n), 2.0);
    res.EA2 = EA2;
 
    QL_REQUIRE(!std::isinf(EA2),
               "moment matching fails (EA2 = inf) - this is possibly due to too high input volatilities.");
 
    // Calculate value
    if (!res.times.empty()) {
        res.tn = res.times.back();
        double s = EA2 / (EA * EA);
        // if future vol = 0 for all dates, then EA2 = EA*EA, but
        // due to numerical precision EA2 can be actually less than EA*EA 
        if (s < 1.0 || QuantLib::close_enough(s, 1)) {
            res.sigma = 0.0;
        } else {
            res.sigma = std::sqrt(std::log(s) / res.tn);
        }
    } else {
        res.tn = 0;
        res.sigma = 0;
    }
    return res;
}

Here is the caller graph for this function: