turnbullwakemanasianengine.cpp

/*
 Copyright (C) 2021 Skandinaviska Enskilda Banken AB (publ)
 
 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/exercise.hpp>
#include <ql/pricingengines/asian/turnbullwakemanasianengine.hpp>
#include <ql/pricingengines/blackcalculator.hpp>
 
using namespace QuantLib;
 
void TurnbullWakemanAsianEngine::calculate() const {
 
    // Enforce a few required things
    QL_REQUIRE(arguments_.exercise->type() == Exercise::European, "not a European Option");
    QL_REQUIRE(arguments_.averageType == Average::Type::Arithmetic,
               "must be Arithmetic Average::Type");
 
    // Calculate the accrued portion
    Size pastFixings = arguments_.pastFixings;
    Size futureFixings = arguments_.fixingDates.size();
    Real accruedAverage = 0;
    if (pastFixings != 0) {
        accruedAverage = arguments_.runningAccumulator / (pastFixings + futureFixings);
    }
    results_.additionalResults["accrued"] = accruedAverage;
 
    Real discount = process_->riskFreeRate()->discount(arguments_.exercise->lastDate());
    results_.additionalResults["discount"] = discount;
 
    ext::shared_ptr<PlainVanillaPayoff> payoff =
        ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
    QL_REQUIRE(payoff, "non-plain payoff given");
 
    // We will read the volatility off the surface at the effective strike
    Real effectiveStrike = payoff->strike() - accruedAverage;
    results_.additionalResults["strike"] = payoff->strike();
    results_.additionalResults["effective_strike"] = effectiveStrike;
 
    // If the effective strike is negative, exercise resp. permanent OTM is guaranteed and the
    // valuation is made easy
    Size m = futureFixings + pastFixings;
    if (effectiveStrike <= 0.0) {
        // For a reference, see "Option Pricing Formulas", Haug, 2nd ed, p. 193
        if (payoff->optionType() == Option::Type::Call) {
            Real spot = process_->stateVariable()->value();
            Real S_A_hat = accruedAverage;
            for (const auto& fd : arguments_.fixingDates) {
                S_A_hat += (spot * process_->dividendYield()->discount(fd) /
                            process_->riskFreeRate()->discount(fd)) /
                           m;
            }
            results_.value = discount * (S_A_hat - payoff->strike());
            results_.delta = discount * (S_A_hat - accruedAverage) / spot;
        } else if (payoff->optionType() == Option::Type::Put) {
            results_.value = 0;
            results_.delta = 0;
        }
        results_.gamma = 0;
        return;
    }
 
    // We should only get this far when the effectiveStrike > 0 but will check anyway
    QL_REQUIRE(effectiveStrike > 0.0, "expected effectiveStrike to be positive");
 
    // Expected value of the non-accrued portion of the average prices
    // In general, m will equal n below if there is no accrued. If accrued, m > n.
    Real EA = 0.0;
    std::vector<Real> forwards;
    std::vector<Time> times;
    std::vector<Real> spotVars;
    std::vector<Real> spotVolsVec; // additional results only
    Real spot = process_->stateVariable()->value();
 
    for (const auto& fd : arguments_.fixingDates) {
        DiscountFactor dividendDiscount = process_->dividendYield()->discount(fd);
        DiscountFactor riskFreeDiscountForFwdEstimation = process_->riskFreeRate()->discount(fd);
 
        forwards.push_back(spot * dividendDiscount / riskFreeDiscountForFwdEstimation);
        times.push_back(process_->blackVolatility()->timeFromReference(fd));
 
        spotVars.push_back(
            process_->blackVolatility()->blackVariance(times.back(), effectiveStrike));
        spotVolsVec.push_back(std::sqrt(spotVars.back() / times.back()));
 
        EA += forwards.back();
    }
    EA /= m;
 
    // Expected value of A^2.
    Real EA2 = 0.0;
    Size n = forwards.size();
 
    for (Size i = 0; i < n; ++i) {
        EA2 += forwards[i] * forwards[i] * exp(spotVars[i]);
        for (Size j = 0; j < i; ++j) {
            EA2 += 2 * forwards[i] * forwards[j] * exp(spotVars[j]);
        }
    }
 
    EA2 /= m * m;
 
    // Calculate value
    Real tn = times.back();
    Real sigma = sqrt(log(EA2 / (EA * EA)) / tn);
 
    // Populate results
    BlackCalculator black(payoff->optionType(), effectiveStrike, EA, sigma * sqrt(tn), discount);
 
    results_.value = black.value();
    results_.delta = black.delta(spot);
    results_.gamma = black.gamma(spot);
 
    results_.additionalResults["forward"] = EA;
    results_.additionalResults["exp_A_2"] = EA2;
    results_.additionalResults["tte"] = tn;
    results_.additionalResults["sigma"] = sigma;
    results_.additionalResults["times"] = times;
    results_.additionalResults["spotVols"] = spotVolsVec;
    results_.additionalResults["forwards"] = forwards;
}