[go: up one dir, main page]
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2018 Klaus Spanderen
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy of the license along with this program; if not, please email
12 <quantlib-dev@lists.sf.net>. The license is also available online at
13 <http://quantlib.org/license.shtml>.
14
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20#include "fdcev.hpp"
21#include "utilities.hpp"
22#include <ql/math/randomnumbers/rngtraits.hpp>
23#include <ql/math/integrals/gausslobattointegral.hpp>
24#include <ql/math/statistics/generalstatistics.hpp>
25#include <ql/pricingengines/vanilla/analyticcevengine.hpp>
26#include <ql/pricingengines/vanilla/fdcevvanillaengine.hpp>
27#include <ql/methods/finitedifferences/utilities/cevrndcalculator.hpp>
28#include <ql/shared_ptr.hpp>
29
30using namespace QuantLib;
31using boost::unit_test_framework::test_suite;
32
33
34namespace {
35 class ExpectationFct {
36 public:
37 ExpectationFct(const CEVRNDCalculator& calculator, Time t)
38 : t_(t), calculator_(calculator) { }
39
40 Real operator()(Real f) const { return f*calculator_.pdf(f, t: t_); }
41
42 private:
43 const Time t_;
44 const CEVRNDCalculator& calculator_;
45 };
46}
47
48void FdCevTest::testLocalMartingale() {
49 BOOST_TEST_MESSAGE(
50 "Testing local martingale property of CEV process with PDF...");
51
52 const Time t = 1.0;
53
54 const Real f0 = 2.1;
55 const Real alpha = 1.75;
56 const Real betas[] = {-2.4, 0.23, 0.9, 1.1, 1.5};
57
58 for (Real beta : betas) {
59 const CEVRNDCalculator rndCalculator(f0, alpha, beta);
60
61 const Real eps = 1e-10;
62 const Real tol = 100*eps;
63
64 const Real upperBound = 10*rndCalculator.invcdf(q: 1-eps, t);
65
66 const Real expectationValue = GaussLobattoIntegral(10000, eps)(
67 ExpectationFct(rndCalculator, t), QL_EPSILON, upperBound);
68
69 const Real diff = expectationValue-f0;
70
71
72 if (beta < 1.0 && std::fabs(x: diff) > tol) {
73 BOOST_ERROR("CEV process should be a martingale for beta < 1.0"
74 << "\n expected: " << f0
75 << std::scientific
76 << "\n difference " << diff
77 << "\n tolerance: " << tol);
78 }
79
80 if (beta > 1.0 && diff > -tol) {
81 BOOST_ERROR("CEV process should only be a local martingale "
82 "for beta > 1.0. Expectation is E[F_t|F_0] < F_0"
83 << "\n E[F_t|F_0]: " << expectationValue
84 << "\n F_0: " << f0);
85 }
86
87 // check local martingale property with Monte-Carlo simulation
88 const Size nSims = 5000;
89
90 const Size nSteps = 2000;
91 const Real dt = t / nSteps;
92 const Real sqrtDt = std::sqrt(x: dt);
93
94 GeneralStatistics stat;
95 const PseudoRandom::rng_type mt(MersenneTwisterUniformRng(42));
96
97 if (beta > 1.2) {
98 for (Size i=0; i < nSims; ++i) {
99 Real f = f0;
100 for (Size j=0; j < nSteps; ++j) {
101 f += alpha * std::pow(x: f, y: beta) * mt.next().value * sqrtDt;
102 f = std::max(a: 0.0, b: f);
103
104 if (f == 0.0) break; // absorbing boundary
105 }
106 stat.add(value: f - f0);
107 }
108
109 const Real calculated = stat.mean();
110 const Real error = stat.errorEstimate();
111
112 if (std::fabs(x: calculated - diff) > 2.35*error) {
113 BOOST_ERROR(
114 "failed to calculate local martingale property "
115 "by Monte-Carlo Simulation for beta > 1.0. "
116 << "\n E[F_t|F_0] : " << expectationValue
117 << "\n E_MC[F_t|F_0]: " << calculated + f0
118 << "\n error_MC : " << error
119 << "\n difference : " << std::fabs(calculated - diff)
120 << "\n tolerance : " << 2.35*error);
121 }
122 }
123 }
124}
125
126void FdCevTest::testFdmCevOp() {
127 BOOST_TEST_MESSAGE(
128 "Testing FDM constant elasticity of variance (CEV) operator...");
129
130 const Date today = Date(22, February, 2018);
131 const DayCounter dc = Actual365Fixed();
132 Settings::instance().evaluationDate() = today;
133
134 const Date maturityDate = today + Period(12, Months);
135 const Real strike = 2.3;
136
137 const Option::Type optionTypes[] = { Option::Call, Option::Put};
138
139 const ext::shared_ptr<Exercise> exercise =
140 ext::make_shared<EuropeanExercise>(args: maturityDate);
141
142 for (auto optionType : optionTypes) {
143 const ext::shared_ptr<PlainVanillaPayoff> payoff =
144 ext::make_shared<PlainVanillaPayoff>(args&: optionType, args: strike);
145
146 const ext::shared_ptr<YieldTermStructure> rTS =
147 flatRate(today, forward: 0.15, dc);
148
149 const Real f0 = 2.1;
150 const Real alpha = 0.75;
151
152 const Real betas[] = { -2.0, -0.5, 0.45, 0.6, 0.9, 1.45 };
153 for (Real beta : betas) {
154
155 VanillaOption option(payoff, exercise);
156 option.setPricingEngine(ext::make_shared<AnalyticCEVEngine>(
157 args: f0, args: alpha, args&: beta, args: Handle<YieldTermStructure>(rTS)));
158
159 const Real analyticNPV = option.NPV();
160
161 const Real eps = 1e-3;
162
163 option.setPricingEngine(ext::make_shared<AnalyticCEVEngine>(
164 args: f0*(1+eps), args: alpha, args&: beta, args: Handle<YieldTermStructure>(rTS)));
165 const Real analyticUpNPV = option.NPV();
166
167 option.setPricingEngine(ext::make_shared<AnalyticCEVEngine>(
168 args: f0*(1-eps), args: alpha, args&: beta, args: Handle<YieldTermStructure>(rTS)));
169 const Real analyticDownNPV = option.NPV();
170
171 const Real analyticDelta = (analyticUpNPV - analyticDownNPV)
172 /(2*eps*f0);
173
174 option.setPricingEngine(ext::make_shared<FdCEVVanillaEngine>(
175 args: f0, args: alpha, args&: beta, args: Handle<YieldTermStructure>(rTS),
176 args: 100, args: 1000, args: 1, args: 1.0, args: 1e-6));
177
178 const Real calculatedNPV = option.NPV();
179 const Real calculatedDelta = option.delta();
180
181 const Real tol = 0.01;
182 if (std::fabs(x: calculatedNPV - analyticNPV) > tol
183 || std::fabs(x: calculatedDelta - analyticDelta) > tol) {
184 BOOST_ERROR(
185 "failed to calculate vanilla option prices/delta "
186 << "\n beta : " << beta
187 << "\n option type : "
188 << ((payoff->optionType() == Option::Call) ? "Call" : "Put")
189 << "\n analytic npv : " << analyticNPV
190 << "\n pde npv : " << calculatedNPV
191 << "\n npv difference : "
192 << std::fabs(calculatedNPV - analyticNPV)
193 << "\n tolerance : " << tol
194 << "\n analytic delta : " << analyticDelta
195 << "\n pde delta : " << calculatedDelta
196 << "\n delta difference: "
197 << std::fabs(calculatedDelta - analyticDelta)
198 << "\n tolerance : " << tol);
199 }
200 }
201 }
202}
203
204
205test_suite* FdCevTest::suite(SpeedLevel speed) {
206 auto* suite = BOOST_TEST_SUITE("Finite Difference CEV tests");
207
208
209 suite->add(QUANTLIB_TEST_CASE(&FdCevTest::testLocalMartingale));
210 suite->add(QUANTLIB_TEST_CASE(&FdCevTest::testFdmCevOp));
211
212 return suite;
213}
214

source code of quantlib/test-suite/fdcev.cpp