juquadraticengine.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2004 Neil Firth
 Copyright (C) 2007 StatPro Italia srl
 Copyright (C) 2013 Fabien Le Floc'h
 
 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/math/distributions/normaldistribution.hpp>
#include <ql/pricingengines/blackcalculator.hpp>
#include <ql/pricingengines/blackformula.hpp>
#include <ql/pricingengines/vanilla/baroneadesiwhaleyengine.hpp>
#include <ql/pricingengines/vanilla/juquadraticengine.hpp>
#include <utility>
 
namespace QuantLib {
 
    /*  An Approximate Formula for Pricing American Options
        Journal of Derivatives Winter 1999
        Ju, N.
    */
 
 
    JuQuadraticApproximationEngine::JuQuadraticApproximationEngine(
        ext::shared_ptr<GeneralizedBlackScholesProcess> process)
    : process_(std::move(process)) {
        registerWith(process_);
    }
 
    void JuQuadraticApproximationEngine::calculate() const {
 
        QL_REQUIRE(arguments_.exercise->type() == Exercise::American,
                   "not an American Option");
 
        ext::shared_ptr<AmericanExercise> ex =
            ext::dynamic_pointer_cast<AmericanExercise>(arguments_.exercise);
        QL_REQUIRE(ex, "non-American exercise given");
        QL_REQUIRE(!ex->payoffAtExpiry(),
                   "payoff at expiry not handled");
 
        ext::shared_ptr<StrikedTypePayoff> payoff =
            ext::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
        QL_REQUIRE(payoff, "non-striked payoff given");
 
        Real variance = process_->blackVolatility()->blackVariance(
            ex->lastDate(), payoff->strike());
        DiscountFactor dividendDiscount = process_->dividendYield()->discount(
            ex->lastDate());
        DiscountFactor riskFreeDiscount = process_->riskFreeRate()->discount(
            ex->lastDate());
        Real spot = process_->stateVariable()->value();
        QL_REQUIRE(spot > 0.0, "negative or null underlying given");
        Real forwardPrice = spot * dividendDiscount / riskFreeDiscount;
        BlackCalculator black(payoff, forwardPrice,
                              std::sqrt(variance), riskFreeDiscount);
 
        if (dividendDiscount>=1.0 && payoff->optionType()==Option::Call) {
            // early exercise never optimal
            results_.value        = black.value();
            results_.delta        = black.delta(spot);
            results_.deltaForward = black.deltaForward();
            results_.elasticity   = black.elasticity(spot);
            results_.gamma        = black.gamma(spot);
 
            DayCounter rfdc  = process_->riskFreeRate()->dayCounter();
            DayCounter divdc = process_->dividendYield()->dayCounter();
            DayCounter voldc = process_->blackVolatility()->dayCounter();
            Time t =
                rfdc.yearFraction(process_->riskFreeRate()->referenceDate(),
                                  arguments_.exercise->lastDate());
            results_.rho = black.rho(t);
 
            t = divdc.yearFraction(process_->dividendYield()->referenceDate(),
                                   arguments_.exercise->lastDate());
            results_.dividendRho = black.dividendRho(t);
 
            t = voldc.yearFraction(process_->blackVolatility()->referenceDate(),
                                   arguments_.exercise->lastDate());
            results_.vega        = black.vega(t);
            results_.theta       = black.theta(spot, t);
            results_.thetaPerDay = black.thetaPerDay(spot, t);
 
            results_.strikeSensitivity  = black.strikeSensitivity();
            results_.itmCashProbability = black.itmCashProbability();
        } else {
            // early exercise can be optimal
            CumulativeNormalDistribution cumNormalDist;
            NormalDistribution normalDist;
 
            Real tolerance = 1e-6;
            Real Sk = BaroneAdesiWhaleyApproximationEngine::criticalPrice(
                payoff, riskFreeDiscount, dividendDiscount, variance,
                tolerance);
 
            Real forwardSk = Sk * dividendDiscount / riskFreeDiscount;
 
            Real alpha = -2.0*std::log(riskFreeDiscount)/(variance);
            Real beta = 2.0*std::log(dividendDiscount/riskFreeDiscount)/
                                                (variance);
            Real h = 1 - riskFreeDiscount;
            Real phi;
            switch (payoff->optionType()) {
                case Option::Call:
                    phi = 1;
                    break;
                case Option::Put:
                    phi = -1;
                    break;
                default:
                  QL_FAIL("unknown option type");
            }
            //it can throw: to be fixed
            Real temp_root = std::sqrt ((beta-1)*(beta-1) + (4*alpha)/h);
            Real lambda = (-(beta-1) + phi * temp_root) / 2;
            Real lambda_prime = - phi * alpha / (h*h * temp_root);
 
            Real black_Sk = blackFormula(payoff->optionType(), payoff->strike(),
                                         forwardSk, std::sqrt(variance)) * riskFreeDiscount;
            Real hA = phi * (Sk - payoff->strike()) - black_Sk;
 
            Real d1_Sk = (std::log(forwardSk/payoff->strike()) + 0.5*variance)
                /std::sqrt(variance);
            Real d2_Sk = d1_Sk - std::sqrt(variance);
            Real part1 = forwardSk * normalDist(d1_Sk) /
                                        (alpha * std::sqrt(variance));
            Real part2 = - phi * forwardSk * cumNormalDist(phi * d1_Sk) *
                      std::log(dividendDiscount) / std::log(riskFreeDiscount);
            Real part3 = + phi * payoff->strike() * cumNormalDist(phi * d2_Sk);
            Real V_E_h = part1 + part2 + part3;
 
            Real b = (1-h) * alpha * lambda_prime / (2*(2*lambda + beta - 1));
            Real c = - ((1 - h) * alpha / (2 * lambda + beta - 1)) *
                (V_E_h / (hA) + 1 / h + lambda_prime / (2*lambda + beta - 1));
            Real temp_spot_ratio = std::log(spot / Sk);
            Real chi = temp_spot_ratio * (b * temp_spot_ratio + c);
 
            if (phi*(Sk-spot) > 0) {
                results_.value = black.value() +
                    hA * std::pow((spot/Sk), lambda) / (1 - chi);
                Real temp_chi_prime = (2 * b / spot) * std::log(spot/Sk);
                Real chi_prime = temp_chi_prime + c / spot;
                Real chi_double_prime = 2*b/(spot*spot)
                    - temp_chi_prime / spot - c / (spot*spot);
                Real d1_S = (std::log(forwardPrice/payoff->strike()) + 0.5*variance)
                    / std::sqrt(variance);
                //There is a typo in the original paper from Ju-Zhong
                //the first term is the Black-Scholes delta/gamma.    
                results_.delta = phi * dividendDiscount * cumNormalDist (phi * d1_S)
                    + (lambda / (spot * (1 - chi)) + chi_prime / ((1 - chi)*(1 - chi))) *
                    (phi * (Sk - payoff->strike()) - black_Sk) * std::pow((spot/Sk), lambda);
 
                results_.gamma = dividendDiscount * normalDist (phi*d1_S) 
                    / (spot * std::sqrt(variance))
                    + (2 * lambda * chi_prime / (spot * (1 - chi) * (1 - chi))
                        + 2 * chi_prime * chi_prime / ((1 - chi) * (1 - chi) * (1 - chi))
                        + chi_double_prime / ((1 - chi) * (1 - chi))
                        + lambda * (lambda - 1) / (spot * spot * (1 - chi)))
                    * (phi * (Sk - payoff->strike()) - black_Sk)
                    * std::pow((spot/Sk), lambda);
            } else {
                results_.value = phi * (spot - payoff->strike());
                results_.delta = phi;
                results_.gamma = 0;
            }
 
        } // end of "early exercise can be optimal"
    }
 
}