suowangdoublebarrierengine.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2013 Yue Tian
 
 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/experimental/barrieroption/suowangdoublebarrierengine.hpp>
#include <ql/instruments/europeanoption.hpp>
#include <ql/pricingengines/vanilla/analyticeuropeanengine.hpp>
#include <utility>
 
namespace QuantLib {
 
    SuoWangDoubleBarrierEngine::SuoWangDoubleBarrierEngine(
        ext::shared_ptr<GeneralizedBlackScholesProcess> process, int series)
    : process_(std::move(process)), series_(series) {
        registerWith(process_);
    }
 
    void SuoWangDoubleBarrierEngine::calculate() const {
 
        ext::shared_ptr<PlainVanillaPayoff> payoff =
            ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
        QL_REQUIRE(payoff, "non-plain payoff given");
        QL_REQUIRE(payoff->strike()>0.0,
                   "strike must be positive");
 
        Real K = payoff->strike();
        Real S = process_->x0();
        QL_REQUIRE(S > 0.0, "negative or null underlying given");
        QL_REQUIRE(!triggered(S), "barrier touched");
 
        DoubleBarrier::Type barrierType = arguments_.barrierType;
        QL_REQUIRE(barrierType == DoubleBarrier::KnockOut || 
                   barrierType == DoubleBarrier::KnockIn,
                   "only KnockIn and KnockOut options supported");
 
        Real L = arguments_.barrier_lo;
        Real H = arguments_.barrier_hi;
        Real K_up = std::min(H, K);
        Real K_down = std::max(L, K);
        Time T = residualTime();
        Real rd = riskFreeRate();
        Real dd = riskFreeDiscount();
        Real rf = dividendYield();
        Real df = dividendDiscount();
        Real vol = volatility();
        Real mu = rd - rf - vol*vol/2.0;
        Real sgn = mu > 0 ? 1.0 :(mu < 0 ? -1.0: 0.0);
        //rebate
        Real R_L = arguments_.rebate;
        Real R_H = arguments_.rebate;
 
        //european option
        EuropeanOption europeanOption(payoff, arguments_.exercise);
        ext::shared_ptr<PricingEngine> analyticEuropeanEngine =
            ext::make_shared<AnalyticEuropeanEngine>(process_);
        europeanOption.setPricingEngine(analyticEuropeanEngine);
        Real european = europeanOption.NPV();
 
        Real barrierOut = 0;
        Real rebateIn = 0;
        for(int n = -series_; n < series_; n++){
            Real d1 = D(S/H*std::pow(L/H, 2.0*n), vol*vol+mu, vol, T);
            Real d2 = d1 - vol*std::sqrt(T);
            Real g1 = D(H/S*std::pow(L/H, 2.0*n - 1.0), vol*vol+mu, vol, T);
            Real g2 = g1 - vol*std::sqrt(T);
            Real h1 = D(S/H*std::pow(L/H, 2.0*n - 1.0), vol*vol+mu, vol, T);
            Real h2 = h1 - vol*std::sqrt(T);
            Real k1 = D(L/S*std::pow(L/H, 2.0*n - 1.0), vol*vol+mu, vol, T);
            Real k2 = k1 - vol*std::sqrt(T);
            Real d1_down = D(S/K_down*std::pow(L/H, 2.0*n), vol*vol+mu, vol, T);
            Real d2_down = d1_down - vol*std::sqrt(T);
            Real d1_up = D(S/K_up*std::pow(L/H, 2.0*n), vol*vol+mu, vol, T);
            Real d2_up = d1_up - vol*std::sqrt(T);
            Real k1_down = D((H*H)/(K_down*S)*std::pow(L/H, 2.0*n), vol*vol+mu, vol, T);
            Real k2_down = k1_down - vol*std::sqrt(T);
            Real k1_up = D((H*H)/(K_up*S)*std::pow(L/H, 2.0*n), vol*vol+mu, vol, T);
            Real k2_up = k1_up - vol*std::sqrt(T);
 
            if( payoff->optionType() == Option::Call) {
                barrierOut += std::pow(L/H, 2.0 * n * mu/(vol*vol))*
                            (df*S*std::pow(L/H, 2.0*n)*(f_(d1_down)-f_(d1))
                            -dd*K*(f_(d2_down)-f_(d2))
                            -df*std::pow(L/H, 2.0*n)*H*H/S*std::pow(H/S, 2.0*mu/(vol*vol))*(f_(k1_down)-f_(k1))
                            +dd*K*std::pow(H/S,2.0*mu/(vol*vol))*(f_(k2_down)-f_(k2)));
            }
            else if(payoff->optionType() == Option::Put){
                barrierOut += std::pow(L/H, 2.0 * n * mu/(vol*vol))*
                            (dd*K*(f_(h2)-f_(d2_up))
                            -df*S*std::pow(L/H, 2.0*n)*(f_(h1)-f_(d1_up))
                            -dd*K*std::pow(H/S,2.0*mu/(vol*vol))*(f_(g2)-f_(k2_up))
                            +df*std::pow(L/H, 2.0*n)*H*H/S*std::pow(H/S, 2.0*mu/(vol*vol))*(f_(g1)-f_(k1_up)));
            }
            else {
                QL_FAIL("option type not recognized");
            }
 
            Real v1 = D(H/S*std::pow(H/L, 2.0*n), -mu, vol, T);
            Real v2 = D(H/S*std::pow(H/L, 2.0*n), mu, vol, T);
            Real v3 = D(S/L*std::pow(H/L, 2.0*n), -mu, vol, T);
            Real v4 = D(S/L*std::pow(H/L, 2.0*n), mu, vol, T);
            rebateIn +=  dd * R_H * sgn * (std::pow(L/H, 2.0*n*mu/(vol*vol)) * f_(sgn * v1) - std::pow(H/S, 2.0*mu/(vol*vol)) * f_(-sgn * v2))
                       + dd * R_L * sgn * (std::pow(L/S, 2.0*mu/(vol*vol)) * f_(-sgn * v3) - std::pow(H/L, 2.0*n*mu/(vol*vol)) * f_(sgn * v4));
        }
 
        //rebate paid at maturity
        if(barrierType == DoubleBarrier::KnockOut)
            results_.value = barrierOut ;
        else
            results_.value = european - barrierOut;
        results_.additionalResults["vanilla"] = european;
        results_.additionalResults["barrierOut"] = barrierOut;
        results_.additionalResults["barrierIn"] = Real(european - barrierOut);
        results_.additionalResults["rebateIn"] = rebateIn;
    }
 
 
    Real SuoWangDoubleBarrierEngine::strike() const {
        ext::shared_ptr<PlainVanillaPayoff> payoff =
            ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
        QL_REQUIRE(payoff, "non-plain payoff given");
        return payoff->strike();
    }
 
    Time SuoWangDoubleBarrierEngine::residualTime() const {
        return process_->time(arguments_.exercise->lastDate());
    }
 
    Volatility SuoWangDoubleBarrierEngine::volatility() const {
        return process_->blackVolatility()->blackVol(residualTime(), strike());
    }
 
    Rate SuoWangDoubleBarrierEngine::riskFreeRate() const {
        return process_->riskFreeRate()->zeroRate(residualTime(), Continuous,
                                                  NoFrequency);
    }
 
    DiscountFactor SuoWangDoubleBarrierEngine::riskFreeDiscount() const {
        return process_->riskFreeRate()->discount(residualTime());
    }
 
    Rate SuoWangDoubleBarrierEngine::dividendYield() const {
        return process_->dividendYield()->zeroRate(residualTime(),
                                                   Continuous, NoFrequency);
    }
 
    DiscountFactor SuoWangDoubleBarrierEngine::dividendDiscount() const {
        return process_->dividendYield()->discount(residualTime());
    }
 
    Real SuoWangDoubleBarrierEngine::D(Real X, Real lambda, Real sigma, Real T) const {
        return (std::log(X) + lambda * T)/(sigma * std::sqrt(T));
    }
 
}