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) 2005, 2006 StatPro Italia srl
6	Copyright (C) 2007 Giorgio Facchinetti
7	Copyright (C) 2009 Dimitri Reiswich
8	Copyright (C) 2014 Peter Caspers
9	Copyright (C) 2018 Klaus Spanderen
10
11	This file is part of QuantLib, a free-software/open-source library
12	for financial quantitative analysts and developers - http://quantlib.org/
13
14	QuantLib is free software: you can redistribute it and/or modify it
15	under the terms of the QuantLib license. You should have received a
16	copy of the license along with this program; if not, please email
17	<quantlib-dev@lists.sf.net>. The license is also available online at
18	<http://quantlib.org/license.shtml>.
19
20	This program is distributed in the hope that it will be useful, but WITHOUT
21	ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
22	FOR A PARTICULAR PURPOSE. See the license for more details.
23	*/
24
25	#include "interpolations.hpp"
26	#include "utilities.hpp"
27	#include <ql/experimental/volatility/noarbsabrinterpolation.hpp>
28	#include <ql/math/bspline.hpp>
29	#include <ql/math/functional.hpp>
30	#include <ql/math/integrals/simpsonintegral.hpp>
31	#include <ql/math/interpolations/backwardflatinterpolation.hpp>
32	#include <ql/math/interpolations/bicubicsplineinterpolation.hpp>
33	#include <ql/math/interpolations/chebyshevinterpolation.hpp>
34	#include <ql/math/interpolations/cubicinterpolation.hpp>
35	#include <ql/math/interpolations/forwardflatinterpolation.hpp>
36	#include <ql/math/interpolations/kernelinterpolation.hpp>
37	#include <ql/math/interpolations/kernelinterpolation2d.hpp>
38	#include <ql/math/interpolations/lagrangeinterpolation.hpp>
39	#include <ql/math/interpolations/linearinterpolation.hpp>
40	#include <ql/math/interpolations/multicubicspline.hpp>
41	#include <ql/math/interpolations/sabrinterpolation.hpp>
42	#include <ql/math/kernelfunctions.hpp>
43	#include <ql/math/optimization/levenbergmarquardt.hpp>
44	#include <ql/math/randomnumbers/sobolrsg.hpp>
45	#include <ql/math/richardsonextrapolation.hpp>
46	#include <ql/tuple.hpp>
47	#include <ql/utilities/dataformatters.hpp>
48	#include <ql/utilities/null.hpp>
49	#include <cmath>
50	#include <utility>
51
52	using namespace QuantLib;
53	using namespace boost::unit_test_framework;
54
55	namespace {
56
57	std::vector<Real> xRange(Real start, Real finish, Size points) {
58	std::vector<Real> x(points);
59	Real dx = (finish-start)/(points-1);
60	for (Size i=0; i<points-1; i++)
61	x[i] = start+i*dx;
62	x[points-1] = finish;
63	return x;
64	}
65
66	std::vector<Real> gaussian(const std::vector<Real>& x) {
67	std::vector<Real> y(x.size());
68	for (Size i=0; i<x.size(); i++)
69	y[i] = std::exp(x: -x[i]*x[i]);
70	return y;
71	}
72
73	std::vector<Real> parabolic(const std::vector<Real>& x) {
74	std::vector<Real> y(x.size());
75	for (Size i=0; i<x.size(); i++)
76	y[i] = -x[i]*x[i];
77	return y;
78	}
79
80	template <class I, class J>
81	void checkValues(const char* type,
82	const CubicInterpolation& cubic,
83	I xBegin, I xEnd, J yBegin) {
84	Real tolerance = 2.0e-15;
85	while (xBegin != xEnd) {
86	Real interpolated = cubic(*xBegin);
87	if (std::fabs(interpolated-*yBegin) > tolerance) {
88	BOOST_ERROR(type << " interpolation failed at x = " << *xBegin
89	<< std::scientific
90	<< "\n interpolated value: " << interpolated
91	<< "\n expected value: " << *yBegin
92	<< "\n error: "
93	<< std::fabs(interpolated-*yBegin));
94	}
95	++xBegin; ++yBegin;
96	}
97	}
98
99	void check1stDerivativeValue(const char* type,
100	const CubicInterpolation& cubic,
101	Real x,
102	Real value) {
103	Real tolerance = 1.0e-14;
104	Real interpolated = cubic.derivative(x);
105	Real error = std::fabs(x: interpolated-value);
106	if (error > tolerance) {
107	BOOST_ERROR(type << " interpolation first derivative failure\n"
108	<< "at x = " << x
109	<< "\n interpolated value: " << interpolated
110	<< "\n expected value: " << value
111	<< std::scientific
112	<< "\n error: " << error);
113	}
114	}
115
116	void check2ndDerivativeValue(const char* type,
117	const CubicInterpolation& cubic,
118	Real x,
119	Real value) {
120	Real tolerance = 1.0e-13;
121	Real interpolated = cubic.secondDerivative(x);
122	Real error = std::fabs(x: interpolated-value);
123	if (error > tolerance) {
124	BOOST_ERROR(type << " interpolation second derivative failure\n"
125	<< "at x = " << x
126	<< "\n interpolated value: " << interpolated
127	<< "\n expected value: " << value
128	<< std::scientific
129	<< "\n error: " << error);
130	}
131	}
132
133	void checkNotAKnotCondition(const char* type,
134	const CubicInterpolation& cubic) {
135	Real tolerance = 1.0e-14;
136	const std::vector<Real>& c = cubic.cCoefficients();
137	if (std::fabs(x: c[0]-c[1]) > tolerance) {
138	BOOST_ERROR(type << " interpolation failure"
139	<< "\n cubic coefficient of the first"
140	<< " polinomial is " << c[0]
141	<< "\n cubic coefficient of the second"
142	<< " polinomial is " << c[1]);
143	}
144	Size n = c.size();
145	if (std::fabs(x: c[n-2]-c[n-1]) > tolerance) {
146	BOOST_ERROR(type << " interpolation failure"
147	<< "\n cubic coefficient of the 2nd to last"
148	<< " polinomial is " << c[n-2]
149	<< "\n cubic coefficient of the last"
150	<< " polinomial is " << c[n-1]);
151	}
152	}
153
154	void checkSymmetry(const char* type,
155	const CubicInterpolation& cubic,
156	Real xMin) {
157	Real tolerance = 1.0e-15;
158	for (Real x = xMin; x < 0.0; x += 0.1) {
159	Real y1 = cubic(x), y2 = cubic(-x);
160	if (std::fabs(x: y1-y2) > tolerance) {
161	BOOST_ERROR(type << " interpolation not symmetric"
162	<< "\n x = " << x
163	<< "\n g(x) = " << y1
164	<< "\n g(-x) = " << y2
165	<< "\n error: " << std::fabs(y1-y2));
166	}
167	}
168	}
169
170	template <class F>
171	class errorFunction {
172	public:
173	errorFunction(F f) : f_(std::move(f)) {}
174	Real operator()(Real x) const {
175	Real temp = f_(x)-std::exp(x: -x*x);
176	return temp*temp;
177	}
178	private:
179	F f_;
180	};
181
182	template <class F>
183	errorFunction<F> make_error_function(const F& f) {
184	return errorFunction<F>(f);
185	}
186
187	Real multif(Real s, Real t, Real u, Real v, Real w) {
188	return std::sqrt(x: s * std::sinh(x: std::log(x: t)) +
189	std::exp(x: std::sin(x: u) * std::sin(x: 3 * v)) +
190	std::sinh(x: std::log(x: v * w)));
191	}
192
193	Real epanechnikovKernel(Real u){
194
195	if(std::fabs(x: u)<=1){
196	return (3.0/4.0)*(1-u*u);
197	}else{
198	return 0.0;
199	}
200	}
201
202	}
203
204
205	/* See J. M. Hyman, "Accurate monotonicity preserving cubic interpolation"
206	SIAM J. of Scientific and Statistical Computing, v. 4, 1983, pp. 645-654.
207	http://math.lanl.gov/~mac/papers/numerics/H83.pdf
208	*/
209	void InterpolationTest::testSplineErrorOnGaussianValues() {
210
211	BOOST_TEST_MESSAGE("Testing spline approximation on Gaussian data sets...");
212
213	Size points[] = { 5, 9, 17, 33 };
214
215	// complete spline data from the original 1983 Hyman paper
216	Real tabulatedErrors[] = { 3.5e-2, 2.0e-3, 4.0e-5, 1.8e-6 };
217	Real toleranceOnTabErr[] = { 0.1e-2, 0.1e-3, 0.1e-5, 0.1e-6 };
218
219	// (complete) MC spline data from the original 1983 Hyman paper
220	// NB: with the improved Hyman filter from the Dougherty, Edelman, and
221	// Hyman 1989 paper the n=17 nonmonotonicity is not filtered anymore
222	// so the error agrees with the non MC method.
223	Real tabulatedMCErrors[] = { 1.7e-2, 2.0e-3, 4.0e-5, 1.8e-6 };
224	Real toleranceOnTabMCErr[] = { 0.1e-2, 0.1e-3, 0.1e-5, 0.1e-6 };
225
226	SimpsonIntegral integral(1e-12, 10000);
227	std::vector<Real> x, y;
228
229	// still unexplained scale factor needed to obtain the numerical
230	// results from the paper
231	Real scaleFactor = 1.9;
232
233	for (Size i=0; i<LENGTH(points); i++) {
234	Size n = points[i];
235	std::vector<Real> x = xRange(start: -1.7, finish: 1.9, points: n);
236	std::vector<Real> y = gaussian(x);
237
238	// Not-a-knot
239	CubicInterpolation f(x.begin(), x.end(), y.begin(),
240	CubicInterpolation::Spline, false,
241	CubicInterpolation::NotAKnot, Null<Real>(),
242	CubicInterpolation::NotAKnot, Null<Real>());
243	f.update();
244	Real result = std::sqrt(x: integral(make_error_function(f), -1.7, 1.9));
245	result /= scaleFactor;
246	if (std::fabs(x: result-tabulatedErrors[i]) > toleranceOnTabErr[i])
247	BOOST_ERROR("Not-a-knot spline interpolation "
248	<< "\n sample points: " << n
249	<< "\n norm of difference: " << result
250	<< "\n it should be: " << tabulatedErrors[i]);
251
252	// MC not-a-knot
253	f = CubicInterpolation(x.begin(), x.end(), y.begin(),
254	CubicInterpolation::Spline, true,
255	CubicInterpolation::NotAKnot, Null<Real>(),
256	CubicInterpolation::NotAKnot, Null<Real>());
257	f.update();
258	result = std::sqrt(x: integral(make_error_function(f), -1.7, 1.9));
259	result /= scaleFactor;
260	if (std::fabs(x: result-tabulatedMCErrors[i]) > toleranceOnTabMCErr[i])
261	BOOST_ERROR("MC Not-a-knot spline interpolation "
262	<< "\n sample points: " << n
263	<< "\n norm of difference: " << result
264	<< "\n it should be: "
265	<< tabulatedMCErrors[i]);
266	}
267
268	}
269
270	/* See J. M. Hyman, "Accurate monotonicity preserving cubic interpolation"
271	SIAM J. of Scientific and Statistical Computing, v. 4, 1983, pp. 645-654.
272	http://math.lanl.gov/~mac/papers/numerics/H83.pdf
273	*/
274	void InterpolationTest::testSplineOnGaussianValues() {
275
276	BOOST_TEST_MESSAGE("Testing spline interpolation on a Gaussian data set...");
277
278	Real interpolated, interpolated2;
279	Size n = 5;
280
281	std::vector<Real> x(n), y(n);
282	Real x1_bad=-1.7, x2_bad=1.7;
283
284	for (Real start = -1.9, j=0; j<2; start+=0.2, j++) {
285	x = xRange(start, finish: start+3.6, points: n);
286	y = gaussian(x);
287
288	// Not-a-knot spline
289	CubicInterpolation f(x.begin(), x.end(), y.begin(),
290	CubicInterpolation::Spline, false,
291	CubicInterpolation::NotAKnot, Null<Real>(),
292	CubicInterpolation::NotAKnot, Null<Real>());
293	f.update();
294	checkValues(type: "Not-a-knot spline", cubic: f,
295	xBegin: x.begin(), xEnd: x.end(), yBegin: y.begin());
296	checkNotAKnotCondition(type: "Not-a-knot spline", cubic: f);
297	// bad performance
298	interpolated = f(x1_bad);
299	interpolated2= f(x2_bad);
300	if (interpolated>0.0 && interpolated2>0.0 ) {
301	BOOST_ERROR("Not-a-knot spline interpolation "
302	<< "bad performance unverified"
303	<< "\nat x = " << x1_bad
304	<< " interpolated value: " << interpolated
305	<< "\nat x = " << x2_bad
306	<< " interpolated value: " << interpolated
307	<< "\n at least one of them was expected to be < 0.0");
308	}
309
310	// MC not-a-knot spline
311	f = CubicInterpolation(x.begin(), x.end(), y.begin(),
312	CubicInterpolation::Spline, true,
313	CubicInterpolation::NotAKnot, Null<Real>(),
314	CubicInterpolation::NotAKnot, Null<Real>());
315	f.update();
316	checkValues(type: "MC not-a-knot spline", cubic: f,
317	xBegin: x.begin(), xEnd: x.end(), yBegin: y.begin());
318	// good performance
319	interpolated = f(x1_bad);
320	if (interpolated<0.0) {
321	BOOST_ERROR("MC not-a-knot spline interpolation "
322	<< "good performance unverified\n"
323	<< "at x = " << x1_bad
324	<< "\ninterpolated value: " << interpolated
325	<< "\nexpected value > 0.0");
326	}
327	interpolated = f(x2_bad);
328	if (interpolated<0.0) {
329	BOOST_ERROR("MC not-a-knot spline interpolation "
330	<< "good performance unverified\n"
331	<< "at x = " << x2_bad
332	<< "\ninterpolated value: " << interpolated
333	<< "\nexpected value > 0.0");
334	}
335	}
336	}
337
338
339	/* See J. M. Hyman, "Accurate monotonicity preserving cubic interpolation"
340	SIAM J. of Scientific and Statistical Computing, v. 4, 1983, pp. 645-654.
341	http://math.lanl.gov/~mac/papers/numerics/H83.pdf
342	*/
343	void InterpolationTest::testSplineOnRPN15AValues() {
344
345	BOOST_TEST_MESSAGE("Testing spline interpolation on RPN15A data set...");
346
347	const Real RPN15A_x[] = {
348	7.99, 8.09, 8.19, 8.7,
349	9.2, 10.0, 12.0, 15.0, 20.0
350	};
351	const Real RPN15A_y[] = {
352	0.0, 2.76429e-5, 4.37498e-5, 0.169183,
353	0.469428, 0.943740, 0.998636, 0.999919, 0.999994
354	};
355
356	Real interpolated;
357
358	// Natural spline
359	CubicInterpolation f = CubicInterpolation(
360	std::begin(arr: RPN15A_x), std::end(arr: RPN15A_x),
361	std::begin(arr: RPN15A_y),
362	CubicInterpolation::Spline, false,
363	CubicInterpolation::SecondDerivative, 0.0,
364	CubicInterpolation::SecondDerivative, 0.0);
365	f.update();
366	checkValues(type: "Natural spline", cubic: f,
367	xBegin: std::begin(arr: RPN15A_x), xEnd: std::end(arr: RPN15A_x), yBegin: std::begin(arr: RPN15A_y));
368	check2ndDerivativeValue(type: "Natural spline", cubic: f,
369	x: *std::begin(arr: RPN15A_x), value: 0.0);
370	check2ndDerivativeValue(type: "Natural spline", cubic: f,
371	x: *(std::end(arr: RPN15A_x)-1), value: 0.0);
372	// poor performance
373	Real x_bad = 11.0;
374	interpolated = f(x_bad);
375	if (interpolated<1.0) {
376	BOOST_ERROR("Natural spline interpolation "
377	<< "poor performance unverified\n"
378	<< "at x = " << x_bad
379	<< "\ninterpolated value: " << interpolated
380	<< "\nexpected value > 1.0");
381	}
382
383
384	// Clamped spline
385	f = CubicInterpolation(std::begin(arr: RPN15A_x), std::end(arr: RPN15A_x), std::begin(arr: RPN15A_y),
386	CubicInterpolation::Spline, false,
387	CubicInterpolation::FirstDerivative, 0.0,
388	CubicInterpolation::FirstDerivative, 0.0);
389	f.update();
390	checkValues(type: "Clamped spline", cubic: f,
391	xBegin: std::begin(arr: RPN15A_x), xEnd: std::end(arr: RPN15A_x), yBegin: std::begin(arr: RPN15A_y));
392	check1stDerivativeValue(type: "Clamped spline", cubic: f,
393	x: *std::begin(arr: RPN15A_x), value: 0.0);
394	check1stDerivativeValue(type: "Clamped spline", cubic: f,
395	x: *(std::end(arr: RPN15A_x)-1), value: 0.0);
396	// poor performance
397	interpolated = f(x_bad);
398	if (interpolated<1.0) {
399	BOOST_ERROR("Clamped spline interpolation "
400	<< "poor performance unverified\n"
401	<< "at x = " << x_bad
402	<< "\ninterpolated value: " << interpolated
403	<< "\nexpected value > 1.0");
404	}
405
406
407	// Not-a-knot spline
408	f = CubicInterpolation(std::begin(arr: RPN15A_x), std::end(arr: RPN15A_x), std::begin(arr: RPN15A_y),
409	CubicInterpolation::Spline, false,
410	CubicInterpolation::NotAKnot, Null<Real>(),
411	CubicInterpolation::NotAKnot, Null<Real>());
412	f.update();
413	checkValues(type: "Not-a-knot spline", cubic: f,
414	xBegin: std::begin(arr: RPN15A_x), xEnd: std::end(arr: RPN15A_x), yBegin: std::begin(arr: RPN15A_y));
415	checkNotAKnotCondition(type: "Not-a-knot spline", cubic: f);
416	// poor performance
417	interpolated = f(x_bad);
418	if (interpolated<1.0) {
419	BOOST_ERROR("Not-a-knot spline interpolation "
420	<< "poor performance unverified\n"
421	<< "at x = " << x_bad
422	<< "\ninterpolated value: " << interpolated
423	<< "\nexpected value > 1.0");
424	}
425
426
427	// MC natural spline values
428	f = CubicInterpolation(std::begin(arr: RPN15A_x), std::end(arr: RPN15A_x),
429	std::begin(arr: RPN15A_y),
430	CubicInterpolation::Spline, true,
431	CubicInterpolation::SecondDerivative, 0.0,
432	CubicInterpolation::SecondDerivative, 0.0);
433	f.update();
434	checkValues(type: "MC natural spline", cubic: f,
435	xBegin: std::begin(arr: RPN15A_x), xEnd: std::end(arr: RPN15A_x), yBegin: std::begin(arr: RPN15A_y));
436	// good performance
437	interpolated = f(x_bad);
438	if (interpolated>1.0) {
439	BOOST_ERROR("MC natural spline interpolation "
440	<< "good performance unverified\n"
441	<< "at x = " << x_bad
442	<< "\ninterpolated value: " << interpolated
443	<< "\nexpected value < 1.0");
444	}
445
446
447	// MC clamped spline values
448	f = CubicInterpolation(std::begin(arr: RPN15A_x), std::end(arr: RPN15A_x), std::begin(arr: RPN15A_y),
449	CubicInterpolation::Spline, true,
450	CubicInterpolation::FirstDerivative, 0.0,
451	CubicInterpolation::FirstDerivative, 0.0);
452	f.update();
453	checkValues(type: "MC clamped spline", cubic: f,
454	xBegin: std::begin(arr: RPN15A_x), xEnd: std::end(arr: RPN15A_x), yBegin: std::begin(arr: RPN15A_y));
455	check1stDerivativeValue(type: "MC clamped spline", cubic: f,
456	x: *std::begin(arr: RPN15A_x), value: 0.0);
457	check1stDerivativeValue(type: "MC clamped spline", cubic: f,
458	x: *(std::end(arr: RPN15A_x)-1), value: 0.0);
459	// good performance
460	interpolated = f(x_bad);
461	if (interpolated>1.0) {
462	BOOST_ERROR("MC clamped spline interpolation "
463	<< "good performance unverified\n"
464	<< "at x = " << x_bad
465	<< "\ninterpolated value: " << interpolated
466	<< "\nexpected value < 1.0");
467	}
468
469
470	// MC not-a-knot spline values
471	f = CubicInterpolation(std::begin(arr: RPN15A_x), std::end(arr: RPN15A_x), std::begin(arr: RPN15A_y),
472	CubicInterpolation::Spline, true,
473	CubicInterpolation::NotAKnot, Null<Real>(),
474	CubicInterpolation::NotAKnot, Null<Real>());
475	f.update();
476	checkValues(type: "MC not-a-knot spline", cubic: f,
477	xBegin: std::begin(arr: RPN15A_x), xEnd: std::end(arr: RPN15A_x), yBegin: std::begin(arr: RPN15A_y));
478	// good performance
479	interpolated = f(x_bad);
480	if (interpolated>1.0) {
481	BOOST_ERROR("MC clamped spline interpolation "
482	<< "good performance unverified\n"
483	<< "at x = " << x_bad
484	<< "\ninterpolated value: " << interpolated
485	<< "\nexpected value < 1.0");
486	}
487	}
488
489	/* Blossey, Frigyik, Farnum "A Note On CubicSpline Splines"
490	Applied Linear Algebra and Numerical Analysis AMATH 352 Lecture Notes
491	http://www.amath.washington.edu/courses/352-winter-2002/spline_note.pdf
492	*/
493	void InterpolationTest::testSplineOnGenericValues() {
494
495	BOOST_TEST_MESSAGE("Testing spline interpolation on generic values...");
496
497	const Real generic_x[] = { 0.0, 1.0, 3.0, 4.0 };
498	const Real generic_y[] = { 0.0, 0.0, 2.0, 2.0 };
499	const Real generic_natural_y2[] = { 0.0, 1.5, -1.5, 0.0 };
500
501	Real interpolated, error;
502	Size i, n = LENGTH(generic_x);
503	std::vector<Real> x35(3);
504
505	// Natural spline
506	CubicInterpolation f(std::begin(arr: generic_x), std::end(arr: generic_x),
507	std::begin(arr: generic_y),
508	CubicInterpolation::Spline, false,
509	CubicInterpolation::SecondDerivative,
510	generic_natural_y2[0],
511	CubicInterpolation::SecondDerivative,
512	generic_natural_y2[n-1]);
513	f.update();
514	checkValues(type: "Natural spline", cubic: f,
515	xBegin: std::begin(arr: generic_x), xEnd: std::end(arr: generic_x), yBegin: std::begin(arr: generic_y));
516	// cached second derivative
517	for (i=0; i<n; i++) {
518	interpolated = f.secondDerivative(x: generic_x[i]);
519	error = interpolated - generic_natural_y2[i];
520	if (std::fabs(x: error)>3e-16) {
521	BOOST_ERROR("Natural spline interpolation "
522	<< "second derivative failed at x=" << generic_x[i]
523	<< "\ninterpolated value: " << interpolated
524	<< "\nexpected value: " << generic_natural_y2[i]
525	<< "\nerror: " << error);
526	}
527	}
528	x35[1] = f(3.5);
529
530
531	// Clamped spline
532	Real y1a = 0.0, y1b = 0.0;
533	f = CubicInterpolation(std::begin(arr: generic_x), std::end(arr: generic_x), std::begin(arr: generic_y),
534	CubicInterpolation::Spline, false,
535	CubicInterpolation::FirstDerivative, y1a,
536	CubicInterpolation::FirstDerivative, y1b);
537	f.update();
538	checkValues(type: "Clamped spline", cubic: f,
539	xBegin: std::begin(arr: generic_x), xEnd: std::end(arr: generic_x), yBegin: std::begin(arr: generic_y));
540	check1stDerivativeValue(type: "Clamped spline", cubic: f,
541	x: *std::begin(arr: generic_x), value: 0.0);
542	check1stDerivativeValue(type: "Clamped spline", cubic: f,
543	x: *(std::end(arr: generic_x)-1), value: 0.0);
544	x35[0] = f(3.5);
545
546
547	// Not-a-knot spline
548	f = CubicInterpolation(std::begin(arr: generic_x), std::end(arr: generic_x), std::begin(arr: generic_y),
549	CubicInterpolation::Spline, false,
550	CubicInterpolation::NotAKnot, Null<Real>(),
551	CubicInterpolation::NotAKnot, Null<Real>());
552	f.update();
553	checkValues(type: "Not-a-knot spline", cubic: f,
554	xBegin: std::begin(arr: generic_x), xEnd: std::end(arr: generic_x), yBegin: std::begin(arr: generic_y));
555	checkNotAKnotCondition(type: "Not-a-knot spline", cubic: f);
556
557	x35[2] = f(3.5);
558
559	if (x35[0]>x35[1] || x35[1]>x35[2]) {
560	BOOST_ERROR("Spline interpolation failure"
561	<< "\nat x = " << 3.5
562	<< "\nclamped spline " << x35[0]
563	<< "\nnatural spline " << x35[1]
564	<< "\nnot-a-knot spline " << x35[2]
565	<< "\nvalues should be in increasing order");
566	}
567	}
568
569
570	void InterpolationTest::testSimmetricEndConditions() {
571
572	BOOST_TEST_MESSAGE("Testing symmetry of spline interpolation "
573	"end-conditions...");
574
575	Size n = 9;
576
577	std::vector<Real> x, y;
578	x = xRange(start: -1.8, finish: 1.8, points: n);
579	y = gaussian(x);
580
581	// Not-a-knot spline
582	CubicInterpolation f(x.begin(), x.end(), y.begin(),
583	CubicInterpolation::Spline, false,
584	CubicInterpolation::NotAKnot, Null<Real>(),
585	CubicInterpolation::NotAKnot, Null<Real>());
586	f.update();
587	checkValues(type: "Not-a-knot spline", cubic: f,
588	xBegin: x.begin(), xEnd: x.end(), yBegin: y.begin());
589	checkNotAKnotCondition(type: "Not-a-knot spline", cubic: f);
590	checkSymmetry(type: "Not-a-knot spline", cubic: f, xMin: x[0]);
591
592
593	// MC not-a-knot spline
594	f = CubicInterpolation(x.begin(), x.end(), y.begin(),
595	CubicInterpolation::Spline, true,
596	CubicInterpolation::NotAKnot, Null<Real>(),
597	CubicInterpolation::NotAKnot, Null<Real>());
598	f.update();
599	checkValues(type: "MC not-a-knot spline", cubic: f,
600	xBegin: x.begin(), xEnd: x.end(), yBegin: y.begin());
601	checkSymmetry(type: "MC not-a-knot spline", cubic: f, xMin: x[0]);
602	}
603
604
605	void InterpolationTest::testDerivativeEndConditions() {
606
607	BOOST_TEST_MESSAGE("Testing derivative end-conditions "
608	"for spline interpolation...");
609
610	Size n = 4;
611
612	std::vector<Real> x, y;
613	x = xRange(start: -2.0, finish: 2.0, points: n);
614	y = parabolic(x);
615
616	// Not-a-knot spline
617	CubicInterpolation f(x.begin(), x.end(), y.begin(),
618	CubicInterpolation::Spline, false,
619	CubicInterpolation::NotAKnot, Null<Real>(),
620	CubicInterpolation::NotAKnot, Null<Real>());
621	f.update();
622	checkValues(type: "Not-a-knot spline", cubic: f,
623	xBegin: x.begin(), xEnd: x.end(), yBegin: y.begin());
624	check1stDerivativeValue(type: "Not-a-knot spline", cubic: f,
625	x: x[0], value: 4.0);
626	check1stDerivativeValue(type: "Not-a-knot spline", cubic: f,
627	x: x[n-1], value: -4.0);
628	check2ndDerivativeValue(type: "Not-a-knot spline", cubic: f,
629	x: x[0], value: -2.0);
630	check2ndDerivativeValue(type: "Not-a-knot spline", cubic: f,
631	x: x[n-1], value: -2.0);
632
633
634	// Clamped spline
635	f = CubicInterpolation(x.begin(), x.end(), y.begin(),
636	CubicInterpolation::Spline, false,
637	CubicInterpolation::FirstDerivative, 4.0,
638	CubicInterpolation::FirstDerivative, -4.0);
639	f.update();
640	checkValues(type: "Clamped spline", cubic: f,
641	xBegin: x.begin(), xEnd: x.end(), yBegin: y.begin());
642	check1stDerivativeValue(type: "Clamped spline", cubic: f,
643	x: x[0], value: 4.0);
644	check1stDerivativeValue(type: "Clamped spline", cubic: f,
645	x: x[n-1], value: -4.0);
646	check2ndDerivativeValue(type: "Clamped spline", cubic: f,
647	x: x[0], value: -2.0);
648	check2ndDerivativeValue(type: "Clamped spline", cubic: f,
649	x: x[n-1], value: -2.0);
650
651
652	// SecondDerivative spline
653	f = CubicInterpolation(x.begin(), x.end(), y.begin(),
654	CubicInterpolation::Spline, false,
655	CubicInterpolation::SecondDerivative, -2.0,
656	CubicInterpolation::SecondDerivative, -2.0);
657	f.update();
658	checkValues(type: "SecondDerivative spline", cubic: f,
659	xBegin: x.begin(), xEnd: x.end(), yBegin: y.begin());
660	check1stDerivativeValue(type: "SecondDerivative spline", cubic: f,
661	x: x[0], value: 4.0);
662	check1stDerivativeValue(type: "SecondDerivative spline", cubic: f,
663	x: x[n-1], value: -4.0);
664	check2ndDerivativeValue(type: "SecondDerivative spline", cubic: f,
665	x: x[0], value: -2.0);
666	check2ndDerivativeValue(type: "SecondDerivative spline", cubic: f,
667	x: x[n-1], value: -2.0);
668
669	// MC Not-a-knot spline
670	f = CubicInterpolation(x.begin(), x.end(), y.begin(),
671	CubicInterpolation::Spline, true,
672	CubicInterpolation::NotAKnot, Null<Real>(),
673	CubicInterpolation::NotAKnot, Null<Real>());
674	f.update();
675	checkValues(type: "MC Not-a-knot spline", cubic: f,
676	xBegin: x.begin(), xEnd: x.end(), yBegin: y.begin());
677	check1stDerivativeValue(type: "MC Not-a-knot spline", cubic: f,
678	x: x[0], value: 4.0);
679	check1stDerivativeValue(type: "MC Not-a-knot spline", cubic: f,
680	x: x[n-1], value: -4.0);
681	check2ndDerivativeValue(type: "MC Not-a-knot spline", cubic: f,
682	x: x[0], value: -2.0);
683	check2ndDerivativeValue(type: "MC Not-a-knot spline", cubic: f,
684	x: x[n-1], value: -2.0);
685
686
687	// MC Clamped spline
688	f = CubicInterpolation(x.begin(), x.end(), y.begin(),
689	CubicInterpolation::Spline, true,
690	CubicInterpolation::FirstDerivative, 4.0,
691	CubicInterpolation::FirstDerivative, -4.0);
692	f.update();
693	checkValues(type: "MC Clamped spline", cubic: f,
694	xBegin: x.begin(), xEnd: x.end(), yBegin: y.begin());
695	check1stDerivativeValue(type: "MC Clamped spline", cubic: f,
696	x: x[0], value: 4.0);
697	check1stDerivativeValue(type: "MC Clamped spline", cubic: f,
698	x: x[n-1], value: -4.0);
699	check2ndDerivativeValue(type: "MC Clamped spline", cubic: f,
700	x: x[0], value: -2.0);
701	check2ndDerivativeValue(type: "MC Clamped spline", cubic: f,
702	x: x[n-1], value: -2.0);
703
704
705	// MC SecondDerivative spline
706	f = CubicInterpolation(x.begin(), x.end(), y.begin(),
707	CubicInterpolation::Spline, true,
708	CubicInterpolation::SecondDerivative, -2.0,
709	CubicInterpolation::SecondDerivative, -2.0);
710	f.update();
711	checkValues(type: "MC SecondDerivative spline", cubic: f,
712	xBegin: x.begin(), xEnd: x.end(), yBegin: y.begin());
713	check1stDerivativeValue(type: "MC SecondDerivative spline", cubic: f,
714	x: x[0], value: 4.0);
715	check1stDerivativeValue(type: "MC SecondDerivative spline", cubic: f,
716	x: x[n-1], value: -4.0);
717	check2ndDerivativeValue(type: "SecondDerivative spline", cubic: f,
718	x: x[0], value: -2.0);
719	check2ndDerivativeValue(type: "MC SecondDerivative spline", cubic: f,
720	x: x[n-1], value: -2.0);
721
722	}
723
724
725	/* See R. L. Dougherty, A. Edelman, J. M. Hyman,
726	"Nonnegativity-, Monotonicity-, or Convexity-Preserving CubicSpline and Quintic
727	Hermite Interpolation"
728	Mathematics Of Computation, v. 52, n. 186, April 1989, pp. 471-494.
729	*/
730	void InterpolationTest::testNonRestrictiveHymanFilter() {
731
732	BOOST_TEST_MESSAGE("Testing non-restrictive Hyman filter...");
733
734	Size n = 4;
735
736	std::vector<Real> x, y;
737	x = xRange(start: -2.0, finish: 2.0, points: n);
738	y = parabolic(x);
739	Real zero=0.0, interpolated, expected=0.0;
740
741	// MC Not-a-knot spline
742	CubicInterpolation f(x.begin(), x.end(), y.begin(),
743	CubicInterpolation::Spline, true,
744	CubicInterpolation::NotAKnot, Null<Real>(),
745	CubicInterpolation::NotAKnot, Null<Real>());
746	f.update();
747	interpolated = f(zero);
748	if (std::fabs(x: interpolated-expected)>1e-15) {
749	BOOST_ERROR("MC not-a-knot spline"
750	<< " interpolation failed at x = " << zero
751	<< "\n interpolated value: " << interpolated
752	<< "\n expected value: " << expected
753	<< "\n error: "
754	<< std::fabs(interpolated-expected));
755	}
756
757
758	// MC Clamped spline
759	f = CubicInterpolation(x.begin(), x.end(), y.begin(),
760	CubicInterpolation::Spline, true,
761	CubicInterpolation::FirstDerivative, 4.0,
762	CubicInterpolation::FirstDerivative, -4.0);
763	f.update();
764	interpolated = f(zero);
765	if (std::fabs(x: interpolated-expected)>1e-15) {
766	BOOST_ERROR("MC clamped spline"
767	<< " interpolation failed at x = " << zero
768	<< "\n interpolated value: " << interpolated
769	<< "\n expected value: " << expected
770	<< "\n error: "
771	<< std::fabs(interpolated-expected));
772	}
773
774
775	// MC SecondDerivative spline
776	f = CubicInterpolation(x.begin(), x.end(), y.begin(),
777	CubicInterpolation::Spline, true,
778	CubicInterpolation::SecondDerivative, -2.0,
779	CubicInterpolation::SecondDerivative, -2.0);
780	f.update();
781	interpolated = f(zero);
782	if (std::fabs(x: interpolated-expected)>1e-15) {
783	BOOST_ERROR("MC SecondDerivative spline"
784	<< " interpolation failed at x = " << zero
785	<< "\n interpolated value: " << interpolated
786	<< "\n expected value: " << expected
787	<< "\n error: "
788	<< std::fabs(interpolated-expected));
789	}
790
791	}
792
793	void InterpolationTest::testMultiSpline() {
794	BOOST_TEST_MESSAGE("Testing N-dimensional cubic spline...");
795
796	std::vector<Size> dim(5);
797	dim[0] = 6; dim[1] = 5; dim[2] = 5; dim[3] = 6; dim[4] = 4;
798
799	std::vector<Real> args(5), offsets(5);
800	offsets[0] = 1.005; offsets[1] = 14.0; offsets[2] = 33.005;
801	offsets[3] = 35.025; offsets[4] = 19.025;
802
803	Real &s = args[0] = offsets[0],
804	&t = args[1] = offsets[1],
805	&u = args[2] = offsets[2],
806	&v = args[3] = offsets[3],
807	&w = args[4] = offsets[4];
808
809	Size i, j, k, l, m;
810
811	SplineGrid grid(5);
812
813	Real r = 0.15;
814
815	for (i = 0; i < 5; ++i) {
816	Real temp = offsets[i];
817	for (j = 0; j < dim[i]; temp += r, ++j)
818	grid[i].push_back(x: temp);
819	}
820
821	MultiCubicSpline<5>::data_table y5(dim);
822
823	for (i = 0; i < dim[0]; ++i)
824	for (j = 0; j < dim[1]; ++j)
825	for (k = 0; k < dim[2]; ++k)
826	for (l = 0; l < dim[3]; ++l)
827	for (m = 0; m < dim[4]; ++m)
828	y5[i][j][k][l][m] =
829	multif(s: grid[0][i], t: grid[1][j], u: grid[2][k],
830	v: grid[3][l], w: grid[4][m]);
831
832	MultiCubicSpline<5> cs(grid, y5);
833	/* it would fail with
834	for (i = 0; i < dim[0]; ++i)
835	for (j = 0; j < dim[1]; ++j)
836	for (k = 0; k < dim[2]; ++k)
837	for (l = 0; l < dim[3]; ++l)
838	for (m = 0; m < dim[4]; ++m) {
839	*/
840	for (i = 1; i < dim[0]-1; ++i)
841	for (j = 1; j < dim[1]-1; ++j)
842	for (k = 1; k < dim[2]-1; ++k)
843	for (l = 1; l < dim[3]-1; ++l)
844	for (m = 1; m < dim[4]-1; ++m) {
845	s = grid[0][i];
846	t = grid[1][j];
847	u = grid[2][k];
848	v = grid[3][l];
849	w = grid[4][m];
850	Real interpolated = cs(args);
851	Real expected = y5[i][j][k][l][m];
852	Real error = std::fabs(x: interpolated-expected);
853	Real tolerance = 1e-16;
854	if (error > tolerance) {
855	BOOST_ERROR(
856	"\n At ("
857	<< s << "," << t << "," << u << ","
858	<< v << "," << w << "):"
859	<< "\n interpolated: " << interpolated
860	<< "\n actual value: " << expected
861	<< "\n error: " << error
862	<< "\n tolerance: " << tolerance);
863	}
864	}
865
866
867	unsigned long seed = 42;
868	SobolRsg rsg(5, seed);
869
870	Real tolerance = 1.7e-4;
871	// actually tested up to 2^21-1=2097151 Sobol draws
872	for (i = 0; i < 1023; ++i) {
873	const std::vector<Real>& next = rsg.nextSequence().value;
874	s = grid[0].front() + next[0]*(grid[0].back()-grid[0].front());
875	t = grid[1].front() + next[1]*(grid[1].back()-grid[1].front());
876	u = grid[2].front() + next[2]*(grid[2].back()-grid[2].front());
877	v = grid[3].front() + next[3]*(grid[3].back()-grid[3].front());
878	w = grid[4].front() + next[4]*(grid[4].back()-grid[4].front());
879	Real interpolated = cs(args), expected = multif(s, t, u, v, w);
880	Real error = std::fabs(x: interpolated-expected);
881	if (error > tolerance) {
882	BOOST_ERROR(
883	"\n At ("
884	<< s << "," << t << "," << u << "," << v << "," << w << "):"
885	<< "\n interpolated: " << interpolated
886	<< "\n actual value: " << expected
887	<< "\n error: " << error
888	<< "\n tolerance: " << tolerance);
889	}
890	}
891	}
892
893	namespace {
894
895	struct NotThrown {};
896
897	}
898
899	void InterpolationTest::testAsFunctor() {
900
901	BOOST_TEST_MESSAGE("Testing use of interpolations as functors...");
902
903	const Real x[] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
904	const Real y[] = { 5.0, 4.0, 3.0, 2.0, 1.0 };
905
906	Interpolation f = LinearInterpolation(std::begin(arr: x), std::end(arr: x), std::begin(arr: y));
907	f.update();
908
909	const Real x2[] = { -2.0, -1.0, 0.0, 1.0, 3.0, 4.0, 5.0, 6.0, 7.0 };
910	Size N = LENGTH(x2);
911	std::vector<Real> y2(N);
912	Real tolerance = 1.0e-12;
913
914	// case 1: extrapolation not allowed
915	try {
916	std::transform(first: std::begin(arr: x2), last: std::end(arr: x2), result: y2.begin(), unary_op: f);
917	throw NotThrown();
918	}
919	catch (Error&) {
920	// as expected; do nothing
921	}
922	catch (NotThrown&) {
923	QL_FAIL("failed to throw exception when trying to extrapolate");
924	}
925
926	// case 2: enable extrapolation
927	f.enableExtrapolation();
928	y2 = std::vector<Real>(N);
929	std::transform(first: std::begin(arr: x2), last: std::end(arr: x2), result: y2.begin(), unary_op: f);
930	for (Size i=0; i<N; i++) {
931	Real expected = 5.0-x2[i];
932	if (std::fabs(x: y2[i]-expected) > tolerance)
933	BOOST_ERROR(
934	"failed to reproduce " << io::ordinal(i+1) << " expected datum"
935	<< std::fixed
936	<< "\n expected: " << expected
937	<< "\n calculated: " << y2[i]
938	<< std::scientific
939	<< "\n error: " << std::fabs(y2[i]-expected));
940	}
941	}
942
943
944	namespace {
945
946	Integer sign(Real y1, Real y2) {
947	return y1 == y2 ? 0 :
948	y1 < y2 ? 1 :
949	-1 ;
950	}
951
952	}
953
954	void InterpolationTest::testFritschButland() {
955
956	BOOST_TEST_MESSAGE("Testing Fritsch-Butland interpolation...");
957
958	const Real x[5] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
959	const Real y[][5] = {{ 1.0, 2.0, 1.0, 1.0, 2.0 },
960	{ 1.0, 2.0, 1.0, 1.0, 1.0 },
961	{ 2.0, 1.0, 0.0, 2.0, 3.0 }};
962
963	for (Size i=0; i<3; ++i) {
964
965	Interpolation f = FritschButlandCubic(std::begin(arr: x), std::end(arr: x), std::begin(arr: y[i]));
966	f.update();
967
968	for (Size j=0; j<4; ++j) {
969	Real left_knot = x[j];
970	Integer expected_sign = sign(y1: y[i][j], y2: y[i][j+1]);
971	for (Size k=0; k<10; ++k) {
972	Real x1 = left_knot + k*0.1, x2 = left_knot + (k+1)*0.1;
973	Real y1 = f(x1), y2 = f(x2);
974	if (std::isnan(x: y1))
975	BOOST_ERROR("NaN detected in case " << i << ":"
976	<< std::fixed
977	<< "\n f(" << x1 << ") = " << y1);
978	else if (sign(y1, y2) != expected_sign)
979	BOOST_ERROR("interpolation is not monotonic "
980	"in case " << i << ":"
981	<< std::fixed
982	<< "\n f(" << x1 << ") = " << y1
983	<< "\n f(" << x2 << ") = " << y2);
984	}
985	}
986	}
987	}
988
989
990	void InterpolationTest::testBackwardFlat() {
991
992	BOOST_TEST_MESSAGE("Testing backward-flat interpolation...");
993
994	const Real x[] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
995	const Real y[] = { 5.0, 4.0, 3.0, 2.0, 1.0 };
996
997	Interpolation f = BackwardFlatInterpolation(std::begin(arr: x), std::end(arr: x), std::begin(arr: y));
998	f.update();
999
1000	Size N = LENGTH(x);
1001	Size i;
1002	Real tolerance = 1.0e-12;
1003
1004	// at original points
1005	for (i=0; i<N; i++) {
1006	Real p = x[i];
1007	Real calculated = f(p);
1008	Real expected = y[i];
1009	if (std::fabs(x: expected-calculated) > tolerance)
1010	BOOST_ERROR(
1011	"failed to reproduce " << io::ordinal(i+1) << " datum"
1012	<< std::fixed
1013	<< "\n expected: " << expected
1014	<< "\n calculated: " << calculated
1015	<< std::scientific
1016	<< "\n error: " << std::fabs(calculated-expected));
1017	}
1018
1019	// at middle points
1020	for (i=0; i<N-1; i++) {
1021	Real p = (x[i]+x[i+1])/2;
1022	Real calculated = f(p);
1023	Real expected = y[i+1];
1024	if (std::fabs(x: expected-calculated) > tolerance)
1025	BOOST_ERROR(
1026	"failed to interpolate correctly at " << p
1027	<< std::fixed
1028	<< "\n expected: " << expected
1029	<< "\n calculated: " << calculated
1030	<< std::scientific
1031	<< "\n error: " << std::fabs(calculated-expected));
1032	}
1033
1034	// outside the original range
1035	f.enableExtrapolation();
1036
1037	Real p = x[0] - 0.5;
1038	Real calculated = f(p);
1039	Real expected = y[0];
1040	if (std::fabs(x: expected-calculated) > tolerance)
1041	BOOST_ERROR(
1042	"failed to extrapolate correctly at " << p
1043	<< std::fixed
1044	<< "\n expected: " << expected
1045	<< "\n calculated: " << calculated
1046	<< std::scientific
1047	<< "\n error: " << std::fabs(calculated-expected));
1048
1049	p = x[N-1] + 0.5;
1050	calculated = f(p);
1051	expected = y[N-1];
1052	if (std::fabs(x: expected-calculated) > tolerance)
1053	BOOST_ERROR(
1054	"failed to extrapolate correctly at " << p
1055	<< std::fixed
1056	<< "\n expected: " << expected
1057	<< "\n calculated: " << calculated
1058	<< std::scientific
1059	<< "\n error: " << std::fabs(calculated-expected));
1060
1061	// primitive at original points
1062	calculated = f.primitive(x: x[0]);
1063	expected = 0.0;
1064	if (std::fabs(x: expected-calculated) > tolerance)
1065	BOOST_ERROR(
1066	"failed to calculate primitive at " << x[0]
1067	<< std::fixed
1068	<< "\n expected: " << expected
1069	<< "\n calculated: " << calculated
1070	<< std::scientific
1071	<< "\n error: " << std::fabs(calculated-expected));
1072
1073	Real sum = 0.0;
1074	for (i=1; i<N; i++) {
1075	sum += (x[i]-x[i-1])*y[i];
1076	Real calculated = f.primitive(x: x[i]);
1077	Real expected = sum;
1078	if (std::fabs(x: expected-calculated) > tolerance)
1079	BOOST_ERROR(
1080	"failed to calculate primitive at " << x[i]
1081	<< std::fixed
1082	<< "\n expected: " << expected
1083	<< "\n calculated: " << calculated
1084	<< std::scientific
1085	<< "\n error: " << std::fabs(calculated-expected));
1086	}
1087
1088	// primitive at middle points
1089	sum = 0.0;
1090	for (i=0; i<N-1; i++) {
1091	Real p = (x[i]+x[i+1])/2;
1092	sum += (x[i+1]-x[i])*y[i+1]/2;
1093	Real calculated = f.primitive(x: p);
1094	Real expected = sum;
1095	sum += (x[i+1]-x[i])*y[i+1]/2;
1096	if (std::fabs(x: expected-calculated) > tolerance)
1097	BOOST_ERROR(
1098	"failed to calculate primitive at " << x[i]
1099	<< std::fixed
1100	<< "\n expected: " << expected
1101	<< "\n calculated: " << calculated
1102	<< std::scientific
1103	<< "\n error: " << std::fabs(calculated-expected));
1104	}
1105
1106	}
1107
1108	void InterpolationTest::testForwardFlat() {
1109
1110	BOOST_TEST_MESSAGE("Testing forward-flat interpolation...");
1111
1112	const Real x[] = { 0.0, 1.0, 2.0, 3.0, 4.0 };
1113	const Real y[] = { 5.0, 4.0, 3.0, 2.0, 1.0 };
1114
1115	Interpolation f = ForwardFlatInterpolation(std::begin(arr: x), std::end(arr: x), std::begin(arr: y));
1116	f.update();
1117
1118	Size N = LENGTH(x);
1119	Size i;
1120	Real tolerance = 1.0e-12;
1121
1122	// at original points
1123	for (i=0; i<N; i++) {
1124	Real p = x[i];
1125	Real calculated = f(p);
1126	Real expected = y[i];
1127	if (std::fabs(x: expected-calculated) > tolerance)
1128	BOOST_ERROR(
1129	"failed to reproduce " << io::ordinal(i+1) << " datum"
1130	<< std::fixed
1131	<< "\n expected: " << expected
1132	<< "\n calculated: " << calculated
1133	<< std::scientific
1134	<< "\n error: " << std::fabs(calculated-expected));
1135	}
1136
1137	// at middle points
1138	for (i=0; i<N-1; i++) {
1139	Real p = (x[i]+x[i+1])/2;
1140	Real calculated = f(p);
1141	Real expected = y[i];
1142	if (std::fabs(x: expected-calculated) > tolerance)
1143	BOOST_ERROR(
1144	"failed to interpolate correctly at " << p
1145	<< std::fixed
1146	<< "\n expected: " << expected
1147	<< "\n calculated: " << calculated
1148	<< std::scientific
1149	<< "\n error: " << std::fabs(calculated-expected));
1150	}
1151
1152	// outside the original range
1153	f.enableExtrapolation();
1154
1155	Real p = x[0] - 0.5;
1156	Real calculated = f(p);
1157	Real expected = y[0];
1158	if (std::fabs(x: expected-calculated) > tolerance)
1159	BOOST_ERROR(
1160	"failed to extrapolate correctly at " << p
1161	<< std::fixed
1162	<< "\n expected: " << expected
1163	<< "\n calculated: " << calculated
1164	<< std::scientific
1165	<< "\n error: " << std::fabs(calculated-expected));
1166
1167	p = x[N-1] + 0.5;
1168	calculated = f(p);
1169	expected = y[N-1];
1170	if (std::fabs(x: expected-calculated) > tolerance)
1171	BOOST_ERROR(
1172	"failed to extrapolate correctly at " << p
1173	<< std::fixed
1174	<< "\n expected: " << expected
1175	<< "\n calculated: " << calculated
1176	<< std::scientific
1177	<< "\n error: " << std::fabs(calculated-expected));
1178
1179	// primitive at original points
1180	calculated = f.primitive(x: x[0]);
1181	expected = 0.0;
1182	if (std::fabs(x: expected-calculated) > tolerance)
1183	BOOST_ERROR(
1184	"failed to calculate primitive at " << x[0]
1185	<< std::fixed
1186	<< "\n expected: " << expected
1187	<< "\n calculated: " << calculated
1188	<< std::scientific
1189	<< "\n error: " << std::fabs(calculated-expected));
1190
1191	Real sum = 0.0;
1192	for (i=1; i<N; i++) {
1193	sum += (x[i]-x[i-1])*y[i-1];
1194	Real calculated = f.primitive(x: x[i]);
1195	Real expected = sum;
1196	if (std::fabs(x: expected-calculated) > tolerance)
1197	BOOST_ERROR(
1198	"failed to calculate primitive at " << x[i]
1199	<< std::fixed
1200	<< "\n expected: " << expected
1201	<< "\n calculated: " << calculated
1202	<< std::scientific
1203	<< "\n error: " << std::fabs(calculated-expected));
1204	}
1205
1206	// primitive at middle points
1207	sum = 0.0;
1208	for (i=0; i<N-1; i++) {
1209	Real p = (x[i]+x[i+1])/2;
1210	sum += (x[i+1]-x[i])*y[i]/2;
1211	Real calculated = f.primitive(x: p);
1212	Real expected = sum;
1213	sum += (x[i+1]-x[i])*y[i]/2;
1214	if (std::fabs(x: expected-calculated) > tolerance)
1215	BOOST_ERROR(
1216	"failed to calculate primitive at " << p
1217	<< std::fixed
1218	<< "\n expected: " << expected
1219	<< "\n calculated: " << calculated
1220	<< std::scientific
1221	<< "\n error: " << std::fabs(calculated-expected));
1222	}
1223	}
1224
1225	void InterpolationTest::testSabrInterpolation(){
1226
1227	BOOST_TEST_MESSAGE("Testing Sabr interpolation...");
1228
1229	// Test SABR function against input volatilities
1230	Real tolerance = 1.0e-12;
1231	std::vector<Real> strikes(31);
1232	std::vector<Real> volatilities(31);
1233	// input strikes
1234	strikes[0] = 0.03 ; strikes[1] = 0.032 ; strikes[2] = 0.034 ;
1235	strikes[3] = 0.036 ; strikes[4] = 0.038 ; strikes[5] = 0.04 ;
1236	strikes[6] = 0.042 ; strikes[7] = 0.044 ; strikes[8] = 0.046 ;
1237	strikes[9] = 0.048 ; strikes[10] = 0.05 ; strikes[11] = 0.052 ;
1238	strikes[12] = 0.054 ; strikes[13] = 0.056 ; strikes[14] = 0.058 ;
1239	strikes[15] = 0.06 ; strikes[16] = 0.062 ; strikes[17] = 0.064 ;
1240	strikes[18] = 0.066 ; strikes[19] = 0.068 ; strikes[20] = 0.07 ;
1241	strikes[21] = 0.072 ; strikes[22] = 0.074 ; strikes[23] = 0.076 ;
1242	strikes[24] = 0.078 ; strikes[25] = 0.08 ; strikes[26] = 0.082 ;
1243	strikes[27] = 0.084 ; strikes[28] = 0.086 ; strikes[29] = 0.088;
1244	strikes[30] = 0.09;
1245	// input volatilities
1246	volatilities[0] = 1.16725837321531 ; volatilities[1] = 1.15226075991385 ; volatilities[2] = 1.13829711098834 ;
1247	volatilities[3] = 1.12524190877505 ; volatilities[4] = 1.11299079244474 ; volatilities[5] = 1.10145609357162 ;
1248	volatilities[6] = 1.09056348513411 ; volatilities[7] = 1.08024942745106 ; volatilities[8] = 1.07045919457758 ;
1249	volatilities[9] = 1.06114533019077 ; volatilities[10] = 1.05226642581503 ; volatilities[11] = 1.04378614411707 ;
1250	volatilities[12] = 1.03567243073732 ; volatilities[13] = 1.0278968727451 ; volatilities[14] = 1.02043417226345 ;
1251	volatilities[15] = 1.01326171139321 ; volatilities[16] = 1.00635919013311 ; volatilities[17] = 0.999708323124949 ;
1252	volatilities[18] = 0.993292584155381 ; volatilities[19] = 0.987096989695393 ; volatilities[20] = 0.98110791455717 ;
1253	volatilities[21] = 0.975312934134512 ; volatilities[22] = 0.969700688771689 ; volatilities[23] = 0.964260766651027;
1254	volatilities[24] = 0.958983602256592 ; volatilities[25] = 0.953860388001395 ; volatilities[26] = 0.948882997029509 ;
1255	volatilities[27] = 0.944043915545469 ; volatilities[28] = 0.939336183299237 ; volatilities[29] = 0.934753341079515 ;
1256	volatilities[30] = 0.930289384251337;
1257
1258	Time expiry = 1.0;
1259	Real forward = 0.039;
1260	// input SABR coefficients (corresponding to the vols above)
1261	Real initialAlpha = 0.3;
1262	Real initialBeta = 0.6;
1263	Real initialNu = 0.02;
1264	Real initialRho = 0.01;
1265	// calculate SABR vols and compare with input vols
1266	for(Size i=0; i< strikes.size(); i++){
1267	Real calculatedVol = sabrVolatility(strike: strikes[i], forward, expiryTime: expiry,
1268	alpha: initialAlpha, beta: initialBeta,
1269	nu: initialNu, rho: initialRho);
1270	if (std::fabs(x: volatilities[i]-calculatedVol) > tolerance)
1271	BOOST_ERROR(
1272	"failed to calculate Sabr function at strike " << strikes[i]
1273	<< "\n expected: " << volatilities[i]
1274	<< "\n calculated: " << calculatedVol
1275	<< "\n error: " << std::fabs(calculatedVol-volatilities[i]));
1276	}
1277
1278	// Test SABR calibration against input parameters
1279	// Use default values (but not null, since then parameters
1280	// will then not be fixed during optimization, see the
1281	// interpolation constructor, thus rendering the test cases
1282	// with fixed parameters non-sensical)
1283	Real alphaGuess = std::sqrt(x: 0.2);
1284	Real betaGuess = 0.5;
1285	Real nuGuess = std::sqrt(x: 0.4);
1286	Real rhoGuess = 0.0;
1287
1288	const bool vegaWeighted[]= {true, false};
1289	const bool isAlphaFixed[]= {true, false};
1290	const bool isBetaFixed[]= {true, false};
1291	const bool isNuFixed[]= {true, false};
1292	const bool isRhoFixed[]= {true, false};
1293
1294	Real calibrationTolerance = 5.0e-8;
1295	// initialize optimization methods
1296	std::vector<ext::shared_ptr<OptimizationMethod>> methods_ = {
1297	ext::shared_ptr<OptimizationMethod>(new Simplex(0.01)),
1298	ext::shared_ptr<OptimizationMethod>(new LevenbergMarquardt(1e-8, 1e-8, 1e-8))
1299	};
1300	// Initialize end criteria
1301	ext::shared_ptr<EndCriteria> endCriteria(new
1302	EndCriteria(100000, 100, 1e-8, 1e-8, 1e-8));
1303	// Test looping over all possibilities
1304	for (auto& method : methods_) {
1305	for (bool i : vegaWeighted) {
1306	for (bool k_a : isAlphaFixed) {
1307	for (bool k_b : isBetaFixed) {
1308	for (bool k_n : isNuFixed) {
1309	for (bool k_r : isRhoFixed) {
1310	// to meet the tough calibration tolerance we need to lower the default
1311	// error threshold for accepting a calibration (to be more specific,
1312	// some of the new test cases arising from fixing a subset of the
1313	// model's parameters do not calibrate with the desired error using the
1314	// initial guess (i.e. optimization runs into a local minimum) - then a
1315	// series of random start values for optimization is chosen until our
1316	// tight custom error threshold is satisfied.
1317	SABRInterpolation sabrInterpolation(
1318	strikes.begin(), strikes.end(), volatilities.begin(), expiry,
1319	forward, k_a ? initialAlpha : alphaGuess,
1320	k_b ? initialBeta : betaGuess, k_n ? initialNu : nuGuess,
1321	k_r ? initialRho : rhoGuess, k_a, k_b, k_n, k_r, i, endCriteria,
1322	method, 1E-10);
1323	sabrInterpolation.update();
1324
1325	// Recover SABR calibration parameters
1326	bool failed = false;
1327	Real calibratedAlpha = sabrInterpolation.alpha();
1328	Real calibratedBeta = sabrInterpolation.beta();
1329	Real calibratedNu = sabrInterpolation.nu();
1330	Real calibratedRho = sabrInterpolation.rho();
1331	Real error;
1332
1333	// compare results: alpha
1334	error = std::fabs(x: initialAlpha - calibratedAlpha);
1335	if (error > calibrationTolerance) {
1336	BOOST_ERROR("\nfailed to calibrate alpha Sabr parameter:"
1337	<< "\n expected: " << initialAlpha
1338	<< "\n calibrated: " << calibratedAlpha
1339	<< "\n error: " << error);
1340	failed = true;
1341	}
1342	// Beta
1343	error = std::fabs(x: initialBeta - calibratedBeta);
1344	if (error > calibrationTolerance) {
1345	BOOST_ERROR("\nfailed to calibrate beta Sabr parameter:"
1346	<< "\n expected: " << initialBeta
1347	<< "\n calibrated: " << calibratedBeta
1348	<< "\n error: " << error);
1349	failed = true;
1350	}
1351	// Nu
1352	error = std::fabs(x: initialNu - calibratedNu);
1353	if (error > calibrationTolerance) {
1354	BOOST_ERROR("\nfailed to calibrate nu Sabr parameter:"
1355	<< "\n expected: " << initialNu
1356	<< "\n calibrated: " << calibratedNu
1357	<< "\n error: " << error);
1358	failed = true;
1359	}
1360	// Rho
1361	error = std::fabs(x: initialRho - calibratedRho);
1362	if (error > calibrationTolerance) {
1363	BOOST_ERROR("\nfailed to calibrate rho Sabr parameter:"
1364	<< "\n expected: " << initialRho
1365	<< "\n calibrated: " << calibratedRho
1366	<< "\n error: " << error);
1367	failed = true;
1368	}
1369
1370	if (failed)
1371	BOOST_FAIL("\nSabr calibration failure:"
1372	<< "\n isAlphaFixed: " << k_a
1373	<< "\n isBetaFixed: " << k_b
1374	<< "\n isNuFixed: " << k_n
1375	<< "\n isRhoFixed: " << k_r
1376	<< "\n vegaWeighted[i]: " << i);
1377	}
1378	}
1379	}
1380	}
1381	}
1382	}
1383	}
1384
1385
1386	void InterpolationTest::testKernelInterpolation() {
1387
1388	BOOST_TEST_MESSAGE("Testing kernel 1D interpolation...");
1389
1390	std::vector<Real> deltaGrid = {0.10, 0.25, 0.50, 0.75, 0.90};
1391
1392	std::vector<Real> yd2(deltaGrid.size()); // test y-values 2
1393
1394	std::vector<Real> yd3(deltaGrid.size()); // test y-values 3
1395
1396	std::vector<std::vector<Real>> yd = {
1397	{11.275, 11.125, 11.250, 11.825, 12.625},
1398	{16.025, 13.450, 11.350, 10.150, 10.075},
1399	{10.300, 9.6375, 9.2000, 9.1125, 9.4000}
1400	};
1401	std::vector<Real> lambdaVec = {0.05, 0.50, 0.75, 1.65, 2.55};
1402
1403	Real tolerance = 2.0e-5;
1404
1405	Real expectedVal;
1406	Real calcVal;
1407
1408	// Check that y-values at knots are exactly the feeded y-values,
1409	// irrespective of kernel parameters
1410	for (Real i : lambdaVec) {
1411	GaussianKernel myKernel(0, i);
1412
1413	for (auto currY : yd) {
1414
1415	KernelInterpolation f(deltaGrid.begin(), deltaGrid.end(), currY.begin(), myKernel
1416	#ifdef __FAST_MATH__
1417	,
1418	1e-6
1419	#endif
1420	);
1421	f.update();
1422
1423	for (Size dIt=0; dIt< deltaGrid.size(); ++dIt) {
1424	expectedVal=currY[dIt];
1425	calcVal=f(deltaGrid[dIt]);
1426
1427	if (std::fabs(x: expectedVal-calcVal)>tolerance) {
1428
1429	BOOST_ERROR("Kernel interpolation failed at x = "
1430	<< deltaGrid[dIt]
1431	<< std::scientific
1432	<< "\n interpolated value: " << calcVal
1433	<< "\n expected value: " << expectedVal
1434	<< "\n error: "
1435	<< std::fabs(expectedVal-calcVal));
1436	}
1437	}
1438	}
1439	}
1440
1441	std::vector<Real> testDeltaGrid = {0.121, 0.279, 0.678, 0.790, 0.980};
1442
1443	// Gaussian Kernel values for testDeltaGrid with a standard
1444	// deviation of 2.05 (the value is arbitrary.) Source: parrallel
1445	// implementation in R, no literature sources found
1446
1447	std::vector<std::vector<Real>> ytd = {
1448	{11.23847, 11.12003, 11.58932, 11.99168, 13.29650},
1449	{15.55922, 13.11088, 10.41615, 10.05153, 10.50741},
1450	{10.17473, 9.557842, 9.09339, 9.149687, 9.779971}
1451	};
1452
1453	GaussianKernel myKernel(0,2.05);
1454
1455	for (Size j=0; j< ytd.size(); ++j) {
1456	std::vector<Real> currY=yd[j];
1457	std::vector<Real> currTY=ytd[j];
1458
1459	// Build interpolation according to original grid + y-values
1460	KernelInterpolation f(deltaGrid.begin(), deltaGrid.end(),
1461	currY.begin(), myKernel);
1462	f.update();
1463
1464	// test values at test Grid
1465	for (Size dIt=0; dIt< testDeltaGrid.size(); ++dIt) {
1466
1467	expectedVal=currTY[dIt];
1468	f.enableExtrapolation();// allow extrapolation
1469
1470	calcVal=f(testDeltaGrid[dIt]);
1471	if (std::fabs(x: expectedVal-calcVal)>tolerance) {
1472
1473	BOOST_ERROR("Kernel interpolation failed at x = "
1474	<< deltaGrid[dIt]
1475	<< std::scientific
1476	<< "\n interpolated value: " << calcVal
1477	<< "\n expected value: " << expectedVal
1478	<< "\n error: "
1479	<< std::fabs(expectedVal-calcVal));
1480	}
1481	}
1482	}
1483	}
1484
1485
1486	void InterpolationTest::testKernelInterpolation2D(){
1487
1488	// No test values known from the literature.
1489	// Testing for consistency of input output data
1490	// at the nodes
1491
1492	BOOST_TEST_MESSAGE("Testing kernel 2D interpolation...");
1493
1494	Real mean=0.0, var=0.18;
1495	GaussianKernel myKernel(mean,var);
1496
1497	std::vector<Real> xVec(10);
1498	xVec[0] = 0.10; xVec[1] = 0.20; xVec[2] = 0.30; xVec[3] = 0.40;
1499	xVec[4] = 0.50; xVec[5] = 0.60; xVec[6] = 0.70; xVec[7] = 0.80;
1500	xVec[8] = 0.90; xVec[9] = 1.00;
1501
1502	std::vector<Real> yVec(3);
1503	yVec[0] = 1.0; yVec[1] = 2.0; yVec[2] = 3.5;
1504
1505	Matrix M(xVec.size(),yVec.size());
1506
1507	M[0][0]=0.25; M[1][0]=0.24; M[2][0]=0.23; M[3][0]=0.20; M[4][0]=0.19;
1508	M[5][0]=0.20; M[6][0]=0.21; M[7][0]=0.22; M[8][0]=0.26; M[9][0]=0.29;
1509
1510	M[0][1]=0.27; M[1][1]=0.26; M[2][1]=0.25; M[3][1]=0.22; M[4][1]=0.21;
1511	M[5][1]=0.22; M[6][1]=0.23; M[7][1]=0.24; M[8][1]=0.28; M[9][1]=0.31;
1512
1513	M[0][2]=0.21; M[1][2]=0.22; M[2][2]=0.27; M[3][2]=0.29; M[4][2]=0.24;
1514	M[5][2]=0.28; M[6][2]=0.25; M[7][2]=0.22; M[8][2]=0.29; M[9][2]=0.30;
1515
1516	KernelInterpolation2D kernel2D(xVec.begin(),xVec.end(),
1517	yVec.begin(),yVec.end(),M,myKernel);
1518
1519	Real calcVal,expectedVal;
1520	Real tolerance = 1.0e-10;
1521
1522	for(Size i=0;i<M.rows();++i){
1523	for(Size j=0;j<M.columns();++j){
1524
1525	calcVal=kernel2D(xVec[i],yVec[j]);
1526	expectedVal=M[i][j];
1527
1528	if(std::fabs(x: expectedVal-calcVal)>tolerance){
1529
1530	BOOST_ERROR("2D Kernel interpolation failed at x = " << xVec[i]
1531	<< ", y = " << yVec[j]
1532	<< "\n interpolated value: " << calcVal
1533	<< "\n expected value: " << expectedVal
1534	<< "\n error: "
1535	<< std::fabs(expectedVal-calcVal));
1536	}
1537	}
1538	}
1539
1540	// alternative data set
1541	std::vector<Real> xVec1(4);
1542	xVec1[0] = 80.0; xVec1[1] = 90.0; xVec1[2] = 100.0; xVec1[3] = 110.0;
1543
1544	std::vector<Real> yVec1(8);
1545	yVec1[0] = 0.5; yVec1[1] = 0.7; yVec1[2] = 1.0; yVec1[3] = 2.0;
1546	yVec1[4] = 3.5; yVec1[5] = 4.5; yVec1[6] = 5.5; yVec1[7] = 6.5;
1547
1548	Matrix M1(xVec1.size(),yVec1.size());
1549	M1[0][0]=10.25; M1[1][0]=12.24;M1[2][0]=14.23;M1[3][0]=17.20;
1550	M1[0][1]=12.25; M1[1][1]=15.24;M1[2][1]=16.23;M1[3][1]=16.20;
1551	M1[0][2]=12.25; M1[1][2]=13.24;M1[2][2]=13.23;M1[3][2]=17.20;
1552	M1[0][3]=13.25; M1[1][3]=15.24;M1[2][3]=12.23;M1[3][3]=19.20;
1553	M1[0][4]=14.25; M1[1][4]=16.24;M1[2][4]=13.23;M1[3][4]=12.20;
1554	M1[0][5]=15.25; M1[1][5]=17.24;M1[2][5]=14.23;M1[3][5]=12.20;
1555	M1[0][6]=16.25; M1[1][6]=13.24;M1[2][6]=15.23;M1[3][6]=10.20;
1556	M1[0][7]=14.25; M1[1][7]=14.24;M1[2][7]=16.23;M1[3][7]=19.20;
1557
1558	// test with function pointer
1559	KernelInterpolation2D kernel2DEp(xVec1.begin(),xVec1.end(),
1560	yVec1.begin(),yVec1.end(),M1,
1561	&epanechnikovKernel);
1562
1563	for(Size i=0;i<M1.rows();++i){
1564	for(Size j=0;j<M1.columns();++j){
1565
1566	calcVal=kernel2DEp(xVec1[i],yVec1[j]);
1567	expectedVal=M1[i][j];
1568
1569	if(std::fabs(x: expectedVal-calcVal)>tolerance){
1570
1571	BOOST_ERROR("2D Epanechnkikov Kernel interpolation failed at x = " << xVec1[i]
1572	<< ", y = " << yVec1[j]
1573	<< "\n interpolated value: " << calcVal
1574	<< "\n expected value: " << expectedVal
1575	<< "\n error: "
1576	<< std::fabs(expectedVal-calcVal));
1577	}
1578	}
1579	}
1580
1581	// test updating mechanism by changing initial variables
1582	xVec1[0] = 60.0; xVec1[1] = 95.0; xVec1[2] = 105.0; xVec1[3] = 135.0;
1583
1584	yVec1[0] = 12.5; yVec1[1] = 13.7; yVec1[2] = 15.0; yVec1[3] = 19.0;
1585	yVec1[4] = 26.5; yVec1[5] = 27.5; yVec1[6] = 29.2; yVec1[7] = 36.5;
1586
1587	kernel2DEp.update();
1588
1589	for(Size i=0;i<M1.rows();++i){
1590	for(Size j=0;j<M1.columns();++j){
1591
1592	calcVal=kernel2DEp(xVec1[i],yVec1[j]);
1593	expectedVal=M1[i][j];
1594
1595	if(std::fabs(x: expectedVal-calcVal)>tolerance){
1596
1597	BOOST_ERROR("2D Epanechnkikov Kernel updated interpolation failed at x = " << xVec1[i]
1598	<< ", y = " << yVec1[j]
1599	<< "\n interpolated value: " << calcVal
1600	<< "\n expected value: " << expectedVal
1601	<< "\n error: "
1602	<< std::fabs(expectedVal-calcVal));
1603	}
1604	}
1605	}
1606	}
1607
1608
1609	void InterpolationTest::testBicubicDerivatives() {
1610	BOOST_TEST_MESSAGE("Testing bicubic spline derivatives...");
1611
1612	std::vector<Real> x(100), y(100);
1613	for (Size i=0; i < 100; ++i) {
1614	x[i] = y[i] = i/20.0;
1615	}
1616
1617	Matrix f(100, 100);
1618	for (Size i=0; i < 100; ++i)
1619	for (Size j=0; j < 100; ++j)
1620	f[i][j] = y[i]/10*std::sin(x: x[j])+std::cos(x: y[i]);
1621
1622	const Real tol=0.005;
1623	BicubicSpline spline(x.begin(), x.end(), y.begin(), y.end(), f);
1624
1625	for (Size i=5; i < 95; i+=10) {
1626	for (Size j=5; j < 95; j+=10) {
1627	Real f_x = spline.derivativeX(x: x[j],y: y[i]);
1628	Real f_xx = spline.secondDerivativeX(x: x[j],y: y[i]);
1629	Real f_y = spline.derivativeY(x: x[j],y: y[i]);
1630	Real f_yy = spline.secondDerivativeY(x: x[j],y: y[i]);
1631	Real f_xy = spline.derivativeXY(x: x[j],y: y[i]);
1632
1633	if (std::fabs(x: f_x - y[i]/10*std::cos(x: x[j])) > tol) {
1634	BOOST_ERROR("Failed to reproduce f_x");
1635	}
1636	if (std::fabs(x: f_xx + y[i]/10*std::sin(x: x[j])) > tol) {
1637	BOOST_ERROR("Failed to reproduce f_xx");
1638	}
1639	if (std::fabs(x: f_y - (std::sin(x: x[j])/10-std::sin(x: y[i]))) > tol) {
1640	BOOST_ERROR("Failed to reproduce f_y");
1641	}
1642	if (std::fabs(x: f_yy + std::cos(x: y[i])) > tol) {
1643	BOOST_ERROR("Failed to reproduce f_yy");
1644	}
1645	if (std::fabs(x: f_xy - std::cos(x: x[j])/10) > tol) {
1646	BOOST_ERROR("Failed to reproduce f_xy");
1647	}
1648	}
1649	}
1650	}
1651
1652
1653	void InterpolationTest::testBicubicUpdate() {
1654	BOOST_TEST_MESSAGE("Testing that bicubic splines actually update...");
1655
1656	Size N=6;
1657	std::vector<Real> x(N), y(N);
1658	for (Size i=0; i < N; ++i) {
1659	x[i] = y[i] = i*0.2;
1660	}
1661
1662	Matrix f(N, N);
1663	for (Size i=0; i < N; ++i)
1664	for (Size j=0; j < N; ++j)
1665	f[i][j] = x[j]*(x[j] + y[i]);
1666
1667	BicubicSpline spline(x.begin(), x.end(), y.begin(), y.end(), f);
1668
1669	Real old_result = spline(x[2]+0.1, y[4]);
1670
1671	// modify input matrix and update.
1672	f[4][3] += 1.0;
1673	spline.update();
1674
1675	Real new_result = spline(x[2]+0.1, y[4]);
1676	if (std::fabs(x: old_result-new_result) < 0.5) {
1677	BOOST_ERROR("Failed to update bicubic spline");
1678	}
1679	}
1680
1681
1682	namespace {
1683	class GF {
1684	public:
1685	GF(Real exponent, Real factor)
1686	: exponent_(exponent), factor_(factor) {}
1687
1688	Real operator()(Real h) const {
1689	return M_PI + factor_*std::pow(x: h, y: exponent_)
1690	+ std::pow(x: factor_*h, y: exponent_ + 1);
1691	}
1692	private:
1693	const Real exponent_, factor_;
1694	};
1695
1696	Real limCos(Real h) {
1697	return -std::cos(x: h);
1698	}
1699	}
1700
1701	void InterpolationTest::testUnknownRichardsonExtrapolation() {
1702	BOOST_TEST_MESSAGE("Testing Richardson extrapolation with "
1703	"unknown order of convergence...");
1704
1705	const Real stepSize = 0.01;
1706
1707	const std::pair<Real, Real> testCases[] = {
1708	std::make_pair(x: 1.0, y: 1.0), std::make_pair(x: 1.0, y: -1.0),
1709	std::make_pair(x: 2.0, y: 0.25), std::make_pair(x: 2.0, y: -1.0),
1710	std::make_pair(x: 3.0, y: 2.0), std::make_pair(x: 3.0, y: -0.5),
1711	std::make_pair(x: 4.0, y: 1.0), std::make_pair(x: 4.0, y: 0.5)
1712	};
1713
1714	for (auto testCase : testCases) {
1715	const RichardsonExtrapolation extrap(GF(testCase.first, testCase.second), stepSize);
1716
1717	const Real calculated = extrap(4.0, 2.0);
1718	const Real diff = std::fabs(M_PI - calculated);
1719
1720	const Real tol = std::pow(x: stepSize, y: testCase.first+1);
1721
1722	if (diff > tol) {
1723	BOOST_ERROR("failed to reproduce Richardson extrapolation "
1724	" with unknown order of convergence");
1725	}
1726	}
1727
1728	const Real highOrder = RichardsonExtrapolation(GF(14.0, 1.0), 0.5)(4.0,2.0);
1729	if (std::fabs(x: highOrder - M_PI) > 1e-12) {
1730	BOOST_ERROR("failed to reproduce Richardson extrapolation "
1731	" with unknown order of convergence");
1732	}
1733
1734	try {
1735	RichardsonExtrapolation(GF(16.0, 1.0), 0.5)(4.0,2.0);
1736	BOOST_ERROR("Richardson extrapolation with order of"
1737	" convergence above 15 should throw exception");
1738	}
1739	catch (...) {}
1740
1741	const Real limCosValue = RichardsonExtrapolation(limCos, 0.01)(4.0,2.0);
1742	if (std::fabs(x: limCosValue + 1.0) > 1e-6)
1743	BOOST_ERROR("failed to reproduce Richardson extrapolation "
1744	" with unknown order of convergence");
1745	}
1746
1747
1748	namespace {
1749	Real f(Real h) {
1750	return std::pow( x: 1.0 + h, y: 1/h);
1751	}
1752	}
1753
1754	void InterpolationTest::testRichardsonExtrapolation() {
1755	BOOST_TEST_MESSAGE("Testing Richardson extrapolation...");
1756
1757	/* example taken from
1758	* http://www.ipvs.uni-stuttgart.de/abteilungen/bv/lehre/
1759	* lehrveranstaltungen/vorlesungen/WS0910/
1760	* NSG_termine/dateien/Richardson.pdf
1761	*/
1762
1763	const Real stepSize = 0.1;
1764	const Real orderOfConvergence = 1.0;
1765	const RichardsonExtrapolation extrap(f, stepSize, orderOfConvergence);
1766
1767
1768	Real tol = 0.00002;
1769	Real expected = 2.71285;
1770
1771	const Real scalingFactor = 2.0;
1772	Real calculated = extrap(scalingFactor);
1773
1774	if (std::fabs(x: expected-calculated) > tol) {
1775	BOOST_ERROR("failed to reproduce Richardson extrapolation");
1776	}
1777
1778	calculated = extrap();
1779	if (std::fabs(x: expected-calculated) > tol) {
1780	BOOST_ERROR("failed to reproduce Richardson extrapolation");
1781	}
1782
1783	expected = 2.721376;
1784	const Real scalingFactor2 = 4.0;
1785	calculated = extrap(scalingFactor2, scalingFactor);
1786
1787	if (std::fabs(x: expected-calculated) > tol) {
1788	BOOST_ERROR("failed to reproduce Richardson extrapolation");
1789	}
1790	}
1791
1792	void InterpolationTest::testNoArbSabrInterpolation(){
1793
1794	BOOST_TEST_MESSAGE("Testing no-arbitrage Sabr interpolation...");
1795
1796	// Test SABR function against input volatilities
1797	#ifndef __FAST_MATH__
1798	Real tolerance = 1.0e-12;
1799	#else
1800	Real tolerance = 1.0e-8;
1801	#endif
1802	std::vector<Real> strikes(31);
1803	std::vector<Real> volatilities(31), volatilities2(31);
1804	// input strikes
1805	strikes[0] = 0.03 ; strikes[1] = 0.032 ; strikes[2] = 0.034 ;
1806	strikes[3] = 0.036 ; strikes[4] = 0.038 ; strikes[5] = 0.04 ;
1807	strikes[6] = 0.042 ; strikes[7] = 0.044 ; strikes[8] = 0.046 ;
1808	strikes[9] = 0.048 ; strikes[10] = 0.05 ; strikes[11] = 0.052 ;
1809	strikes[12] = 0.054 ; strikes[13] = 0.056 ; strikes[14] = 0.058 ;
1810	strikes[15] = 0.06 ; strikes[16] = 0.062 ; strikes[17] = 0.064 ;
1811	strikes[18] = 0.066 ; strikes[19] = 0.068 ; strikes[20] = 0.07 ;
1812	strikes[21] = 0.072 ; strikes[22] = 0.074 ; strikes[23] = 0.076 ;
1813	strikes[24] = 0.078 ; strikes[25] = 0.08 ; strikes[26] = 0.082 ;
1814	strikes[27] = 0.084 ; strikes[28] = 0.086 ; strikes[29] = 0.088;
1815	strikes[30] = 0.09;
1816	// input volatilities for noarb sabr (other than above
1817	// alpha is 0.2 here due to the restriction sigmaI <= 1.0 !)
1818	volatilities[0] = 0.773729077752926;
1819	volatilities[1] = 0.763916242454194;
1820	volatilities[2] = 0.754773878663612;
1821	volatilities[3] = 0.746222305031368;
1822	volatilities[4] = 0.738193023523582;
1823	volatilities[5] = 0.730629785825930;
1824	volatilities[6] = 0.723484825471685;
1825	volatilities[7] = 0.716716812668892;
1826	volatilities[8] = 0.710290301049393;
1827	volatilities[9] = 0.704174528906769;
1828	volatilities[10] = 0.698342635400901;
1829	volatilities[11] = 0.692771033345972;
1830	volatilities[12] = 0.687438902593476;
1831	volatilities[13] = 0.682327777297265;
1832	volatilities[14] = 0.677421206991904;
1833	volatilities[15] = 0.672704476238547;
1834	volatilities[16] = 0.668164371832768;
1835	volatilities[17] = 0.663788984329375;
1836	volatilities[18] = 0.659567547226380;
1837	volatilities[19] = 0.655490294349232;
1838	volatilities[20] = 0.651548341349061;
1839	volatilities[21] = 0.647733583657137;
1840	volatilities[22] = 0.644038608699086;
1841	volatilities[23] = 0.640456620061898;
1842	volatilities[24] = 0.636981371712714;
1843	volatilities[25] = 0.633607110719560;
1844	volatilities[26] = 0.630328527192861;
1845	volatilities[27] = 0.627140710386248;
1846	volatilities[28] = 0.624039110072250;
1847	volatilities[29] = 0.621019502453590;
1848	volatilities[30] = 0.618077959983455;
1849
1850	Time expiry = 1.0;
1851	Real forward = 0.039;
1852	// input SABR coefficients (corresponding to the vols above)
1853	Real initialAlpha = 0.2;
1854	Real initialBeta = 0.6;
1855	Real initialNu = 0.02;
1856	Real initialRho = 0.01;
1857	// calculate SABR vols and compare with input vols
1858	NoArbSabrSmileSection noarbSabr(expiry, forward,
1859	{initialAlpha, initialBeta, initialNu, initialRho});
1860	for (Size i = 0; i < strikes.size(); i++) {
1861	Real calculatedVol = noarbSabr.volatility(strike: strikes[i]);
1862	if (std::fabs(x: volatilities[i]-calculatedVol) > tolerance)
1863	BOOST_ERROR(
1864	"failed to calculate noarb-Sabr function at strike " << strikes[i]
1865	<< "\n expected: " << volatilities[i]
1866	<< "\n calculated: " << calculatedVol
1867	<< "\n error: " << std::fabs(calculatedVol-volatilities[i]));
1868	}
1869
1870	// Test SABR calibration against input parameters
1871	Real betaGuess = 0.5;
1872	Real alphaGuess = 0.2 / std::pow(x: forward,y: betaGuess-1.0); // new default value for alpha
1873	Real nuGuess = std::sqrt(x: 0.4);
1874	Real rhoGuess = 0.0;
1875
1876	const bool vegaWeighted[]= {true, false};
1877	const bool isAlphaFixed[]= {true, false};
1878	const bool isBetaFixed[]= {true, false};
1879	const bool isNuFixed[]= {true, false};
1880	const bool isRhoFixed[]= {true, false};
1881
1882	Real calibrationTolerance = 5.0e-6;
1883	// initialize optimization methods
1884	std::vector<ext::shared_ptr<OptimizationMethod>> methods_ = {
1885	ext::shared_ptr<OptimizationMethod>(new Simplex(0.01)),
1886	ext::shared_ptr<OptimizationMethod>(new LevenbergMarquardt(1e-8, 1e-8, 1e-8))
1887	};
1888	// Initialize end criteria
1889	ext::shared_ptr<EndCriteria> endCriteria(new
1890	EndCriteria(100000, 100, 1e-8, 1e-8, 1e-8));
1891	// Test looping over all possibilities
1892	for (Size j=1; j<methods_.size(); ++j) { // skip simplex (gets caught in some cases)
1893	for (bool i : vegaWeighted) {
1894	for (bool k_a : isAlphaFixed) {
1895	for (Size k_b=0; k_b<1/*LENGTH(isBetaFixed)*/; ++k_b) { // keep beta fixed (all 4 params free is a problem for this kind of test)
1896	for (bool k_n : isNuFixed) {
1897	for (bool k_r : isRhoFixed) {
1898	NoArbSabrInterpolation noarbSabrInterpolation(
1899	strikes.begin(), strikes.end(), volatilities.begin(), expiry,
1900	forward, k_a ? initialAlpha : alphaGuess,
1901	isBetaFixed[k_b] ? initialBeta : betaGuess,
1902	k_n ? initialNu : nuGuess, k_r ? initialRho : rhoGuess, k_a,
1903	isBetaFixed[k_b], k_n, k_r, i, endCriteria, methods_[j], 1E-10);
1904	noarbSabrInterpolation.update();
1905
1906	// Recover SABR calibration parameters
1907	bool failed = false;
1908	Real calibratedAlpha = noarbSabrInterpolation.alpha();
1909	Real calibratedBeta = noarbSabrInterpolation.beta();
1910	Real calibratedNu = noarbSabrInterpolation.nu();
1911	Real calibratedRho = noarbSabrInterpolation.rho();
1912	Real error;
1913
1914	// compare results: alpha
1915	error = std::fabs(x: initialAlpha-calibratedAlpha);
1916	if (error > calibrationTolerance) {
1917	BOOST_ERROR("\nfailed to calibrate alpha Sabr parameter:" <<
1918	"\n expected: " << initialAlpha <<
1919	"\n calibrated: " << calibratedAlpha <<
1920	"\n error: " << error);
1921	failed = true;
1922	}
1923	// Beta
1924	error = std::fabs(x: initialBeta-calibratedBeta);
1925	if (error > calibrationTolerance) {
1926	BOOST_ERROR("\nfailed to calibrate beta Sabr parameter:" <<
1927	"\n expected: " << initialBeta <<
1928	"\n calibrated: " << calibratedBeta <<
1929	"\n error: " << error);
1930	failed = true;
1931	}
1932	// Nu
1933	error = std::fabs(x: initialNu-calibratedNu);
1934	if (error > calibrationTolerance) {
1935	BOOST_ERROR("\nfailed to calibrate nu Sabr parameter:" <<
1936	"\n expected: " << initialNu <<
1937	"\n calibrated: " << calibratedNu <<
1938	"\n error: " << error);
1939	failed = true;
1940	}
1941	// Rho
1942	error = std::fabs(x: initialRho-calibratedRho);
1943	if (error > calibrationTolerance) {
1944	BOOST_ERROR("\nfailed to calibrate rho Sabr parameter:" <<
1945	"\n expected: " << initialRho <<
1946	"\n calibrated: " << calibratedRho <<
1947	"\n error: " << error);
1948	failed = true;
1949	}
1950
1951	if (failed)
1952	BOOST_TEST_MESSAGE("\nnoarb-Sabr calibration failure:"
1953	<< "\n isAlphaFixed: " << k_a
1954	<< "\n isBetaFixed: " << isBetaFixed[k_b]
1955	<< "\n isNuFixed: " << k_n
1956	<< "\n isRhoFixed: " << k_r
1957	<< "\n vegaWeighted[i]: " << i);
1958	}
1959	}
1960	}
1961	}
1962	}
1963	}
1964
1965	}
1966
1967
1968	void InterpolationTest::testSabrSingleCases() {
1969
1970	BOOST_TEST_MESSAGE("Testing Sabr calibration single cases...");
1971
1972	// case #1
1973	// this fails with an exception thrown in 1.4, fixed in 1.5
1974
1975	std::vector<Real> strikes = { 0.01, 0.01125, 0.0125, 0.01375, 0.0150 };
1976	std::vector<Real> vols = { 0.1667, 0.2020, 0.2785, 0.3279, 0.3727 };
1977
1978	Real tte = 0.3833;
1979	Real forward = 0.011025;
1980
1981	SABRInterpolation s0(strikes.begin(), strikes.end(), vols.begin(), tte, forward,
1982	Null<Real>(), 0.25, Null<Real>(), Null<Real>(),
1983	false, true, false, false);
1984	s0.update();
1985
1986	if (s0.maxError() > 0.01 || s0.rmsError() > 0.01) {
1987	BOOST_ERROR("Sabr case #1 failed with max error ("
1988	<< s0.maxError() << ") and rms error (" << s0.rmsError()
1989	<< "), both should be < 0.01");
1990	}
1991
1992	}
1993
1994	void InterpolationTest::testTransformations() {
1995
1996	BOOST_TEST_MESSAGE("Testing Sabr and no-arbitrage Sabr transformation functions...");
1997
1998	Real size = 25.0; // test inputs from [-size,size]^4
1999
2000	Size N = 100000;
2001
2002	Array x(4), y(4), z(4);
2003	std::vector<Real> s;
2004	std::vector<bool> fixed(4, false);
2005	std::vector<Real> params(4, 0.0);
2006	Real forward = 0.03;
2007
2008	HaltonRsg h(4, 42, false, false);
2009
2010	for (Size i = 0; i < 1E6; ++i) {
2011
2012	s = h.nextSequence().value;
2013	for (Size j = 0; j < 4; ++j)
2014	x[j] = 2.0 * size * s[j] - size;
2015
2016	// sabr
2017	y = QuantLib::detail::SABRSpecs().direct(x, fixed, params, forward);
2018	validateSabrParameters(alpha: y[0], beta: y[1], nu: y[2], rho: y[3]);
2019	z = QuantLib::detail::SABRSpecs().inverse(y, fixed, params, forward);
2020	z = QuantLib::detail::SABRSpecs().direct(x: z, fixed, params, forward);
2021	if (!close(x: z[0], y: y[0], n: N) || !close(x: z[1], y: y[1], n: N) || !close(x: z[2], y: y[2], n: N) ||
2022	!close(x: z[3], y: y[3], n: N))
2023	BOOST_ERROR("SabrInterpolation: direct(inverse("
2024	<< y[0] << "," << y[1] << "," << y[2] << "," << y[3]
2025	<< ")) = (" << z[0] << "," << z[1] << "," << z[2] << ","
2026	<< z[3] << "), difference is (" << z[0] - y[0] << ","
2027	<< z[1] - y[1] << "," << z[2] - y[2] << ","
2028	<< z[3] - y[3] << ")");
2029
2030	// noarb sabr
2031	y = QuantLib::detail::NoArbSabrSpecs().direct(x, paramIsFixed: fixed, params, forward);
2032
2033	// we can not invoke the constructor, this would be too slow, so
2034	// we copy the parameter check here ...
2035	Real alpha = y[0];
2036	Real beta = y[1];
2037	Real nu = y[2];
2038	Real rho = y[3];
2039	QL_REQUIRE(beta >= QuantLib::detail::NoArbSabrModel::beta_min &&
2040	beta <= QuantLib::detail::NoArbSabrModel::beta_max,
2041	"beta (" << beta << ") out of bounds");
2042	Real sigmaI = alpha * std::pow(x: forward, y: beta - 1.0);
2043	QL_REQUIRE(sigmaI >= QuantLib::detail::NoArbSabrModel::sigmaI_min &&
2044	sigmaI <= QuantLib::detail::NoArbSabrModel::sigmaI_max,
2045	"sigmaI = alpha*forward^(beta-1.0) ("
2046	<< sigmaI << ") out of bounds, alpha=" << alpha
2047	<< " beta=" << beta << " forward=" << forward);
2048	QL_REQUIRE(nu >= QuantLib::detail::NoArbSabrModel::nu_min &&
2049	nu <= QuantLib::detail::NoArbSabrModel::nu_max,
2050	"nu (" << nu << ") out of bounds");
2051	QL_REQUIRE(rho >= QuantLib::detail::NoArbSabrModel::rho_min &&
2052	rho <= QuantLib::detail::NoArbSabrModel::rho_max,
2053	"rho (" << rho << ") out of bounds");
2054
2055	z = QuantLib::detail::NoArbSabrSpecs().inverse(y, paramIsFixed: fixed, params, forward);
2056	z = QuantLib::detail::NoArbSabrSpecs().direct(x: z, paramIsFixed: fixed, params, forward);
2057	if (!close(x: z[0], y: y[0], n: N) || !close(x: z[1], y: y[1], n: N) || !close(x: z[2], y: y[2], n: N) ||
2058	!close(x: z[3], y: y[3], n: N))
2059	BOOST_ERROR("NoArbSabrInterpolation: direct(inverse("
2060	<< y[0] << "," << y[1] << "," << y[2] << "," << y[3]
2061	<< ")) = (" << z[0] << "," << z[1] << "," << z[2] << ","
2062	<< z[3] << "), difference is (" << z[0] - y[0] << ","
2063	<< z[1] - y[1] << "," << z[2] - y[2] << ","
2064	<< z[3] - y[3] << ")");
2065	}
2066	}
2067
2068	void InterpolationTest::testFlochKennedySabrIsSmoothAroundATM() {
2069	BOOST_TEST_MESSAGE("Testing that Andersen SABR formula is smooth "
2070	"close to the ATM level...");
2071
2072	const Real f0 = 1.1;
2073	const Real alpha = 0.35;
2074	const Real nu = 1.1;
2075	const Real rho = 0.25;
2076	const Real beta = 0.3;
2077	const Real strike= f0;
2078	const Time t = 2.1;
2079
2080	const Real vol = sabrFlochKennedyVolatility(strike, forward: f0, expiryTime: t, alpha, beta, nu, rho);
2081
2082	const Real expected = 0.3963883944;
2083	const Real tol = 1e-8;
2084	const Real diff = std::fabs(x: expected - vol);
2085	if (diff > tol) {
2086	BOOST_ERROR("\nfailed to get ATM value :" <<
2087	"\n expected: " << expected <<
2088	"\n calculated: " << vol <<
2089	"\n diff: " << diff);
2090	}
2091
2092	Real k = 0.996*strike;
2093	Real v = sabrFlochKennedyVolatility(strike: k, forward: f0, expiryTime: t, alpha, beta, nu, rho);
2094
2095	for (; k < 1.004*strike; k += 0.0001*strike) {
2096	const Real vt = sabrFlochKennedyVolatility(strike: k, forward: f0, expiryTime: t, alpha, beta, nu, rho);
2097
2098	const Real diff = std::fabs(x: v - vt);
2099
2100	if (diff > 1e-5) {
2101	BOOST_ERROR("\nSabr vol spike around ATM :" <<
2102	"\n volatility at " << k-0.0001*strike <<
2103	" is " << v <<
2104	"\n volatility at " << k << " is " << vt <<
2105	"\n difference: " << diff <<
2106	"\n tolerance : " << 1e-5);
2107	}
2108	v = vt;
2109	}
2110	}
2111
2112	void InterpolationTest::testLeFlochKennedySabrExample() {
2113	BOOST_TEST_MESSAGE("Testing Le Floc'h Kennedy SABR Example...");
2114
2115	/*
2116	Example is taken from F. Le Floc'h, G. Kennedy:
2117	Explicit SABR Calibration through Simple Expansions.
2118	https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2467231
2119	*/
2120
2121	const Real f0 = 1.0;
2122	const Real alpha = 0.35;
2123	const Real nu = 1.0;
2124	const Real rho = 0.25;
2125	const Real beta = 0.25;
2126	const Real strikes[]= {1.0, 1.5, 0.5};
2127	const Time t = 2.0;
2128
2129	const Real expected[] = {0.408702473958, 0.428489933046, 0.585701651161};
2130
2131	for (Size i=0; i < LENGTH(strikes); ++i) {
2132	const Real strike = strikes[i];
2133	const Real vol =
2134	sabrFlochKennedyVolatility(strike, forward: f0, expiryTime: t, alpha, beta, nu, rho);
2135
2136	const Real tol = 1e-8;
2137	const Real diff = std::fabs(x: expected[i] - vol);
2138
2139	if (diff > tol) {
2140	BOOST_ERROR("\nfailed to reproduce reference examples :" <<
2141	"\n expected: " << expected[i] <<
2142	"\n calculated: " << vol <<
2143	"\n diff: " << diff);
2144	}
2145	}
2146	}
2147
2148	namespace {
2149	Real lagrangeTestFct(Real x) {
2150	return std::fabs(x: x) + 0.5*x - x*x;
2151	}
2152	}
2153
2154	void InterpolationTest::testLagrangeInterpolation() {
2155
2156	BOOST_TEST_MESSAGE("Testing Lagrange interpolation...");
2157
2158	const std::vector<Real> x = {-1.0 , -0.5, -0.25, 0.1, 0.4, 0.75, 0.96};
2159	Array y(x.size());
2160	std::transform(first: x.begin(), last: x.end(), result: y.begin(), unary_op: &lagrangeTestFct);
2161
2162	LagrangeInterpolation interpl(x.begin(), x.end(), y.begin());
2163
2164	// reference results are taken from R package pracma
2165	const Real references[] = {
2166	-0.5000000000000000,-0.5392414024347419,-0.5591485962711904,
2167	-0.5629199661387594,-0.5534414777017116,-0.5333043347921566,
2168	-0.5048221831582063,-0.4700478608272949,-0.4307896950846587,
2169	-0.3886273460669714,-0.3449271969711449,-0.3008572908782903,
2170	-0.2574018141928359,-0.2153751266968088,-0.1754353382192734,
2171	-0.1380974319209344,-0.1037459341938971,-0.0726471311765894,
2172	-0.0449608318838433,-0.0207516779521373,0.0000000000000000,
2173	0.0173877793964286,0.0315691961126723,0.0427562482700356,
2174	0.0512063534145595,0.0572137590808174,0.0611014067405497,
2175	0.0632132491361394,0.0639070209989264,0.0635474631523613,
2176	0.0625000000000000,0.0611248703983366,0.0597717119144768,
2177	0.0587745984686508,0.0584475313615655,0.0590803836865967,
2178	0.0609352981268212,0.0642435381368876,0.0692027925097279,
2179	0.0759749333281079,0.0846842273010179,0.0954160004849021,
2180	0.1082157563897290,0.1230887474699003,0.1400000000000001,
2181	0.1588747923353829,0.1795995865576031,0.2020234135046815,
2182	0.2259597111862140,0.2511886165833182,0.2774597108334206,
2183	0.3044952177998833,0.3319936560264689,0.3596339440766487,
2184	0.3870799592577457,0.4139855497299214,0.4400000000000001,
2185	0.4647739498001331,0.4879657663513030,0.5092483700116673,
2186	0.5283165133097421,0.5448945133624253,0.5587444376778583,
2187	0.5696747433431296,0.5775493695968156,0.5822972837863635,
2188	0.5839224807103117,0.5825144353453510,0.5782590089582251,
2189	0.5714498086024714,0.5625000000000000,0.5519545738075141,
2190	0.5405030652677689,0.5289927272456703,0.5184421566492137,
2191	0.5100553742352614,0.5052363578001620,0.5056040287552059,
2192	0.5130076920869246
2193	};
2194
2195	constexpr double tol = 50*QL_EPSILON;
2196	for (Size i=0; i < 79; ++i) {
2197	const Real xx = -1.0 + i*0.025;
2198	const Real calculated = interpl(xx);
2199	if ( std::isnan(x: calculated)
2200	|| std::fabs(x: references[i] - calculated) > tol) {
2201	BOOST_FAIL("failed to reproduce the Lagrange interpolation"
2202	<< "\n x : " << xx
2203	<< "\n calculated: " << calculated
2204	<< "\n expected : " << references[i]);
2205	}
2206	}
2207	}
2208
2209	void InterpolationTest::testLagrangeInterpolationAtSupportPoint() {
2210	BOOST_TEST_MESSAGE(
2211	"Testing Lagrange interpolation at supporting points...");
2212
2213	const Size n=5;
2214	Array x(n), y(n);
2215	for (Size i=0; i < n; ++i) {
2216	x[i] = i/Real(n);
2217	y[i] = 1.0/(1.0 - x[i]);
2218	}
2219	LagrangeInterpolation interpl(x.begin(), x.end(), y.begin());
2220
2221	const Real relTol = 5e-12;
2222
2223	for (Size i=1; i < n-1; ++i) {
2224	for (Real z = x[i] - 100*QL_EPSILON;
2225	z < x[i] + 100*QL_EPSILON; z+=2*QL_EPSILON) {
2226	const Real expected = 1.0/(1.0 - x[i]);
2227	const Real calculated = interpl(z);
2228
2229	if ( std::isnan(x: calculated)
2230	|| std::fabs(x: expected - calculated) > relTol) {
2231	BOOST_FAIL("failed to reproduce the Lagrange interplation"
2232	<< "\n x : " << z
2233	<< "\n calculated: " << calculated
2234	<< "\n expected : " << expected);
2235	}
2236	}
2237	}
2238	}
2239
2240	void InterpolationTest::testLagrangeInterpolationDerivative() {
2241	BOOST_TEST_MESSAGE(
2242	"Testing Lagrange interpolation derivatives...");
2243
2244	Array x(5), y(5);
2245	x[0] = -1.0; y[0] = 2.0;
2246	x[1] = -0.3; y[1] = 3.0;
2247	x[2] = 0.1; y[2] = 6.0;
2248	x[3] = 0.3; y[3] = 3.0;
2249	x[4] = 0.9; y[4] =-1.0;
2250
2251	LagrangeInterpolation interpl(x.begin(), x.end(), y.begin());
2252
2253	const Real eps = std::sqrt(QL_EPSILON);
2254	for (Real x=-1.0; x <= 0.9; x+=0.01) {
2255	const Real calculated = interpl.derivative(x, allowExtrapolation: true);
2256	const Real expected = (interpl(x+eps, true)
2257	- interpl(x-eps, true))/(2*eps);
2258
2259	if ( std::isnan(x: calculated)
2260	|| std::fabs(x: expected - calculated) > 25*eps) {
2261	BOOST_FAIL("failed to reproduce the Lagrange"
2262	" interplation derivative"
2263	<< "\n x : " << x
2264	<< "\n calculated: " << calculated
2265	<< "\n expected : " << expected);
2266	}
2267	}
2268	}
2269
2270	void InterpolationTest::testLagrangeInterpolationOnChebyshevPoints() {
2271	BOOST_TEST_MESSAGE(
2272	"Testing Lagrange interpolation on Chebyshev nodes...");
2273
2274	// Test example taken from
2275	// J.P. Berrut, L.N. Trefethen, Barycentric Lagrange Interpolation
2276	// https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
2277
2278	const Size n=50;
2279	Array x(n+1), y(n+1);
2280	for (Size i=0; i <= n; ++i) {
2281	// Chebyshev nodes
2282	x[i] = std::cos( x: (2*i+1)*M_PI/(2*n+2) );
2283	y[i] = std::exp(x: x[i])/std::cos(x: x[i]);
2284	}
2285
2286	LagrangeInterpolation interpl(x.begin(), x.end(), y.begin());
2287
2288	const Real tol = 1e-13;
2289	const Real tolDeriv = 1e-11;
2290
2291	for (Real x=-1.0; x <= 1.0; x+=0.03) {
2292	const Real calculated = interpl(x, true);
2293	const Real expected = std::exp(x: x)/std::cos(x: x);
2294
2295	const Real diff = std::fabs(x: expected - calculated);
2296	if (std::isnan(x: calculated) || diff > tol) {
2297	BOOST_FAIL("failed to reproduce the Lagrange"
2298	" interpolation on Chebyshev nodes"
2299	<< "\n x : " << x
2300	<< "\n calculated: " << calculated
2301	<< "\n expected : " << expected
2302	<< std::scientific
2303	<< "\n difference: " << diff);
2304	}
2305
2306	const Real calculatedDeriv = interpl.derivative(x, allowExtrapolation: true);
2307	const Real expectedDeriv = std::exp(x: x)*(std::cos(x: x) + std::sin(x: x)) / squared(x: std::cos(x: x));
2308
2309	const Real diffDeriv = std::fabs(x: expectedDeriv - calculatedDeriv);
2310	if (std::isnan(x: calculated) || diffDeriv > tolDeriv) {
2311	BOOST_FAIL("failed to reproduce the Lagrange"
2312	" interpolation derivative on Chebyshev nodes"
2313	<< "\n x : " << x
2314	<< "\n calculated: " << calculatedDeriv
2315	<< "\n expected : " << expectedDeriv
2316	<< std::scientific
2317	<< "\n difference: " << diffDeriv);
2318	}
2319	}
2320	}
2321
2322	void InterpolationTest::testBSplines() {
2323	BOOST_TEST_MESSAGE("Testing B-Splines...");
2324
2325	// reference values have been generate with the R package splines2
2326	// https://cran.r-project.org/web/packages/splines2/splines2.pdf
2327
2328	std::vector<Real> knots = { -1.0, 0.5, 0.75, 1.2, 3.0, 4.0, 5.0 };
2329
2330	const Natural p = 2;
2331	const BSpline bspline(p, knots.size()-p-2, knots);
2332
2333	std::vector<ext::tuple<Natural, Real, Real>> referenceValues = {
2334	ext::make_tuple(args: 0, args: -0.95, args: 9.5238095238e-04),
2335	ext::make_tuple(args: 0, args: -0.01, args: 0.37337142857),
2336	ext::make_tuple(args: 0, args: 0.49, args: 0.84575238095),
2337	ext::make_tuple(args: 0, args: 1.21, args: 0.0),
2338	ext::make_tuple(args: 1, args: 1.49, args: 0.562987654321),
2339	ext::make_tuple(args: 1, args: 1.59, args: 0.490888888889),
2340	ext::make_tuple(args: 2, args: 1.99, args: 0.62429409171),
2341	ext::make_tuple(args: 3, args: 1.19, args: 0.0),
2342	ext::make_tuple(args: 3, args: 1.99, args: 0.12382936508),
2343	ext::make_tuple(args: 3, args: 3.59, args: 0.765914285714)
2344	};
2345
2346	const Real tol = 1e-10;
2347	for (auto& referenceValue : referenceValues) {
2348	const Natural idx = ext::get<0>(t&: referenceValue);
2349	const Real x = ext::get<1>(t&: referenceValue);
2350	const Real expected = ext::get<2>(t&: referenceValue);
2351
2352	const Real calculated = bspline(idx, x);
2353
2354	if ( std::isnan(x: calculated)
2355	|| std::fabs(x: calculated - expected) > tol) {
2356	BOOST_FAIL("failed to reproduce the B-Spline value"
2357	<< "\n i : " << idx
2358	<< "\n x : " << x
2359	<< "\n calculated: " << calculated
2360	<< "\n expected : " << expected
2361	<< "\n difference: " << std::fabs(calculated-expected)
2362	<< "\n tolerance : " << tol);
2363	}
2364	}
2365	}
2366
2367	void InterpolationTest::testBackwardFlatOnSinglePoint() {
2368	BOOST_TEST_MESSAGE("Testing piecewise constant interpolation on a "
2369	"single point...");
2370	const std::vector<Real> knots(1, 1.0), values(1, 2.5);
2371
2372	const Interpolation impl(BackwardFlat().interpolate(
2373	xBegin: knots.begin(), xEnd: knots.end(), yBegin: values.begin()));
2374
2375	const Real x[] = { -1.0, 1.0, 2.0, 3.0 };
2376
2377	for (Real i : x) {
2378	const Real calculated = impl(i, true);
2379	const Real expected = values[0];
2380
2381	if (!close_enough(x: calculated, y: expected)) {
2382	BOOST_FAIL("failed to reproduce a piecewise constant "
2383	"interpolation on a single point "
2384	<< "\n x : " << i << "\n expected : " << expected
2385	<< "\n calculated: " << calculated);
2386	}
2387
2388	const Real expectedPrimitive = values[0] * (i - knots[0]);
2389	const Real calculatedPrimitive = impl.primitive(x: i, allowExtrapolation: true);
2390
2391	if (!close_enough(x: calculatedPrimitive, y: expectedPrimitive)) {
2392	BOOST_FAIL("failed to reproduce primitive on a piecewise constant "
2393	"interpolation for a single point "
2394	<< "\n x : " << i << "\n expected : " << expectedPrimitive
2395	<< "\n calculated: " << calculatedPrimitive);
2396	}
2397	}
2398	}
2399
2400	void InterpolationTest::testChebyshevInterpolation() {
2401	BOOST_TEST_MESSAGE("Testing Chebyshev interpolation...");
2402
2403	const auto fcts =
2404	std::vector<std::pair<std::function<Real(Real)>, std::string> >{
2405	{[](Real x) { return std::sin(x: x); }, "sin"},
2406	{[](Real x) { return std::cos(x: x); }, "cos"},
2407	{[](Real x) { return std::exp(x: -x*x); }, "e^(-x*x)"}
2408	};
2409
2410	const auto tests = std::vector<std::pair<Size, Real> >{
2411	{11, 1e-5},
2412	{20, 1e-11}
2413	};
2414
2415	for (const auto& t: tests) {
2416	for (const auto& fct: fcts) {
2417	ChebyshevInterpolation interp(t.first, fct.first);
2418
2419	for (Real x=-0.99; x < 1.0; x+=0.01) {
2420	const Real expected = fct.first(x);
2421	const Real calculated = interp(x);
2422	const Real diff = std::fabs(x: expected-calculated);
2423	const Real tol = t.second;
2424
2425	if ( std::isnan(x: calculated)
2426	|| std::fabs(x: calculated - expected) > tol) {
2427	BOOST_FAIL("failed to reproduce the Chebyshev interpolation values"
2428	<< "\n x : " << x
2429	<< "\n fct : " << fct.second
2430	<< "\n calculated: " << calculated
2431	<< "\n expected : " << expected
2432	<< "\n difference: " << diff
2433	<< "\n tolerance : " << tol);
2434	}
2435	}
2436	}
2437	}
2438	}
2439
2440	void InterpolationTest::testChebyshevInterpolationOnNodes() {
2441	BOOST_TEST_MESSAGE("Testing Chebyshev interpolation on and around nodes...");
2442
2443	constexpr double tol = 10*QL_EPSILON;
2444	const auto testFct = [](Real x) { return std::sin(x: x);};
2445
2446	const Size nrNodes = 7;
2447	Array y(nrNodes);
2448
2449	for (auto pointType: {ChebyshevInterpolation::FirstKind,
2450	ChebyshevInterpolation::SecondKind}) {
2451
2452	const Array nodes = ChebyshevInterpolation::nodes(n: nrNodes, pointsType: pointType);
2453	std::transform(first: std::begin(cont: nodes), last: std::end(cont: nodes), result: std::begin(cont&: y), unary_op: testFct);
2454
2455	const ChebyshevInterpolation interp(y, pointType);
2456	for (auto node: nodes) {
2457	// test on Chebyshev node
2458	const Real expected = testFct(node);
2459	const Real calculated = interp(node);
2460	const Real diff = std::abs(x: expected - calculated);
2461
2462	if (diff > tol) {
2463	BOOST_ERROR("failed to reproduce the node values"
2464	<< std::setprecision(16)
2465	<< "\n node : " << node
2466	<< "\n calculated: " << calculated
2467	<< "\n expected : " << expected
2468	<< "\n difference: " << diff
2469	<< "\n tolerance : " << tol);
2470	}
2471
2472
2473	// check around Chebyshev node
2474	for (Integer i=-50; i < 50; ++i) {
2475	const Real x = node + i*QL_EPSILON;
2476	const Real expected = testFct(x);
2477	const Real calculated = interp(x, true);
2478	const Real diff = std::abs(x: expected - calculated);
2479
2480	if (diff > tol) {
2481	BOOST_ERROR("failed to reproduce values around nodes"
2482	<< std::setprecision(16)
2483	<< "\n node : " << node
2484	<< "\n epsilon : " << x - node
2485	<< "\n calculated: " << calculated
2486	<< "\n expected : " << expected
2487	<< "\n difference: " << diff
2488	<< "\n tolerance : " << tol);
2489	}
2490	}
2491	}
2492	}
2493	}
2494
2495	void InterpolationTest::testChebyshevInterpolationUpdateY() {
2496	BOOST_TEST_MESSAGE("Testing Y update for Chebyshev interpolation...");
2497
2498	Array y({1, 4, 7, 4});
2499	ChebyshevInterpolation interp(y);
2500
2501	Array yd({6, 4, 5, 6});
2502	interp.updateY(y: yd);
2503
2504	constexpr double tol = 10*QL_EPSILON;
2505
2506	for (Size i=0; i < y.size(); ++i) {
2507	const Real expected = yd[i];
2508	const Real calculated = interp(interp.nodes()[i], true);
2509	const Real diff = std::abs(x: calculated - expected);
2510
2511	if (diff > tol) {
2512	BOOST_ERROR("failed to reproduce updated node values"
2513	<< std::setprecision(16)
2514	<< "\n node : " << i
2515	<< "\n expected : " << expected
2516	<< "\n calculated: " << calculated
2517	<< "\n difference: " << diff
2518	<< "\n tolerance : " << tol);
2519	}
2520	}
2521	}
2522
2523
2524	test_suite* InterpolationTest::suite(SpeedLevel speed) {
2525	auto* suite = BOOST_TEST_SUITE("Interpolation tests");
2526
2527	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testSplineOnGenericValues));
2528	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testSimmetricEndConditions));
2529	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testDerivativeEndConditions));
2530	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testNonRestrictiveHymanFilter));
2531	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testSplineOnRPN15AValues));
2532	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testSplineOnGaussianValues));
2533	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testSplineErrorOnGaussianValues));
2534	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testMultiSpline));
2535	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testAsFunctor));
2536	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testFritschButland));
2537	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testBackwardFlat));
2538	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testForwardFlat));
2539	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testSabrInterpolation));
2540	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testFlochKennedySabrIsSmoothAroundATM));
2541	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testLeFlochKennedySabrExample));
2542	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testKernelInterpolation));
2543	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testKernelInterpolation2D));
2544	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testBicubicDerivatives));
2545	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testBicubicUpdate));
2546	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testUnknownRichardsonExtrapolation));
2547	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testRichardsonExtrapolation));
2548	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testSabrSingleCases));
2549	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testTransformations));
2550	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testLagrangeInterpolation));
2551	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testLagrangeInterpolationAtSupportPoint));
2552	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testLagrangeInterpolationDerivative));
2553	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testLagrangeInterpolationOnChebyshevPoints));
2554	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testBSplines));
2555	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testBackwardFlatOnSinglePoint));
2556	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testChebyshevInterpolation));
2557	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testChebyshevInterpolationOnNodes));
2558	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testChebyshevInterpolationUpdateY));
2559	if (speed <= Fast) {
2560	suite->add(QUANTLIB_TEST_CASE(&InterpolationTest::testNoArbSabrInterpolation));
2561	}
2562
2563	return suite;
2564	}
2565