Logo
Fully annotated reference manual - version 1.8.12
Loading...
Searching...
No Matches
crcirpp.cpp
Go to the documentation of this file.
1/*
2 Copyright (C) 2020 Quaternion Risk Management Ltd
3 All rights reserved.
4
5 This file is part of ORE, a free-software/open-source library
6 for transparent pricing and risk analysis - http://opensourcerisk.org
7
8 ORE is free software: you can redistribute it and/or modify it
9 under the terms of the Modified BSD License. You should have received a
10 copy of the license along with this program.
11 The license is also available online at <http://opensourcerisk.org>
12
13 This program is distributed on the basis that it will form a useful
14 contribution to risk analytics and model standardisation, but WITHOUT
15 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the license for more details.
17*/
18
19#include "utilities.hpp"
20#include "toplevelfixture.hpp"
21#include <boost/test/unit_test.hpp>
22#include <ql/currencies/europe.hpp>
23#include <ql/experimental/callablebonds/callablebond.hpp>
24#include <ql/indexes/ibor/euribor.hpp>
25#include <ql/instruments/callabilityschedule.hpp>
26#include <ql/instruments/creditdefaultswap.hpp>
27#include <ql/math/optimization/levenbergmarquardt.hpp>
28#include <ql/math/randomnumbers/rngtraits.hpp>
29#include <ql/math/statistics/incrementalstatistics.hpp>
30#include <ql/methods/montecarlo/multipathgenerator.hpp>
31#include <ql/methods/montecarlo/pathgenerator.hpp>
32#include <ql/quotes/simplequote.hpp>
33#include <ql/termstructures/credit/flathazardrate.hpp>
34#include <ql/termstructures/yield/flatforward.hpp>
35#include <ql/time/calendars/target.hpp>
36#include <ql/time/daycounters/actual360.hpp>
37#include <ql/time/daycounters/actualactual.hpp>
38#include <ql/time/daycounters/thirty360.hpp>
39#include <ql/pricingengines/credit/midpointcdsengine.hpp>
40
41#include <qle/instruments/cdsoption.hpp>
42#include <qle/methods/multipathgeneratorbase.hpp>
43#include <qle/models/crossassetmodel.hpp>
44#include <qle/models/cdsoptionhelper.hpp>
45#include <qle/models/crcirpp.hpp>
46#include <qle/models/cirppconstantfellerparametrization.hpp>
47#include <qle/models/crossassetmodel.hpp>
48#include <qle/processes/crcirppstateprocess.hpp>
49
50#include <boost/make_shared.hpp>
51
52#include <fstream>
53#include <iostream>
54
55// fix for boost 1.64, see https://lists.boost.org/Archives/boost/2016/11/231756.php
56#if BOOST_VERSION >= 106400
57#include <boost/serialization/array_wrapper.hpp>
58#endif
59#include <boost/accumulators/accumulators.hpp>
60#include <boost/accumulators/statistics/covariance.hpp>
61#include <boost/accumulators/statistics/density.hpp>
62#include <boost/accumulators/statistics/error_of_mean.hpp>
63#include <boost/accumulators/statistics/mean.hpp>
64#include <boost/accumulators/statistics/stats.hpp>
65#include <boost/accumulators/statistics/variates/covariate.hpp>
66#include <boost/make_shared.hpp>
67
68using namespace QuantLib;
69using namespace QuantExt;
70
71using namespace boost::accumulators;
72
73namespace {
74struct CreditModelTestData_flat: public qle::test::TopLevelFixture {
75 CreditModelTestData_flat()
76 : referenceDate(29, July, 2017), dts(QuantLib::ext::make_shared<FlatHazardRate>(referenceDate, 0.04, ActualActual(ActualActual::ISDA))),
77 yts(QuantLib::ext::make_shared<FlatForward>(referenceDate, 0.02, ActualActual(ActualActual::ISDA))) {
78
79 Settings::instance().evaluationDate() = referenceDate;
80
81 kappa = 0.206;
82 theta = 0.04;
83 sigma = sqrt(2 * kappa * theta) - 1E-10;
84 y0 = theta;
85 shifted = true;
86
87 recoveryRate = 0.4;
88 // QL_REQUIRE(2 * eurKappa * eurTheta / eurSigma / eurSigma > 1 , "Feller Condition is not satisfied!");
89 cirParametrization = QuantLib::ext::make_shared<CrCirppConstantWithFellerParametrization>(
90 EURCurrency(), dts, kappa, theta, sigma, y0, shifted);
91 QL_REQUIRE(cirParametrization != NULL, "CrCirppConstantWithFellerParametrization has null pointer!");
92
93 model = QuantLib::ext::make_shared<CrCirpp>(cirParametrization);
94 BOOST_TEST_MESSAGE("CIR++ parameters: ");
95 BOOST_TEST_MESSAGE("Kappa: \t" << model->parametrization()->kappa(0));
96 BOOST_TEST_MESSAGE("Theta: \t" << model->parametrization()->theta(0));
97 BOOST_TEST_MESSAGE("Sigma: \t" << model->parametrization()->sigma(0));
98 BOOST_TEST_MESSAGE("y0: \t" << model->parametrization()->y0(0));
99 BOOST_TEST_MESSAGE("Feller condition is (>1 ok) "
100 << 2.0 * model->parametrization()->kappa(0) * model->parametrization()->theta(0) /
101 model->parametrization()->sigma(0) / model->parametrization()->sigma(0));
102 }
103
104 SavedSettings backup;
105 Date referenceDate;
106 Handle<DefaultProbabilityTermStructure> dts;
107 Handle<YieldTermStructure> yts;
108 Real kappa, theta, sigma, y0;
109 bool shifted;
110 Real recoveryRate;
111 QuantLib::ext::shared_ptr<CrCirppConstantWithFellerParametrization> cirParametrization;
112 QuantLib::ext::shared_ptr<CrCirpp> model;
113}; // IrTSModelTestData
114} // namespace
115
116BOOST_FIXTURE_TEST_SUITE(QuantExtTestSuite, CreditModelTestData_flat)
117
118BOOST_AUTO_TEST_SUITE(CrCirppModelTest)
119
120BOOST_AUTO_TEST_CASE(testMartingaleProperty) {
121
122 BOOST_TEST_MESSAGE("Testing martingale property in credit-CIR++ model for Brigo-Alfonsi discretizations...");
123
124 QuantLib::ext::shared_ptr<StochasticProcess> process = model->stateProcess();
125 QL_REQUIRE(process != NULL, "process has null pointer!");
126
127 // CIRplusplusStateProcess::Discretization discType = CIRplusplusStateProcess::Discretization::Reflection;
128 // CIRplusplusStateProcess::Discretization discType = CIRplusplusStateProcess::Discretization::PartialTruncation;
129 // CIRplusplusStateProcess::Discretization discType = CIRplusplusStateProcess::Discretization::FullTruncation;
130 // QuantLib::ext::shared_ptr<StochasticProcess> process = QuantLib::ext::make_shared<CIRplusplusStateProcess>(model.get(),
131 // discType); BOOST_TEST_MESSAGE("Simulation type of negative variance process "<<discType);
132 Size n = 10000; // number of paths
133 Size seed = 42; // rng seed
134 // Time T = 25.0; // maturity of payoff
135 // Time T2 = 40.0; // zerobond maturity
136 Time T = 10.0; // maturity of payoff
137 Time T2 = 20.0; // zerobond maturity
138 // QL_REQUIRE(T2 == model->maxSimulationHorizon(), "Forward measure horizon must be equalt to the maturity of the
139 // zero bond");
140 Size steps = static_cast<Size>(T * 52); // number of steps taken (euler)
141
142 TimeGrid grid(T, steps);
143 // LowDiscrepancy::rsg_type sg = LowDiscrepancy::make_sequence_generator( 9 * steps, seed);
144 // MultiPathGenerator<LowDiscrepancy::rsg_type> pg(process, grid, sg, false);
145 MultiPathGeneratorMersenneTwister pg(process, grid, seed, true);
146 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean>>> meanTest_y;
147 accumulator_set<double, stats<tag::variance, tag::error_of<tag::mean>>> varTest_y;
148
149 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean>>> sp;
150 accumulator_set<double, stats<tag::mean, tag::error_of<tag::mean>>> numeraire;
151
152 accumulator_set<double, stats<tag::density>> histXAcc(tag::density::cache_size = 10000,
153 tag::density::num_bins = 50);
154
155 // typedef boost::iterator_range<std::vector<std::pair<double, double>>::iterator> histogram_type;
156
157 // //create an accumulator
158 // acc histXAcc( tag::density::num_bins = 20, tag::density::cache_size = 10);
159
160 std::ofstream of;
161 // of.open("y_paths.txt");
162 for (Size j = 0; j < n; ++j) {
163 Sample<MultiPath> path = pg.next();
164 Size l = path.value[0].length() - 1;
165 Real y = path.value[0][l];
166 Real num = path.value[1][l];
167 sp(model->survivalProbability(T, T2, y) * num);
168 numeraire(num);
169 meanTest_y(y);
170 varTest_y(y);
171 histXAcc(y);
172 }
173 // histogram_type histX = density(histXAcc);
174 // double total = 0.0;
175 // std::ofstream ofHistX;
176 // ofHistX.open("y_hist.txt");
177
178 // Real binSize = histX[1].first - histX[0].first;
179 // for (Size i = 0; i < histX.size(); i++) {
180 // ofHistX << histX[i].first << "," << histX[i].second / binSize << std::endl;
181 // total += histX[i].second;
182 // }
183 // BOOST_TEST_MESSAGE("Total cdf(x): " << total); // should be 1 (and it is)
184
185 // total = 0.0;
186 // std::ofstream ofPDFX;
187 // ofPDFX.open("y_pdf.txt");
188 // // for (Real x = - 5.0; x <= 7.0; x += 0.1)
189 // for (Real x = 0.00018; x <= 0.01; x += 0.001) {
190 // ofPDFX << x << "," << model->density(x, T) << std::endl;
191 // total += model->density(x, T);
192 // }
193 // BOOST_TEST_MESSAGE("Total cdf(x): " << total); // should be 1 (and it is)
194
195 // of.close();
196 // ofPDFX.close();
197 // ofHistX.close();
198
199 // Real kappa = model->parametrization()->kappa(0);
200 // Real theta = model->parametrization()->theta(0);
201 // Real sigma = model->parametrization()->sigma(0);
202 // Real y0 = model->parametrization()->y0(0);
203 // Real sigma2 = sigma * sigma;
204
205 BOOST_TEST_MESSAGE("\nBrigo-Alfonsi:");
206 // BOOST_TEST_MESSAGE("Mean = " << mean(meanTest_y) << " +- " << error_of<tag::mean>(meanTest_y) << " analytic "
207 // <<y0 * exp(-kappa*T) + theta * (1- exp(-kappa*T))); BOOST_TEST_MESSAGE("Variance = " << variance(varTest_y) << "
208 // +- " << error_of<tag::mean>(varTest_y)<<" analytic " <<y0*sigma2/theta*(exp(-kappa*T) - exp(-2*kappa*T)) +
209 // theta*sigma2/2/kappa*(1- exp(-2*kappa *T))*(1- exp(-2*kappa *T)) );
210 BOOST_TEST_MESSAGE("SP = " << mean(sp) << " +- " << error_of<tag::mean>(sp) << " vs analytical "
211 << dts->survivalProbability(T2));
212 BOOST_TEST_MESSAGE("Num = " << mean(numeraire) << " +- " << error_of<tag::mean>(numeraire) << " vs analytical "
213 << dts->survivalProbability(T));
214
215 Real tol2 = 12.0E-4;
216 Real expectedSP = dts->survivalProbability(T);
217 Real expectedCondSP = dts->survivalProbability(T2);
218 if (std::abs(mean(numeraire) - expectedSP) > tol2)
219 BOOST_FAIL("Martingale test failed for SP(t) (Brigo-Alfonsi discr.), expected "
220 << expectedSP << ", got " << mean(numeraire) << ", tolerance " << tol2);
221 if (std::abs(mean(sp) - expectedCondSP) > tol2)
222 BOOST_FAIL("Martingale test failed for SP(t,T) (Brigo-Alfonsi discr.), expected "
223 << expectedCondSP << ", got " << mean(sp) << ", tolerance " << tol2);
224
225} // testIrTSMartingaleProperty
226
227
228BOOST_AUTO_TEST_SUITE_END()
229
230BOOST_AUTO_TEST_SUITE_END()
cdsoption.hpp
CDS option, removed requirements (rec must knock out, no upfront amount), that should be taken care o...
cdsoptionhelper.hpp
cds option calibration helper
cirppconstantfellerparametrization.hpp
constant CIR ++ parametrization
QuantExt::MultiPathGeneratorMersenneTwister
Instantiation of MultiPathGenerator with standard PseudoRandom traits.
Definition: multipathgeneratorbase.hpp:59
QuantExt::MultiPathGeneratorMersenneTwister::next
const Sample< MultiPath > & next() const override
Definition: multipathgeneratorbase.cpp:44
crcirpp.hpp
CIR++ credit model class.
crcirppstateprocess.hpp
CIR++ model state process.
crossassetmodel.hpp
cross asset model
multipathgeneratorbase.hpp
base class for multi path generators
QuantExt::sqrt
RandomVariable sqrt(RandomVariable x)
Definition: randomvariable.cpp:771
steps
std::vector< Size > steps
Definition: commodityschwartzmodel.cpp:109
BOOST_AUTO_TEST_CASE
BOOST_AUTO_TEST_CASE(testMartingaleProperty)
Definition: crcirpp.cpp:120
toplevelfixture.hpp
Fixture that can be used at top level.
utilities.hpp
helper macros and methods for tests
Generated by Doxygen 1.9.5