1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2004 Neil Firth
5 Copyright (C) 2007 StatPro Italia srl
6 Copyright (C) 2013 Fabien Le Floc'h
7
8 This file is part of QuantLib, a free-software/open-source library
9 for financial quantitative analysts and developers - http://quantlib.org/
10
11 QuantLib is free software: you can redistribute it and/or modify it
12 under the terms of the QuantLib license. You should have received a
13 copy of the license along with this program; if not, please email
14 <quantlib-dev@lists.sf.net>. The license is also available online at
15 <http://quantlib.org/license.shtml>.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the license for more details.
20*/
21
22#include <ql/exercise.hpp>
23#include <ql/math/distributions/normaldistribution.hpp>
24#include <ql/pricingengines/blackcalculator.hpp>
25#include <ql/pricingengines/blackformula.hpp>
26#include <ql/pricingengines/vanilla/baroneadesiwhaleyengine.hpp>
27#include <ql/pricingengines/vanilla/juquadraticengine.hpp>
28#include <utility>
29
30namespace QuantLib {
31
32 /* An Approximate Formula for Pricing American Options
33 Journal of Derivatives Winter 1999
34 Ju, N.
35 */
36
37
38 JuQuadraticApproximationEngine::JuQuadraticApproximationEngine(
39 ext::shared_ptr<GeneralizedBlackScholesProcess> process)
40 : process_(std::move(process)) {
41 registerWith(h: process_);
42 }
43
44 void JuQuadraticApproximationEngine::calculate() const {
45
46 QL_REQUIRE(arguments_.exercise->type() == Exercise::American,
47 "not an American Option");
48
49 ext::shared_ptr<AmericanExercise> ex =
50 ext::dynamic_pointer_cast<AmericanExercise>(r: arguments_.exercise);
51 QL_REQUIRE(ex, "non-American exercise given");
52 QL_REQUIRE(!ex->payoffAtExpiry(),
53 "payoff at expiry not handled");
54
55 ext::shared_ptr<StrikedTypePayoff> payoff =
56 ext::dynamic_pointer_cast<StrikedTypePayoff>(r: arguments_.payoff);
57 QL_REQUIRE(payoff, "non-striked payoff given");
58
59 Real variance = process_->blackVolatility()->blackVariance(
60 d: ex->lastDate(), strike: payoff->strike());
61 DiscountFactor dividendDiscount = process_->dividendYield()->discount(
62 d: ex->lastDate());
63 DiscountFactor riskFreeDiscount = process_->riskFreeRate()->discount(
64 d: ex->lastDate());
65 Real spot = process_->stateVariable()->value();
66 QL_REQUIRE(spot > 0.0, "negative or null underlying given");
67 Real forwardPrice = spot * dividendDiscount / riskFreeDiscount;
68 BlackCalculator black(payoff, forwardPrice,
69 std::sqrt(x: variance), riskFreeDiscount);
70
71 if (dividendDiscount>=1.0 && payoff->optionType()==Option::Call) {
72 // early exercise never optimal
73 results_.value = black.value();
74 results_.delta = black.delta(spot);
75 results_.deltaForward = black.deltaForward();
76 results_.elasticity = black.elasticity(spot);
77 results_.gamma = black.gamma(spot);
78
79 DayCounter rfdc = process_->riskFreeRate()->dayCounter();
80 DayCounter divdc = process_->dividendYield()->dayCounter();
81 DayCounter voldc = process_->blackVolatility()->dayCounter();
82 Time t =
83 rfdc.yearFraction(d1: process_->riskFreeRate()->referenceDate(),
84 d2: arguments_.exercise->lastDate());
85 results_.rho = black.rho(maturity: t);
86
87 t = divdc.yearFraction(d1: process_->dividendYield()->referenceDate(),
88 d2: arguments_.exercise->lastDate());
89 results_.dividendRho = black.dividendRho(maturity: t);
90
91 t = voldc.yearFraction(d1: process_->blackVolatility()->referenceDate(),
92 d2: arguments_.exercise->lastDate());
93 results_.vega = black.vega(maturity: t);
94 results_.theta = black.theta(spot, maturity: t);
95 results_.thetaPerDay = black.thetaPerDay(spot, maturity: t);
96
97 results_.strikeSensitivity = black.strikeSensitivity();
98 results_.itmCashProbability = black.itmCashProbability();
99 } else {
100 // early exercise can be optimal
101 CumulativeNormalDistribution cumNormalDist;
102 NormalDistribution normalDist;
103
104 Real tolerance = 1e-6;
105 Real Sk = BaroneAdesiWhaleyApproximationEngine::criticalPrice(
106 payoff, riskFreeDiscount, dividendDiscount, variance,
107 tolerance);
108
109 Real forwardSk = Sk * dividendDiscount / riskFreeDiscount;
110
111 Real alpha = -2.0*std::log(x: riskFreeDiscount)/(variance);
112 Real beta = 2.0*std::log(x: dividendDiscount/riskFreeDiscount)/
113 (variance);
114 Real h = 1 - riskFreeDiscount;
115 Real phi;
116 switch (payoff->optionType()) {
117 case Option::Call:
118 phi = 1;
119 break;
120 case Option::Put:
121 phi = -1;
122 break;
123 default:
124 QL_FAIL("unknown option type");
125 }
126 //it can throw: to be fixed
127 Real temp_root = std::sqrt (x: (beta-1)*(beta-1) + (4*alpha)/h);
128 Real lambda = (-(beta-1) + phi * temp_root) / 2;
129 Real lambda_prime = - phi * alpha / (h*h * temp_root);
130
131 Real black_Sk = blackFormula(optionType: payoff->optionType(), strike: payoff->strike(),
132 forward: forwardSk, stdDev: std::sqrt(x: variance)) * riskFreeDiscount;
133 Real hA = phi * (Sk - payoff->strike()) - black_Sk;
134
135 Real d1_Sk = (std::log(x: forwardSk/payoff->strike()) + 0.5*variance)
136 /std::sqrt(x: variance);
137 Real d2_Sk = d1_Sk - std::sqrt(x: variance);
138 Real part1 = forwardSk * normalDist(d1_Sk) /
139 (alpha * std::sqrt(x: variance));
140 Real part2 = - phi * forwardSk * cumNormalDist(phi * d1_Sk) *
141 std::log(x: dividendDiscount) / std::log(x: riskFreeDiscount);
142 Real part3 = + phi * payoff->strike() * cumNormalDist(phi * d2_Sk);
143 Real V_E_h = part1 + part2 + part3;
144
145 Real b = (1-h) * alpha * lambda_prime / (2*(2*lambda + beta - 1));
146 Real c = - ((1 - h) * alpha / (2 * lambda + beta - 1)) *
147 (V_E_h / (hA) + 1 / h + lambda_prime / (2*lambda + beta - 1));
148 Real temp_spot_ratio = std::log(x: spot / Sk);
149 Real chi = temp_spot_ratio * (b * temp_spot_ratio + c);
150
151 if (phi*(Sk-spot) > 0) {
152 results_.value = black.value() +
153 hA * std::pow(x: (spot/Sk), y: lambda) / (1 - chi);
154 Real temp_chi_prime = (2 * b / spot) * std::log(x: spot/Sk);
155 Real chi_prime = temp_chi_prime + c / spot;
156 Real chi_double_prime = 2*b/(spot*spot)
157 - temp_chi_prime / spot - c / (spot*spot);
158 Real d1_S = (std::log(x: forwardPrice/payoff->strike()) + 0.5*variance)
159 / std::sqrt(x: variance);
160 //There is a typo in the original paper from Ju-Zhong
161 //the first term is the Black-Scholes delta/gamma.
162 results_.delta = phi * dividendDiscount * cumNormalDist (phi * d1_S)
163 + (lambda / (spot * (1 - chi)) + chi_prime / ((1 - chi)*(1 - chi))) *
164 (phi * (Sk - payoff->strike()) - black_Sk) * std::pow(x: (spot/Sk), y: lambda);
165
166 results_.gamma = dividendDiscount * normalDist (phi*d1_S)
167 / (spot * std::sqrt(x: variance))
168 + (2 * lambda * chi_prime / (spot * (1 - chi) * (1 - chi))
169 + 2 * chi_prime * chi_prime / ((1 - chi) * (1 - chi) * (1 - chi))
170 + chi_double_prime / ((1 - chi) * (1 - chi))
171 + lambda * (lambda - 1) / (spot * spot * (1 - chi)))
172 * (phi * (Sk - payoff->strike()) - black_Sk)
173 * std::pow(x: (spot/Sk), y: lambda);
174 } else {
175 results_.value = phi * (spot - payoff->strike());
176 results_.delta = phi;
177 results_.gamma = 0;
178 }
179
180 } // end of "early exercise can be optimal"
181 }
182
183}
184

source code of quantlib/ql/pricingengines/vanilla/juquadraticengine.cpp