/src/quantlib/ql/pricingengines/basket/stulzengine.cpp

Source
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2004 Ferdinando Ametrano
 Copyright (C) 2004 Neil Firth
 Copyright (C) 2007 StatPro Italia srl

 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
 <https://www.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/bivariatenormaldistribution.hpp>
#include <ql/math/distributions/normaldistribution.hpp>
#include <ql/pricingengines/basket/stulzengine.hpp>
#include <ql/pricingengines/blackcalculator.hpp>
#include <ql/pricingengines/blackformula.hpp>
#include <utility>

namespace QuantLib {

    namespace {

        // calculate the value of euro min basket call
        Real euroTwoAssetMinBasketCall(Real forward1, Real forward2,
                                       Real strike,
                                       DiscountFactor riskFreeDiscount,
                                       Real variance1, Real variance2,
                                       Real rho) {

            Real stdDev1 = std::sqrt(variance1);
            Real stdDev2 = std::sqrt(variance2);

            Real variance = variance1 + variance2 - 2*rho*stdDev1*stdDev2;
            Real stdDev = std::sqrt(variance);

            Real modRho1 = (rho * stdDev2 - stdDev1) / stdDev;
            Real modRho2 = (rho * stdDev1 - stdDev2) / stdDev;

            Real D1 = (std::log(forward1/forward2) + 0.5*variance) / stdDev;

            Real alfa, beta, gamma;
            if (strike != 0.0) {
                auto bivCNorm =
                    BivariateCumulativeNormalDistribution(rho);
                auto bivCNormMod2 =
                    BivariateCumulativeNormalDistribution(modRho2);
                auto bivCNormMod1 =
                    BivariateCumulativeNormalDistribution(modRho1);

                Real D1_1 =
                    (std::log(forward1/strike) + 0.5*variance1) / stdDev1;
                Real D1_2 =
                    (std::log(forward2/strike) + 0.5*variance2) / stdDev2;
                alfa = bivCNormMod1(D1_1, -D1);
                beta = bivCNormMod2(D1_2, D1 - stdDev);
                gamma = bivCNorm(D1_1 - stdDev1, D1_2 - stdDev2);
            } else {
                CumulativeNormalDistribution cum;
                alfa = cum(-D1);
                beta = cum(D1 - stdDev);
                gamma = 1.0;
            }

            return riskFreeDiscount *
                (forward1*alfa + forward2*beta - strike*gamma);

        }

        // calculate the value of euro max basket call
        Real euroTwoAssetMaxBasketCall(Real forward1, Real forward2,
                                       Real strike,
                                       DiscountFactor riskFreeDiscount,
                                       Real variance1, Real variance2,
                                       Real rho) {

            auto payoff = ext::make_shared<PlainVanillaPayoff>(Option::Call, strike);

            Real black1 = blackFormula(payoff->optionType(), payoff->strike(),
                forward1, std::sqrt(variance1)) * riskFreeDiscount;

            Real black2 = blackFormula(payoff->optionType(), payoff->strike(),
                forward2, std::sqrt(variance2)) * riskFreeDiscount;

            return black1 + black2 -
                euroTwoAssetMinBasketCall(forward1, forward2, strike,
                                          riskFreeDiscount,
                                          variance1, variance2, rho);
        }
    }

    StulzEngine::StulzEngine(ext::shared_ptr<GeneralizedBlackScholesProcess> process1,
                             ext::shared_ptr<GeneralizedBlackScholesProcess> process2,
                             Real correlation)
    : process1_(std::move(process1)), process2_(std::move(process2)), rho_(correlation) {
        registerWith(process1_);
        registerWith(process2_);
    }

    void StulzEngine::calculate() const {

        QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
                   "not an European Option");

        ext::shared_ptr<EuropeanExercise> exercise =
            ext::dynamic_pointer_cast<EuropeanExercise>(arguments_.exercise);
        QL_REQUIRE(exercise, "not an European Option");

        ext::shared_ptr<BasketPayoff> basket_payoff =
            ext::dynamic_pointer_cast<BasketPayoff>(arguments_.payoff);

        ext::shared_ptr<MinBasketPayoff> min_basket =
            ext::dynamic_pointer_cast<MinBasketPayoff>(arguments_.payoff);

        ext::shared_ptr<MaxBasketPayoff> max_basket =
            ext::dynamic_pointer_cast<MaxBasketPayoff>(arguments_.payoff);
        QL_REQUIRE(min_basket || max_basket, "unknown basket type");

        ext::shared_ptr<PlainVanillaPayoff> payoff =
            ext::dynamic_pointer_cast<PlainVanillaPayoff>(basket_payoff->basePayoff());
        QL_REQUIRE(payoff, "non-plain payoff given");

        Real strike = payoff->strike();

        Real variance1 = process1_->blackVolatility()->blackVariance(
                                                exercise->lastDate(), strike);
        Real variance2 = process2_->blackVolatility()->blackVariance(
                                                exercise->lastDate(), strike);

        DiscountFactor riskFreeDiscount =
            process1_->riskFreeRate()->discount(exercise->lastDate());

        // cannot handle non zero dividends, so don't believe this...
        DiscountFactor dividendDiscount1 =
            process1_->dividendYield()->discount(exercise->lastDate());
        DiscountFactor dividendDiscount2 =
            process2_->dividendYield()->discount(exercise->lastDate());

        Real forward1 = process1_->stateVariable()->value() *
            dividendDiscount1 / riskFreeDiscount;
        Real forward2 = process2_->stateVariable()->value() *
            dividendDiscount2 / riskFreeDiscount;

        if (max_basket != nullptr) {
            switch (payoff->optionType()) {
              // euro call on a two asset max basket
              case Option::Call:
                results_.value =
                    euroTwoAssetMaxBasketCall(forward1, forward2, strike,
                                              riskFreeDiscount,
                                              variance1, variance2,
                                              rho_);

                break;
              // euro put on a two asset max basket
              case Option::Put:
                results_.value = strike * riskFreeDiscount -
                    euroTwoAssetMaxBasketCall(forward1, forward2, 0.0,
                                              riskFreeDiscount,
                                              variance1, variance2, rho_) +
                    euroTwoAssetMaxBasketCall(forward1, forward2, strike,
                                              riskFreeDiscount,
                                              variance1, variance2, rho_);
                break;
              default:
                QL_FAIL("unknown option type");
            }
        } else if (min_basket != nullptr) {
            switch (payoff->optionType()) {
              // euro call on a two asset min basket
              case Option::Call:
                results_.value =
                    euroTwoAssetMinBasketCall(forward1, forward2, strike,
                                              riskFreeDiscount,
                                              variance1, variance2,
                                              rho_);
                break;
              // euro put on a two asset min basket
              case Option::Put:
                results_.value = strike * riskFreeDiscount -
                    euroTwoAssetMinBasketCall(forward1, forward2, 0.0,
                                              riskFreeDiscount,
                                              variance1, variance2, rho_) +
                    euroTwoAssetMinBasketCall(forward1, forward2, strike,
                                              riskFreeDiscount,
                                              variance1, variance2, rho_);
                break;
              default:
                QL_FAIL("unknown option type");
            }
        } else {
            QL_FAIL("unknown type");
        }
    }

}

Coverage Report

Created: 2026-08-14 07:10

Line	Count	Source
1		/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3		/*
4		Copyright (C) 2004 Ferdinando Ametrano
5		Copyright (C) 2004 Neil Firth
6		Copyright (C) 2007 StatPro Italia srl
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		<https://www.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/bivariatenormaldistribution.hpp>
24		#include <ql/math/distributions/normaldistribution.hpp>
25		#include <ql/pricingengines/basket/stulzengine.hpp>
26		#include <ql/pricingengines/blackcalculator.hpp>
27		#include <ql/pricingengines/blackformula.hpp>
28		#include <utility>
29
30		namespace QuantLib {
31
32		namespace {
33
34		// calculate the value of euro min basket call
35		Real euroTwoAssetMinBasketCall(Real forward1, Real forward2,
36		Real strike,
37		DiscountFactor riskFreeDiscount,
38		Real variance1, Real variance2,
39	634	Real rho) {
40
41	634	Real stdDev1 = std::sqrt(variance1);
42	634	Real stdDev2 = std::sqrt(variance2);
43
44	634	Real variance = variance1 + variance2 - 2rhostdDev1*stdDev2;
45	634	Real stdDev = std::sqrt(variance);
46
47	634	Real modRho1 = (rho * stdDev2 - stdDev1) / stdDev;
48	634	Real modRho2 = (rho * stdDev1 - stdDev2) / stdDev;
49
50	634	Real D1 = (std::log(forward1/forward2) + 0.5*variance) / stdDev;
51
52	634	Real alfa, beta, gamma;
53	634	if (strike != 0.0) {
54	385	auto bivCNorm =
55	385	BivariateCumulativeNormalDistribution(rho);
56	385	auto bivCNormMod2 =
57	385	BivariateCumulativeNormalDistribution(modRho2);
58	385	auto bivCNormMod1 =
59	385	BivariateCumulativeNormalDistribution(modRho1);
60
61	385	Real D1_1 =
62	385	(std::log(forward1/strike) + 0.5*variance1) / stdDev1;
63	385	Real D1_2 =
64	385	(std::log(forward2/strike) + 0.5*variance2) / stdDev2;
65	385	alfa = bivCNormMod1(D1_1, -D1);
66	385	beta = bivCNormMod2(D1_2, D1 - stdDev);
67	385	gamma = bivCNorm(D1_1 - stdDev1, D1_2 - stdDev2);
68	385	} else {
69	249	CumulativeNormalDistribution cum;
70	249	alfa = cum(-D1);
71	249	beta = cum(D1 - stdDev);
72	249	gamma = 1.0;
73	249	}
74
75	634	return riskFreeDiscount *
76	634	(forward1alfa + forward2beta - strike*gamma);
77
78	634	}
79
80		// calculate the value of euro max basket call
81		Real euroTwoAssetMaxBasketCall(Real forward1, Real forward2,
82		Real strike,
83		DiscountFactor riskFreeDiscount,
84		Real variance1, Real variance2,
85	506	Real rho) {
86
87	506	auto payoff = ext::make_shared<PlainVanillaPayoff>(Option::Call, strike);
88
89	506	Real black1 = blackFormula(payoff->optionType(), payoff->strike(),
90	506	forward1, std::sqrt(variance1)) * riskFreeDiscount;
91
92	506	Real black2 = blackFormula(payoff->optionType(), payoff->strike(),
93	506	forward2, std::sqrt(variance2)) * riskFreeDiscount;
94
95	506	return black1 + black2 -
96	506	euroTwoAssetMinBasketCall(forward1, forward2, strike,
97	506	riskFreeDiscount,
98	506	variance1, variance2, rho);
99	506	}
100		}
101
102		StulzEngine::StulzEngine(ext::shared_ptr<GeneralizedBlackScholesProcess> process1,
103		ext::shared_ptr<GeneralizedBlackScholesProcess> process2,
104		Real correlation)
105	408	: process1_(std::move(process1)), process2_(std::move(process2)), rho_(correlation) {
106	408	registerWith(process1_);
107	408	registerWith(process2_);
108	408	}
109
110	408	void StulzEngine::calculate() const {
111
112	408	QL_REQUIRE(arguments_.exercise->type() == Exercise::European,
113	408	"not an European Option");
114
115	408	ext::shared_ptr<EuropeanExercise> exercise =
116	408	ext::dynamic_pointer_cast<EuropeanExercise>(arguments_.exercise);
117	408	QL_REQUIRE(exercise, "not an European Option");
118
119	408	ext::shared_ptr<BasketPayoff> basket_payoff =
120	408	ext::dynamic_pointer_cast<BasketPayoff>(arguments_.payoff);
121
122	408	ext::shared_ptr<MinBasketPayoff> min_basket =
123	408	ext::dynamic_pointer_cast<MinBasketPayoff>(arguments_.payoff);
124
125	408	ext::shared_ptr<MaxBasketPayoff> max_basket =
126	408	ext::dynamic_pointer_cast<MaxBasketPayoff>(arguments_.payoff);
127	408	QL_REQUIRE(min_basket \|\| max_basket, "unknown basket type");
128
129	385	ext::shared_ptr<PlainVanillaPayoff> payoff =
130	385	ext::dynamic_pointer_cast<PlainVanillaPayoff>(basket_payoff->basePayoff());
131	385	QL_REQUIRE(payoff, "non-plain payoff given");
132
133	385	Real strike = payoff->strike();
134
135	385	Real variance1 = process1_->blackVolatility()->blackVariance(
136	385	exercise->lastDate(), strike);
137	385	Real variance2 = process2_->blackVolatility()->blackVariance(
138	385	exercise->lastDate(), strike);
139
140	385	DiscountFactor riskFreeDiscount =
141	385	process1_->riskFreeRate()->discount(exercise->lastDate());
142
143		// cannot handle non zero dividends, so don't believe this...
144	385	DiscountFactor dividendDiscount1 =
145	385	process1_->dividendYield()->discount(exercise->lastDate());
146	385	DiscountFactor dividendDiscount2 =
147	385	process2_->dividendYield()->discount(exercise->lastDate());
148
149	385	Real forward1 = process1_->stateVariable()->value() *
150	385	dividendDiscount1 / riskFreeDiscount;
151	385	Real forward2 = process2_->stateVariable()->value() *
152	385	dividendDiscount2 / riskFreeDiscount;
153
154	385	if (max_basket != nullptr) {
155	307	switch (payoff->optionType()) {
156		// euro call on a two asset max basket
157	108	case Option::Call:
158	108	results_.value =
159	108	euroTwoAssetMaxBasketCall(forward1, forward2, strike,
160	108	riskFreeDiscount,
161	108	variance1, variance2,
162	108	rho_);
163
164	108	break;
165		// euro put on a two asset max basket
166	199	case Option::Put:
167	199	results_.value = strike * riskFreeDiscount -
168	199	euroTwoAssetMaxBasketCall(forward1, forward2, 0.0,
169	199	riskFreeDiscount,
170	199	variance1, variance2, rho_) +
171	199	euroTwoAssetMaxBasketCall(forward1, forward2, strike,
172	199	riskFreeDiscount,
173	199	variance1, variance2, rho_);
174	199	break;
175	0	default:
176	0	QL_FAIL("unknown option type");
177	307	}
178	307	} else if (min_basket != nullptr) {
179	78	switch (payoff->optionType()) {
180		// euro call on a two asset min basket
181	28	case Option::Call:
182	28	results_.value =
183	28	euroTwoAssetMinBasketCall(forward1, forward2, strike,
184	28	riskFreeDiscount,
185	28	variance1, variance2,
186	28	rho_);
187	28	break;
188		// euro put on a two asset min basket
189	50	case Option::Put:
190	50	results_.value = strike * riskFreeDiscount -
191	50	euroTwoAssetMinBasketCall(forward1, forward2, 0.0,
192	50	riskFreeDiscount,
193	50	variance1, variance2, rho_) +
194	50	euroTwoAssetMinBasketCall(forward1, forward2, strike,
195	50	riskFreeDiscount,
196	50	variance1, variance2, rho_);
197	50	break;
198	0	default:
199	0	QL_FAIL("unknown option type");
200	78	}
201	78	} else {
202		QL_FAIL("unknown type");
203	0	}
204	385	}
205
206		}