[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) 2010, 2011, 2012 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 "vpp.hpp"
21#include "utilities.hpp"
22#include <ql/experimental/finitedifferences/dynprogvppintrinsicvalueengine.hpp>
23#include <ql/experimental/finitedifferences/fdklugeextouspreadengine.hpp>
24#include <ql/experimental/finitedifferences/fdmklugeextouop.hpp>
25#include <ql/experimental/finitedifferences/fdmspreadpayoffinnervalue.hpp>
26#include <ql/experimental/finitedifferences/fdmvppstepconditionfactory.hpp>
27#include <ql/experimental/finitedifferences/fdsimpleextoustorageengine.hpp>
28#include <ql/experimental/finitedifferences/fdsimpleklugeextouvppengine.hpp>
29#include <ql/experimental/finitedifferences/vanillavppoption.hpp>
30#include <ql/experimental/processes/extendedornsteinuhlenbeckprocess.hpp>
31#include <ql/experimental/processes/extouwithjumpsprocess.hpp>
32#include <ql/experimental/processes/gemanroncoroniprocess.hpp>
33#include <ql/experimental/processes/klugeextouprocess.hpp>
34#include <ql/instruments/basketoption.hpp>
35#include <ql/instruments/vanillastorageoption.hpp>
36#include <ql/instruments/vanillaswingoption.hpp>
37#include <ql/math/generallinearleastsquares.hpp>
38#include <ql/math/randomnumbers/rngtraits.hpp>
39#include <ql/math/statistics/generalstatistics.hpp>
40#include <ql/math/functional.hpp>
41#include <ql/methods/finitedifferences/meshers/exponentialjump1dmesher.hpp>
42#include <ql/methods/finitedifferences/meshers/fdmmeshercomposite.hpp>
43#include <ql/methods/finitedifferences/meshers/fdmsimpleprocess1dmesher.hpp>
44#include <ql/methods/finitedifferences/meshers/uniform1dmesher.hpp>
45#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
46#include <ql/methods/finitedifferences/utilities/fdminnervaluecalculator.hpp>
47#include <ql/methods/montecarlo/lsmbasissystem.hpp>
48#include <ql/methods/montecarlo/multipathgenerator.hpp>
49#include <ql/processes/ornsteinuhlenbeckprocess.hpp>
50#include <ql/processes/stochasticprocessarray.hpp>
51#include <ql/quotes/simplequote.hpp>
52#include <ql/termstructures/yield/zerocurve.hpp>
53#include <ql/time/daycounters/actualactual.hpp>
54#include <deque>
55#include <utility>
56
57using namespace QuantLib;
58using namespace boost::unit_test_framework;
59
60namespace vpp_test {
61
62 ext::function<Real(Real)> constant_b(Real b) {
63 return [=](Real x){ return b; };
64 }
65
66 ext::shared_ptr<ExtOUWithJumpsProcess> createKlugeProcess() {
67 Array x0(2);
68 x0[0] = 3.0; x0[1] = 0.0;
69
70 const Real beta = 5.0;
71 const Real eta = 2.0;
72 const Real jumpIntensity = 1.0;
73 const Real speed = 1.0;
74 const Real volatility = 2.0;
75
76 ext::shared_ptr<ExtendedOrnsteinUhlenbeckProcess> ouProcess(
77 new ExtendedOrnsteinUhlenbeckProcess(speed, volatility, x0[0],
78 constant_b(b: x0[0])));
79 return ext::make_shared<ExtOUWithJumpsProcess>(
80 args&: ouProcess, args&: x0[1], args: beta,
81 args: jumpIntensity, args: eta);
82 }
83
84 class linear {
85 Real alpha, beta;
86 public:
87 linear(Real alpha, Real beta) : alpha(alpha), beta(beta) {}
88 Real operator()(Real x) const {
89 return alpha + beta*x;
90 }
91 };
92
93}
94
95
96void VPPTest::testGemanRoncoroniProcess() {
97
98 BOOST_TEST_MESSAGE("Testing Geman-Roncoroni process...");
99
100 /* Example induced by H. Geman, A. Roncoroni,
101 "Understanding the Fine Structure of Electricity Prices",
102 http://papers.ssrn.com/sol3/papers.cfm?abstract_id=638322
103 Results are verified against the authors MatLab-Code.
104 http://semeq.unipmn.it/files/Ch19_spark_spread.zip
105 */
106
107 using namespace vpp_test;
108
109 const Date today = Date(18, December, 2011);
110 Settings::instance().evaluationDate() = today;
111 const DayCounter dc = ActualActual(ActualActual::ISDA);
112
113 ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, forward: 0.03, dc);
114
115 const Real x0 = 3.3;
116 const Real beta = 0.05;
117 const Real alpha = 3.1;
118 const Real gamma = -0.09;
119 const Real delta = 0.07;
120 const Real eps = -0.40;
121 const Real zeta = -0.40;
122 const Real d = 1.6;
123 const Real k = 1.0;
124 const Real tau = 0.5;
125 const Real sig2 = 10.0;
126 const Real a =-7.0;
127 const Real b =-0.3;
128 const Real theta1 = 35.0;
129 const Real theta2 = 9.0;
130 const Real theta3 = 0.10;
131 const Real psi = 1.9;
132
133 ext::shared_ptr<GemanRoncoroniProcess> grProcess(
134 new GemanRoncoroniProcess(x0, alpha, beta, gamma, delta,
135 eps, zeta, d, k, tau, sig2, a, b,
136 theta1, theta2, theta3, psi));
137
138
139 const Real speed = 5.0;
140 const Volatility vol = std::sqrt(x: 1.4);
141 const Real betaG = 0.08;
142 const Real alphaG = 1.0;
143 const Real x0G = 1.1;
144
145 ext::function<Real (Real)> f = vpp_test::linear(alphaG, betaG);
146
147 ext::shared_ptr<StochasticProcess1D> eouProcess(
148 new ExtendedOrnsteinUhlenbeckProcess(speed, vol, x0G, f,
149 ExtendedOrnsteinUhlenbeckProcess::Trapezodial));
150
151 std::vector<ext::shared_ptr<StochasticProcess1D> > processes = {grProcess, eouProcess};
152
153 Matrix correlation(2, 2, 1.0);
154 correlation[0][1] = correlation[1][0] = 0.25;
155
156 ext::shared_ptr<StochasticProcess> pArray(
157 new StochasticProcessArray(processes, correlation));
158
159 const Time T = 10.0;
160 const Size stepsPerYear = 250;
161 const Size steps = Size(T*Real(stepsPerYear));
162
163 TimeGrid grid(T, steps);
164
165 typedef PseudoRandom::rsg_type rsg_type;
166 typedef MultiPathGenerator<rsg_type>::sample_type sample_type;
167 rsg_type rsg = PseudoRandom::make_sequence_generator(
168 dimension: pArray->size()*(grid.size()-1), seed: BigNatural(421));
169
170 GeneralStatistics npv, onTime;
171 MultiPathGenerator<rsg_type> generator(pArray, grid, rsg, false);
172
173 const Real heatRate = 8.0;
174 const Size nrTrails = 250;
175
176 for (Size n=0; n < nrTrails; ++n) {
177 Real plantValue = 0.0;
178 sample_type path = generator.next();
179
180 for (Size i=1; i <= steps; ++i) {
181 const Time t = Real(i)/stepsPerYear;
182 const DiscountFactor df = rTS->discount(t);
183
184 const Real fuelPrice = std::exp(x: path.value[1][i]);
185 const Real electricityPrice = std::exp(x: path.value[0][i]);
186
187 const Real sparkSpread = electricityPrice - heatRate*fuelPrice;
188 plantValue += std::max(a: 0.0, b: sparkSpread)*df;
189 onTime.add(value: (sparkSpread > 0.0) ? 1.0 : 0.0);
190 }
191
192 npv.add(value: plantValue);
193 }
194
195 const Real expectedNPV = 12500;
196 const Real calculatedNPV = npv.mean();
197 const Real errorEstimateNPV = npv.errorEstimate();
198
199 if (std::fabs(x: calculatedNPV - expectedNPV) > 3.0*errorEstimateNPV) {
200 BOOST_ERROR("Failed to reproduce cached price with MC engine"
201 << "\n calculated: " << calculatedNPV
202 << "\n expected: " << expectedNPV
203 << " +/- " << errorEstimateNPV);
204 }
205
206 const Real expectedOnTime = 0.43;
207 const Real calculatedOnTime = onTime.mean();
208 const Real errorEstimateOnTime
209 = std::sqrt(x: calculatedOnTime*(1-calculatedOnTime))/nrTrails;
210
211 if (std::fabs(x: calculatedOnTime - expectedOnTime)>3.0*errorEstimateOnTime) {
212 BOOST_ERROR("Failed to reproduce cached price with MC engine"
213 << "\n calculated: " << calculatedNPV
214 << "\n expected: " << expectedNPV
215 << " +/- " << errorEstimateNPV);
216 }
217}
218
219void VPPTest::testSimpleExtOUStorageEngine() {
220
221 BOOST_TEST_MESSAGE("Testing simple-storage option based on ext. OU model...");
222
223 using namespace vpp_test;
224
225 Date settlementDate = Date(18, December, 2011);
226 Settings::instance().evaluationDate() = settlementDate;
227 DayCounter dayCounter = ActualActual(ActualActual::ISDA);
228 Date maturityDate = settlementDate + Period(12, Months);
229
230 std::vector<Date> exerciseDates(1, settlementDate+Period(1, Days));
231 while (exerciseDates.back() < maturityDate) {
232 exerciseDates.push_back(x: exerciseDates.back()+Period(1, Days));
233 }
234 ext::shared_ptr<BermudanExercise> bermudanExercise(
235 new BermudanExercise(exerciseDates));
236
237 const Real x0 = 3.0;
238 const Real speed = 1.0;
239 const Real volatility = 0.5;
240 const Rate irRate = 0.1;
241
242 ext::shared_ptr<ExtendedOrnsteinUhlenbeckProcess> ouProcess(
243 new ExtendedOrnsteinUhlenbeckProcess(speed, volatility, x0,
244 constant_b(b: x0)));
245
246 ext::shared_ptr<YieldTermStructure> rTS(
247 flatRate(today: settlementDate, forward: irRate, dc: dayCounter));
248
249 ext::shared_ptr<PricingEngine> storageEngine(
250 new FdSimpleExtOUStorageEngine(ouProcess, rTS, 1, 25));
251
252 VanillaStorageOption storageOption(bermudanExercise, 50, 0, 1);
253
254 storageOption.setPricingEngine(storageEngine);
255
256 const Real expected = 69.5755;
257 const Real calculated = storageOption.NPV();
258
259 if (std::fabs(x: expected - calculated) > 5e-2) {
260 BOOST_ERROR("Failed to reproduce cached values" <<
261 "\n calculated: " << calculated <<
262 "\n expected: " << expected);
263 }
264}
265
266
267void VPPTest::testKlugeExtOUSpreadOption() {
268
269 BOOST_TEST_MESSAGE("Testing simple Kluge ext-Ornstein-Uhlenbeck spread option...");
270
271 using namespace vpp_test;
272
273 Date settlementDate = Date(18, December, 2011);
274 Settings::instance().evaluationDate() = settlementDate;
275
276 DayCounter dayCounter = ActualActual(ActualActual::ISDA);
277 Date maturityDate = settlementDate + Period(1, Years);
278 Time maturity = dayCounter.yearFraction(d1: settlementDate, d2: maturityDate);
279
280 const Real speed = 1.0;
281 const Volatility vol = std::sqrt(x: 1.4);
282 const Real betaG = 0.0;
283 const Real alphaG = 3.0;
284 const Real x0G = 3.0;
285
286 const Rate irRate = 0.0;
287 const Real heatRate = 2.0;
288 const Real rho = 0.5;
289
290 ext::shared_ptr<ExtOUWithJumpsProcess>
291 klugeProcess = createKlugeProcess();
292 ext::function<Real (Real)> f = vpp_test::linear(alphaG, betaG);
293
294 ext::shared_ptr<ExtendedOrnsteinUhlenbeckProcess> extOUProcess(
295 new ExtendedOrnsteinUhlenbeckProcess(speed, vol, x0G, f,
296 ExtendedOrnsteinUhlenbeckProcess::Trapezodial));
297
298 ext::shared_ptr<YieldTermStructure> rTS(
299 flatRate(today: settlementDate, forward: irRate, dc: dayCounter));
300
301 ext::shared_ptr<KlugeExtOUProcess> klugeOUProcess(
302 new KlugeExtOUProcess(rho, klugeProcess, extOUProcess));
303
304
305 ext::shared_ptr<Payoff> payoff(new PlainVanillaPayoff(Option::Call, 0.0));
306
307 Array spreadFactors(2);
308 spreadFactors[0] = 1.0; spreadFactors[1] = -heatRate;
309 ext::shared_ptr<BasketPayoff> basketPayoff(
310 new AverageBasketPayoff(payoff, spreadFactors));
311
312 ext::shared_ptr<Exercise> exercise(new EuropeanExercise(maturityDate));
313
314 BasketOption option(basketPayoff, exercise);
315 option.setPricingEngine(ext::shared_ptr<PricingEngine>(
316 new FdKlugeExtOUSpreadEngine(klugeOUProcess, rTS,
317 5, 200, 50, 20)));
318
319 TimeGrid grid(maturity, 50);
320 typedef PseudoRandom::rsg_type rsg_type;
321 typedef MultiPathGenerator<rsg_type>::sample_type sample_type;
322
323 rsg_type rsg = PseudoRandom::make_sequence_generator(
324 dimension: klugeOUProcess->factors() * (grid.size() - 1), seed: 1234UL);
325
326 MultiPathGenerator<rsg_type> generator(klugeOUProcess, grid, rsg, false);
327
328
329 GeneralStatistics npv;
330 const Size nTrails = 20000;
331 for (Size i=0; i < nTrails; ++i) {
332 const sample_type& path = generator.next();
333
334 Array p(2);
335 p[0] = path.value[0].back() + path.value[1].back();
336 p[1] = path.value[2].back();
337 npv.add(value: (*basketPayoff)(Exp(v: p)));
338 }
339
340 const Real calculated = option.NPV();
341 const Real expectedMC = npv.mean();
342 const Real mcError = npv.errorEstimate();
343 if (std::fabs(x: expectedMC - calculated) > 3*mcError) {
344 BOOST_ERROR("Failed to reproduce referenc values"
345 << "\n calculated: " << calculated
346 << "\n expected(MC): " << expectedMC
347 << "\n mc error : " << mcError);
348
349 }
350}
351
352namespace vpp_test {
353 // for a "real" gas and power forward curve
354 // please see. e.g. http://www.kyos.com/?content=64
355 const std::vector<Real> fuelPrices = {
356 20.74,21.65,20.78,21.58,21.43,20.82,22.02,21.52,
357 21.02,21.46,21.75,20.69,22.16,20.38,20.82,20.68,
358 20.57,21.92,22.04,20.45,20.75,21.92,20.53,20.67,
359 20.88,21.02,20.82,21.67,21.82,22.12,20.45,20.74,
360 22.39,20.95,21.71,20.70,20.94,21.59,22.33,21.13,
361 21.50,21.42,20.56,21.23,21.37,21.90,20.62,21.17,
362 21.86,22.04,22.05,21.00,20.70,21.12,21.26,22.40,
363 21.31,22.24,21.96,21.02,21.71,20.48,21.36,21.75,
364 21.90,20.44,21.26,22.29,20.34,21.79,21.66,21.50,
365 20.76,20.27,20.84,20.24,21.97,20.52,20.98,21.40,
366 20.39,20.71,20.78,20.30,21.56,21.72,20.27,21.57,
367 21.82,20.57,21.33,20.51,22.32,21.99,20.57,22.11,
368 21.56,22.24,20.62,21.70,21.11,21.19,21.79,20.46,
369 22.21,20.82,20.52,22.29,20.71,21.45,22.40,20.63,
370 20.95,21.97,22.20,20.67,21.01,22.25,20.76,21.33,
371 20.49,20.33,21.94,20.64,20.99,21.09,20.97,22.17,
372 20.72,22.06,20.86,21.40,21.75,20.78,21.79,20.47,
373 21.19,21.60,20.75,21.36,21.61,20.37,21.67,20.28,
374 22.33,21.37,21.33,20.87,21.25,22.01,22.08,20.81,
375 20.70,21.84,21.82,21.68,21.24,22.36,20.83,20.64,
376 21.03,20.57,22.34,20.96,21.54,21.26,21.43,22.39};
377
378 const std::vector<Real> powerPrices = {
379 40.40,36.71,31.87,25.81,31.61,35.00,46.22,60.68,
380 42.45,38.01,33.84,29.79,31.84,38.53,49.23,59.92,
381 43.85,37.47,34.89,29.99,30.85,29.19,29.25,38.67,
382 36.90,25.93,22.12,20.19,17.19,19.29,13.51,18.14,
383 33.76,30.48,25.63,18.01,23.86,32.41,48.56,64.69,
384 38.42,39.31,32.73,29.97,31.41,35.02,46.85,58.12,
385 39.14,35.42,32.61,28.76,29.41,35.83,46.73,61.41,
386 61.01,59.43,60.43,66.29,62.79,62.66,57.66,51.63,
387 62.18,60.53,61.94,64.86,59.57,58.15,53.74,48.36,
388 45.64,51.21,51.54,50.79,54.50,49.92,41.58,39.81,
389 28.86,37.42,39.78,42.36,45.67,36.84,33.91,28.75,
390 62.97,63.84,62.91,68.77,64.33,61.95,59.12,54.89,
391 63.62,60.90,66.57,69.51,64.71,59.89,57.28,57.10,
392 65.09,63.82,67.52,70.51,65.59,59.36,58.22,54.64,
393 52.17,53.02,57.12,53.50,53.16,49.21,52.21,40.96,
394 49.01,47.94,49.89,53.83,52.96,50.33,51.72,46.99,
395 39.06,47.99,47.91,52.35,48.51,47.39,50.45,43.66,
396 25.62,35.76,42.76,46.51,45.62,46.79,48.76,41.00,
397 52.65,55.57,57.67,56.79,55.15,54.74,50.31,47.49,
398 53.72,55.62,55.89,58.11,54.46,52.92,49.61,44.68,
399 51.59,57.44,56.50,55.12,57.22,54.61,49.92,45.20};
400}
401
402void VPPTest::testVPPIntrinsicValue() {
403
404 BOOST_TEST_MESSAGE("Testing VPP step condition...");
405
406 using namespace vpp_test;
407
408 const Date today = Date(18, December, 2011);
409 const DayCounter dc = ActualActual(ActualActual::ISDA);
410 Settings::instance().evaluationDate() = today;
411
412 // vpp parameters
413 const Real pMin = 8;
414 const Real pMax = 40;
415 const Size tMinUp = 2;
416 const Size tMinDown = 2;
417 const Real startUpFuel = 20;
418 const Real startUpFixCost = 100;
419 const Real fuelCostAddon = 3.0;
420
421 const ext::shared_ptr<SwingExercise> exercise(new SwingExercise(today, today + 6, 3600U));
422
423 // Expected values are calculated using mixed integer programming
424 // based on the gnu linear programming toolkit. For details please see:
425 // http://spanderen.de/
426 // 2011/06/23/vpp-pricing-ii-mixed-integer-linear-programming/
427 const Real efficiency[] = { 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.9 };
428 const Real expected[] = { 0.0, 2056.04, 11145.577778, 26452.04,
429 44512.461818, 62000.626667, 137591.911111};
430
431 for (Size i=0; i < LENGTH(efficiency); ++i) {
432 const Real heatRate = 1.0/efficiency[i];
433
434 VanillaVPPOption option(heatRate, pMin, pMax, tMinUp, tMinDown,
435 startUpFuel, startUpFixCost, exercise);
436
437 option.setPricingEngine(ext::shared_ptr<PricingEngine>(
438 new DynProgVPPIntrinsicValueEngine(fuelPrices, powerPrices,
439 fuelCostAddon, flatRate(forward: 0.0, dc))));
440
441 const Real calculated = option.NPV();
442
443 if (std::fabs(x: expected[i] - calculated) > 1e-4) {
444 BOOST_ERROR("Failed to reproduce reference values"
445 << "\n calculated: " << calculated
446 << "\n expected: " << expected[i]);
447
448 }
449 }
450}
451
452namespace vpp_test {
453
454 class PathFuelPrice : public FdmInnerValueCalculator {
455 public:
456 typedef FdSimpleKlugeExtOUVPPEngine::Shape Shape;
457
458 PathFuelPrice(const MultiPathGenerator<PseudoRandom>::sample_type::value_type& path,
459 ext::shared_ptr<Shape> shape)
460 : path_(path), shape_(std::move(shape)) {}
461 Real innerValue(const FdmLinearOpIterator&, Time t) override {
462 QL_REQUIRE(t-std::sqrt(QL_EPSILON) <= shape_->back().first,
463 "invalid time");
464
465 const Size i = Size(t * 365U * 24U);
466 const Real f = std::lower_bound(first: shape_->begin(), last: shape_->end(),
467 val: std::pair<Time, Real>(t-std::sqrt(QL_EPSILON), 0.0))->second;
468
469 return std::exp(x: path_[2][i] + f);
470 }
471 Real avgInnerValue(const FdmLinearOpIterator& iter, Time t) override {
472 return innerValue(iter, t);
473 }
474
475 private:
476 const MultiPathGenerator<PseudoRandom>::sample_type::value_type& path_;
477 const ext::shared_ptr<Shape> shape_;
478 };
479
480 class PathSparkSpreadPrice : public FdmInnerValueCalculator {
481 public:
482 typedef FdSimpleKlugeExtOUVPPEngine::Shape Shape;
483
484 PathSparkSpreadPrice(Real heatRate,
485 const MultiPathGenerator<PseudoRandom>::sample_type::value_type& path,
486 ext::shared_ptr<Shape> fuelShape,
487 ext::shared_ptr<Shape> powerShape)
488 : heatRate_(heatRate), path_(path), fuelShape_(std::move(fuelShape)),
489 powerShape_(std::move(powerShape)) {}
490
491 Real innerValue(const FdmLinearOpIterator&, Time t) override {
492 QL_REQUIRE(t-std::sqrt(QL_EPSILON) <= powerShape_->back().first,
493 "invalid time");
494
495 const Size i = Size(t * 365U * 24U);
496 const Real f = std::lower_bound(
497 first: powerShape_->begin(), last: powerShape_->end(),
498 val: std::pair<Time, Real>(t-std::sqrt(QL_EPSILON), 0.0))->second;
499 const Real g = std::lower_bound(
500 first: fuelShape_->begin(),last: fuelShape_->end(),
501 val: std::pair<Time, Real>(t-std::sqrt(QL_EPSILON), 0.0))->second;
502
503 return std::exp(x: f + path_[0][i]+path_[1][i])
504 - heatRate_*std::exp(x: g + path_[2][i]);
505 }
506 Real avgInnerValue(const FdmLinearOpIterator& iter, Time t) override {
507 return innerValue(iter, t);
508 }
509
510 private:
511 const Real heatRate_;
512 const MultiPathGenerator<PseudoRandom>::sample_type::value_type& path_;
513 const ext::shared_ptr<Shape> fuelShape_;
514 const ext::shared_ptr<Shape> powerShape_;
515 };
516}
517
518
519namespace vpp_test {
520 ext::shared_ptr<KlugeExtOUProcess> createKlugeExtOUProcess() {
521 // model definition
522 const Real beta = 200;
523 const Real eta = 1.0/0.2;
524 const Real lambda = 4.0;
525 const Real alpha = 7.0;
526 const Real volatility_x = 1.4;
527 const Real kappa = 4.45;
528 const Real volatility_u = std::sqrt(x: 1.3);
529 const Real rho = 0.7;
530
531 Array x0(2);
532 x0[0] = 0.0; x0[1] = 0.0;
533
534 const ext::shared_ptr<ExtendedOrnsteinUhlenbeckProcess> ouProcess(
535 new ExtendedOrnsteinUhlenbeckProcess(alpha, volatility_x, x0[0],
536 constant_b(b: x0[0])));
537 const ext::shared_ptr<ExtOUWithJumpsProcess> lnPowerProcess(
538 new ExtOUWithJumpsProcess(ouProcess, x0[1], beta, lambda, eta));
539
540 const Real u=0.0;
541 const ext::shared_ptr<ExtendedOrnsteinUhlenbeckProcess> lnGasProcess(
542 new ExtendedOrnsteinUhlenbeckProcess(kappa, volatility_u, u,
543 constant_b(b: u)));
544
545 ext::shared_ptr<KlugeExtOUProcess> klugeOUProcess(
546 new KlugeExtOUProcess(rho, lnPowerProcess, lnGasProcess));
547
548 return klugeOUProcess;
549 }
550}
551
552
553void VPPTest::testVPPPricing() {
554 BOOST_TEST_MESSAGE("Testing VPP pricing using perfect foresight or FDM...");
555
556 using namespace vpp_test;
557
558 const Date today = Date(18, December, 2011);
559 const DayCounter dc = ActualActual(ActualActual::ISDA);
560 Settings::instance().evaluationDate() = today;
561
562 // vpp parameter
563 const Real heatRate = 2.5;
564 const Real pMin = 8;
565 const Real pMax = 40;
566 const Size tMinUp = 6;
567 const Size tMinDown = 2;
568 const Real startUpFuel = 20;
569 const Real startUpFixCost = 100;
570
571 const ext::shared_ptr<SwingExercise> exercise(new SwingExercise(today, today + 6, 3600U));
572
573 VanillaVPPOption vppOption(heatRate, pMin, pMax, tMinUp, tMinDown,
574 startUpFuel, startUpFixCost, exercise);
575
576 const ext::shared_ptr<KlugeExtOUProcess> klugeOUProcess
577 = createKlugeExtOUProcess();
578 const ext::shared_ptr<ExtOUWithJumpsProcess> lnPowerProcess
579 = klugeOUProcess->getKlugeProcess();
580 const ext::shared_ptr<ExtendedOrnsteinUhlenbeckProcess> ouProcess
581 = lnPowerProcess->getExtendedOrnsteinUhlenbeckProcess();
582 const ext::shared_ptr<ExtendedOrnsteinUhlenbeckProcess> lnGasProcess
583 = klugeOUProcess->getExtOUProcess();
584
585 const Real beta = lnPowerProcess->beta();
586 const Real eta = lnPowerProcess->eta();
587 const Real lambda = lnPowerProcess->jumpIntensity();
588 const Real alpha = ouProcess->speed();
589 const Real volatility_x = ouProcess->volatility();
590 const Real kappa = lnGasProcess->speed();
591 const Real volatility_u = lnGasProcess->volatility();
592
593 const Rate irRate = 0.00;
594 const Real fuelCostAddon = 3.0;
595
596 const ext::shared_ptr<YieldTermStructure> rTS
597 = flatRate(today, forward: irRate, dc);
598
599 const Size nHours = powerPrices.size();
600
601 typedef FdSimpleKlugeExtOUVPPEngine::Shape Shape;
602 ext::shared_ptr<Shape> fuelShape(new Shape(nHours));
603 ext::shared_ptr<Shape> powerShape(new Shape(nHours));
604
605 for (Size i=0; i < nHours; ++i) {
606 const Time t = (i+1)/(365*24.);
607
608 const Real fuelPrice = fuelPrices[i];
609 const Real gs = std::log(x: fuelPrice)-squared(x: volatility_u)
610 /(4*kappa)*(1-std::exp(x: -2*kappa*t));
611 (*fuelShape)[i] = Shape::value_type(t, gs);
612
613 const Real powerPrice = powerPrices[i];
614 const Real ps = std::log(x: powerPrice)-squared(x: volatility_x)
615 /(4*alpha)*(1-std::exp(x: -2*alpha*t))
616 -lambda/beta*std::log(x: (eta-std::exp(x: -beta*t))/(eta-1.0));
617
618 (*powerShape)[i] = Shape::value_type(t, ps);
619 }
620
621 // Test: intrinsic value
622 vppOption.setPricingEngine(ext::shared_ptr<PricingEngine>(
623 new DynProgVPPIntrinsicValueEngine(fuelPrices, powerPrices,
624 fuelCostAddon, flatRate(forward: 0.0, dc))));
625
626 const Real intrinsic = vppOption.NPV();
627 const Real expectedIntrinsic = 2056.04;
628 if (std::fabs(x: intrinsic - expectedIntrinsic) > 0.1) {
629 BOOST_ERROR("Failed to reproduce intrinsic value"
630 << "\n calculated: " << intrinsic
631 << "\n expected : " << expectedIntrinsic);
632 }
633
634 // Test: finite difference price
635 const ext::shared_ptr<PricingEngine> engine(
636 new FdSimpleKlugeExtOUVPPEngine(klugeOUProcess, rTS,
637 fuelShape, powerShape, fuelCostAddon,
638 1, 25, 11, 10));
639
640 vppOption.setPricingEngine(engine);
641
642 const Real fdmPrice = vppOption.NPV();
643 const Real expectedFdmPrice = 5217.68;
644 if (std::fabs(x: fdmPrice - expectedFdmPrice) > 0.1) {
645 BOOST_ERROR("Failed to reproduce finite difference price"
646 << "\n calculated: " << fdmPrice
647 << "\n expected : " << expectedFdmPrice);
648 }
649
650 // Test: Monte-Carlo perfect foresight price
651 VanillaVPPOption::arguments args;
652 vppOption.setupArguments(&args);
653
654 const FdmVPPStepConditionFactory stepConditionFactory(args);
655
656 const ext::shared_ptr<FdmMesher> oneDimMesher(new FdmMesherComposite(
657 stepConditionFactory.stateMesher()));
658 const Size nStates = oneDimMesher->layout()->dim()[0];
659
660 const FdmVPPStepConditionMesher vppMesh = {.stateDirection: 0U, .mesher: oneDimMesher};
661
662 const TimeGrid grid(dc.yearFraction(d1: today, d2: exercise->lastDate()+1),
663 exercise->dates().size());
664 typedef PseudoRandom::rsg_type rsg_type;
665 typedef MultiPathGenerator<rsg_type>::sample_type sample_type;
666
667 rsg_type rsg = PseudoRandom::make_sequence_generator(
668 dimension: klugeOUProcess->factors() * (grid.size() - 1), seed: 1234UL);
669 MultiPathGenerator<rsg_type> generator(klugeOUProcess, grid, rsg, false);
670
671 GeneralStatistics npv;
672 const Size nTrails = 2500;
673
674 for (Size i=0; i < nTrails; ++i) {
675 const sample_type& path = generator.next();
676 const ext::shared_ptr<FdmVPPStepCondition> stepCondition(
677 stepConditionFactory.build(
678 mesh: vppMesh, fuelCostAddon,
679 fuel: ext::shared_ptr<FdmInnerValueCalculator>(
680 new PathFuelPrice(path.value, fuelShape)),
681 spark: ext::shared_ptr<FdmInnerValueCalculator>(
682 new PathSparkSpreadPrice(heatRate, path.value,
683 fuelShape, powerShape))));
684
685 Array state(nStates, 0.0);
686 for (Size j=exercise->dates().size(); j > 0; --j) {
687 stepCondition->applyTo(a&: state, t: grid.at(i: j));
688 state*=rTS->discount(t: grid.at(i: j))/rTS->discount(t: grid.at(i: j-1));
689 }
690
691 npv.add(value: state.back());
692 }
693 Real npvMC = npv.mean();
694 Real errorMC = npv.errorEstimate();
695 const Real expectedMC = 5250.0;
696 if (std::fabs(x: npvMC-expectedMC) > 3*errorMC) {
697 BOOST_ERROR("Failed to reproduce Monte-Carlo price"
698 << "\n calculated: " << npvMC
699 << "\n error ; " << errorMC
700 << "\n expected : " << expectedMC);
701 }
702 npv.reset();
703
704 // Test: Longstaff Schwartz least squares Monte-Carlo
705 // implementation is not strictly correct but saves some coding
706 const Size nCalibrationTrails = 1000U;
707 std::vector<sample_type> calibrationPaths;
708 std::vector<ext::shared_ptr<FdmVPPStepCondition> > stepConditions;
709 std::vector<ext::shared_ptr<FdmInnerValueCalculator> > sparkSpreads;
710
711 sparkSpreads.reserve(n: nCalibrationTrails);
712 stepConditions.reserve(n: nCalibrationTrails);
713 calibrationPaths.reserve(n: nCalibrationTrails);
714
715 for (Size i=0; i < nCalibrationTrails; ++i) {
716 calibrationPaths.push_back(x: generator.next());
717
718 sparkSpreads.push_back(x: ext::shared_ptr<FdmInnerValueCalculator>(
719 new PathSparkSpreadPrice(heatRate, calibrationPaths.back().value,
720 fuelShape, powerShape)));
721 stepConditions.push_back(x: stepConditionFactory.build(
722 mesh: vppMesh, fuelCostAddon,
723 fuel: ext::shared_ptr<FdmInnerValueCalculator>(
724 new PathFuelPrice(calibrationPaths.back().value, fuelShape)),
725 spark: sparkSpreads.back()));
726 }
727
728
729 const FdmLinearOpIterator iter = oneDimMesher->layout()->begin();
730
731 // prices of all calibration paths for all states
732 std::vector<Array> prices(nCalibrationTrails, Array(nStates, 0.0));
733
734 // regression coefficients for all states and all exercise dates
735 std::vector<std::vector<Array> > coeff(
736 nStates, std::vector<Array>(exercise->dates().size(), Array()));
737
738 // regression functions
739 const Size dim = 1U;
740 std::vector<ext::function<Real(Array)> > v(
741 LsmBasisSystem::multiPathBasisSystem(dim, order: 5U, type: LsmBasisSystem::Monomial));
742
743 for (Size i = exercise->dates().size(); i > 0U; --i) {
744 const Time t = grid.at(i);
745
746 std::vector<Array> x(nCalibrationTrails, Array(dim));
747
748 for (Size j=0; j < nCalibrationTrails; ++j) {
749 x[j][0] = sparkSpreads[j]->innerValue(iter, t);
750 }
751
752 for (Size k=0; k < nStates; ++k) {
753 std::vector<Real> y(nCalibrationTrails);
754
755 for (Size j=0; j < nCalibrationTrails; ++j) {
756 y[j] = prices[j][k];
757 }
758 coeff[k][i-1] = GeneralLinearLeastSquares(x, y, v).coefficients();
759
760 for (Size j=0; j < nCalibrationTrails; ++j) {
761 prices[j][k] = 0.0;
762 for (Size l=0; l < v.size(); ++l) {
763 prices[j][k] += coeff[k][i-1][l]*v[l](x[j]);
764 }
765 }
766 }
767
768 for (Size j=0; j < nCalibrationTrails; ++j) {
769 stepConditions[j]->applyTo(a&: prices[j], t: grid.at(i));
770 }
771 }
772
773 Real tmpValue = 0.0;
774 for (Size i=0; i < nTrails; ++i) {
775 Array x(dim), state(nStates, 0.0), contState(nStates, 0.0);
776
777 const sample_type& path = (i % 2) != 0U ? generator.antithetic() : generator.next();
778
779 const ext::shared_ptr<FdmInnerValueCalculator> fuelPrices(
780 new PathFuelPrice(path.value, fuelShape));
781
782 const ext::shared_ptr<FdmInnerValueCalculator> sparkSpreads(
783 new PathSparkSpreadPrice(heatRate, path.value,
784 fuelShape, powerShape));
785
786 for (Size j = exercise->dates().size(); j > 0U; --j) {
787 const Time t = grid.at(i: j);
788 const Real fuelPrice = fuelPrices->innerValue(iter, t);
789 const Real sparkSpread = sparkSpreads->innerValue(iter, t);
790 const Real startUpCost
791 = startUpFixCost + (fuelPrice + fuelCostAddon)*startUpFuel;
792
793 x[0] = sparkSpread;
794 for (Size k=0; k < nStates; ++k) {
795 contState[k] = 0.0;
796 for (Size l=0; l < v.size(); ++l) {
797 contState[k] += coeff[k][j-1][l]*v[l](x);
798 }
799 }
800
801 const Real pMinFlow = pMin*(sparkSpread - heatRate*fuelCostAddon);
802 const Real pMaxFlow = pMax*(sparkSpread - heatRate*fuelCostAddon);
803
804 // rollback continuation states and the path states
805 for (Size i=0; i < 2*tMinUp; ++i) {
806 if (i < tMinUp) {
807 state[i] += pMinFlow;
808 contState[i]+= pMinFlow;
809 }
810 else {
811 state[i] += pMaxFlow;
812 contState[i]+= pMaxFlow;
813 }
814 }
815
816 // dynamic programming using the continuation values
817 Array retVal(nStates);
818 for (Size i=0; i < tMinUp-1; ++i) {
819 retVal[i] = retVal[tMinUp + i]
820 = (contState[i+1] > contState[tMinUp + i+1])?
821 state[i+1] : state[tMinUp + i+1];
822 }
823
824 if (contState[2*tMinUp] >=
825 std::max(a: contState[tMinUp-1], b: contState[2*tMinUp-1])) {
826 retVal[tMinUp-1] = retVal[2*tMinUp-1] = state[2*tMinUp];
827 }
828 else if (contState[tMinUp-1] >= contState[2*tMinUp-1]) {
829 retVal[tMinUp-1] = retVal[2*tMinUp-1] = state[tMinUp-1];
830 }
831 else {
832 retVal[tMinUp-1] = retVal[2*tMinUp-1] = state[2*tMinUp-1];
833 }
834
835 for (Size i=0; i < tMinDown-1; ++i) {
836 retVal[2*tMinUp + i] = state[2*tMinUp + i+1];
837 }
838
839 if (contState.back() >=
840 std::max(a: contState.front(), b: contState[tMinUp]) - startUpCost) {
841 retVal.back() = state.back();
842 }
843 else if (contState.front() > contState[tMinUp]) {
844 retVal.back() = state.front()-startUpCost;
845 }
846 else {
847 retVal.back() = state[tMinUp]-startUpCost;
848 }
849 state = retVal;
850 }
851 tmpValue+=0.5*state.back();
852 if ((i % 2) != 0U) {
853 npv.add(value: tmpValue, weight: 1.0);
854 tmpValue = 0.0;
855 }
856 }
857
858 npvMC = npv.mean();
859 errorMC = npv.errorEstimate();
860 if (std::fabs(x: npvMC-fdmPrice) > 3*errorMC) {
861 BOOST_ERROR("Failed to reproduce Least Square Monte-Carlo price"
862 << "\n calculated : " << npvMC
863 << "\n error : " << errorMC
864 << "\n expected FDM : " << fdmPrice);
865 }
866}
867
868void VPPTest::testKlugeExtOUMatrixDecomposition() {
869 BOOST_TEST_MESSAGE("Testing KlugeExtOU matrix decomposition...");
870
871 using namespace vpp_test;
872
873 const Date today = Date(18, December, 2011);
874 Settings::instance().evaluationDate() = today;
875
876 const ext::shared_ptr<KlugeExtOUProcess> klugeOUProcess
877 = createKlugeExtOUProcess();
878
879 const Size xGrid = 50;
880 const Size yGrid = 20;
881 const Size uGrid = 20;
882 const Time maturity = 1;
883
884 const ext::shared_ptr<ExtOUWithJumpsProcess> klugeProcess
885 = klugeOUProcess->getKlugeProcess();
886 const ext::shared_ptr<StochasticProcess1D> ouProcess
887 = klugeProcess->getExtendedOrnsteinUhlenbeckProcess();
888
889 const ext::shared_ptr<FdmMesher> mesher(
890 new FdmMesherComposite(
891 ext::shared_ptr<Fdm1dMesher>(
892 new FdmSimpleProcess1dMesher(xGrid, ouProcess, maturity)),
893 ext::shared_ptr<Fdm1dMesher>(
894 new ExponentialJump1dMesher(yGrid,
895 klugeProcess->beta(),
896 klugeProcess->jumpIntensity(),
897 klugeProcess->eta())),
898 ext::shared_ptr<Fdm1dMesher>(
899 new FdmSimpleProcess1dMesher(uGrid,
900 klugeOUProcess->getExtOUProcess(),
901 maturity))));
902
903 const ext::shared_ptr<FdmLinearOpComposite> op(
904 new FdmKlugeExtOUOp(mesher, klugeOUProcess,
905 flatRate(today, forward: 0.0, dc: ActualActual(ActualActual::ISDA)),
906 FdmBoundaryConditionSet(), 16));
907 op->setTime(t1: 0.1, t2: 0.2);
908
909 Array x(mesher->layout()->size());
910
911 PseudoRandom::rng_type rng(PseudoRandom::urng_type(12345UL));
912 for (Real& i : x) {
913 i = rng.next().value;
914 }
915
916 const Real tol = 1e-9;
917 const Array applyExpected = op->apply(r: x);
918 const Array applyExpectedMixed = op->apply_mixed(r: x);
919
920 const std::vector<SparseMatrix> matrixDecomp(op->toMatrixDecomp());
921 const Array applyCalculated = prod(A: op->toMatrix(), x);
922 const Array applyCalculatedMixed = prod(A: matrixDecomp.back(), x);
923
924 for (Size i=0; i < x.size(); ++i) {
925 const Real diffApply = std::fabs(x: applyExpected[i]-applyCalculated[i]);
926 if (diffApply > tol && diffApply > std::fabs(x: applyExpected[i])*tol) {
927 BOOST_ERROR("Failed to reproduce apply operation" <<
928 "\n expected : " << applyExpected[i] <<
929 "\n calculated: " << applyCalculated[i] <<
930 "\n diff : " << diffApply);
931 }
932
933 const Real diffMixed = std::fabs(x: applyExpectedMixed[i]-applyCalculatedMixed[i]);
934 if (diffMixed > tol && diffMixed > std::fabs(x: applyExpected[i])*tol) {
935 BOOST_ERROR("Failed to reproduce apply operation" <<
936 "\n expected : " << applyExpectedMixed[i] <<
937 "\n calculated: " << applyCalculatedMixed[i] <<
938 "\n diff : " << diffMixed);
939 }
940 }
941
942
943 for (Size i=0; i < 3; ++i) {
944 const Array applyExpectedDir = op->apply_direction(direction: i, r: x);
945 const Array applyCalculatedDir = prod(A: matrixDecomp[i], x);
946
947 for (Size j=0; j < x.size(); ++j) {
948 const Real diff
949 = std::fabs(x: (applyExpectedDir[j] - applyCalculatedDir[j]));
950
951 if (diff > tol && diff > std::fabs(x: applyExpectedDir[j]*tol)) {
952 BOOST_ERROR("Failed to reproduce apply operation" <<
953 "\n expected : " << applyExpectedDir[i] <<
954 "\n calculated: " << applyCalculatedDir[i] <<
955 "\n diff : " << diff);
956 }
957 }
958 }
959}
960
961
962test_suite* VPPTest::suite(SpeedLevel speed) {
963 auto* suite = BOOST_TEST_SUITE("VPP Test");
964
965 suite->add(QUANTLIB_TEST_CASE(&VPPTest::testGemanRoncoroniProcess));
966 suite->add(QUANTLIB_TEST_CASE(&VPPTest::testSimpleExtOUStorageEngine));
967 suite->add(QUANTLIB_TEST_CASE(&VPPTest::testKlugeExtOUSpreadOption));
968 suite->add(QUANTLIB_TEST_CASE(&VPPTest::testVPPIntrinsicValue));
969 suite->add(QUANTLIB_TEST_CASE(
970 &VPPTest::testKlugeExtOUMatrixDecomposition));
971
972 if (speed == Slow) {
973 suite->add(QUANTLIB_TEST_CASE(&VPPTest::testVPPPricing));
974 }
975
976 return suite;
977}
978
979
980

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