blackscholes 0.24.0

Black-Scholes option pricing model calculator
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

//
// This source code resides at www.jaeckel.org/LetsBeRational.7z .
//
// ======================================================================================
// Copyright © 2013-2014 Peter Jäckel.
//
// Permission to use, copy, modify, and distribute this software is freely granted,
// provided that this notice is preserved.
//
// WARRANTY DISCLAIMER
// The Software is provided "as is" without warranty of any kind, either express or implied,
// including without limitation any implied warranties of condition, uninterrupted use,
// merchantability, fitness for a particular purpose, or non-infringement.
// ======================================================================================
//

#include "rationalcubic.h"

#if defined(_MSC_VER)
#define NOMINMAX // to suppress MSVC's definitions of min() and max()
// These four pragmas are the equivalent to /fp:fast.
// YOU NEED THESE FOR THE SAKE OF *ACCURACY* WHEN |x| IS LARGE, say, |x|>50.
// This is because they effectively enable the evaluation of certain
// expressions in 80 bit registers without loss of intermediate accuracy.
#pragma float_control(except, off)
#pragma float_control(precise, off)
#pragma fp_contract(on)
#pragma fenv_access(off)
#endif

#include <algorithm>
#include <cmath>
#include <float.h>

// Based on
//
//    “Shape preserving piecewise rational interpolation”, R. Delbourgo, J.A. Gregory - SIAM journal on scientific and statistical computing, 1985 - SIAM.
//    http://dspace.brunel.ac.uk/bitstream/2438/2200/1/TR_10_83.pdf  [caveat emptor: there are some typographical errors in that draft version]
//

namespace {
const double minimum_rational_cubic_control_parameter_value = -(1 - sqrt(DBL_EPSILON));
const double maximum_rational_cubic_control_parameter_value = 2 / (DBL_EPSILON * DBL_EPSILON);
inline bool is_zero(double x) {
    return fabs(x) < DBL_MIN;
}
} // namespace

double rational_cubic_interpolation(double x, double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double r) {
    const double h = (x_r - x_l);
    if (fabs(h) <= 0)
        return 0.5 * (y_l + y_r);
    // r should be greater than -1. We do not use  assert(r > -1)  here in order to allow values such as NaN to be propagated as they should.
    const double t = (x - x_l) / h;
    if (!(r >= maximum_rational_cubic_control_parameter_value)) {
        const double t = (x - x_l) / h, omt = 1 - t, t2 = t * t, omt2 = omt * omt;
        // Formula (2.4) divided by formula (2.5)
        return (y_r * t2 * t + (r * y_r - h * d_r) * t2 * omt + (r * y_l + h * d_l) * t * omt2 + y_l * omt2 * omt) / (1 + (r - 3) * t * omt);
    }
    // Linear interpolation without over-or underflow.
    return y_r * t + y_l * (1 - t);
}

double rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double second_derivative_l) {
    const double h = (x_r - x_l), numerator = 0.5 * h * second_derivative_l + (d_r - d_l);
    if (is_zero(numerator))
        return 0;
    const double denominator = (y_r - y_l) / h - d_l;
    if (is_zero(denominator))
        return numerator > 0 ? maximum_rational_cubic_control_parameter_value : minimum_rational_cubic_control_parameter_value;
    return numerator / denominator;
}

double rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double second_derivative_r) {
    const double h = (x_r - x_l), numerator = 0.5 * h * second_derivative_r + (d_r - d_l);
    if (is_zero(numerator))
        return 0;
    const double denominator = d_r - (y_r - y_l) / h;
    if (is_zero(denominator))
        return numerator > 0 ? maximum_rational_cubic_control_parameter_value : minimum_rational_cubic_control_parameter_value;
    return numerator / denominator;
}

double minimum_rational_cubic_control_parameter(double d_l, double d_r, double s, bool preferShapePreservationOverSmoothness) {
    const bool monotonic = d_l * s >= 0 && d_r * s >= 0, convex = d_l <= s && s <= d_r, concave = d_l >= s && s >= d_r;
    if (!monotonic && !convex && !concave) // If 3==r_non_shape_preserving_target, this means revert to standard cubic.
        return minimum_rational_cubic_control_parameter_value;
    const double d_r_m_d_l = d_r - d_l, d_r_m_s = d_r - s, s_m_d_l = s - d_l;
    double r1 = -DBL_MAX, r2 = r1;
    // If monotonicity on this interval is possible, set r1 to satisfy the monotonicity condition (3.8).
    if (monotonic) {
        if (!is_zero(s))                                         // (3.8), avoiding division by zero.
            r1 = (d_r + d_l) / s;                                // (3.8)
        else if (preferShapePreservationOverSmoothness)          // If division by zero would occur, and shape preservation is preferred, set value to enforce linear interpolation.
            r1 = maximum_rational_cubic_control_parameter_value; // This value enforces linear interpolation.
    }
    if (convex || concave) {
        if (!(is_zero(s_m_d_l) || is_zero(d_r_m_s))) // (3.18), avoiding division by zero.
            r2 = std::max(fabs(d_r_m_d_l / d_r_m_s), fabs(d_r_m_d_l / s_m_d_l));
        else if (preferShapePreservationOverSmoothness)
            r2 = maximum_rational_cubic_control_parameter_value; // This value enforces linear interpolation.
    } else if (monotonic && preferShapePreservationOverSmoothness)
        r2 = maximum_rational_cubic_control_parameter_value; // This enforces linear interpolation along segments that are inconsistent with the slopes on the boundaries, e.g., a perfectly horizontal segment that has negative slopes on either edge.
    return std::max(minimum_rational_cubic_control_parameter_value, std::max(r1, r2));
}

double convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double second_derivative_l, bool preferShapePreservationOverSmoothness) {
    const double r = rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(x_l, x_r, y_l, y_r, d_l, d_r, second_derivative_l);
    const double r_min = minimum_rational_cubic_control_parameter(d_l, d_r, (y_r - y_l) / (x_r - x_l), preferShapePreservationOverSmoothness);
    return std::max(r, r_min);
}

double convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(double x_l, double x_r, double y_l, double y_r, double d_l, double d_r, double second_derivative_r, bool preferShapePreservationOverSmoothness) {
    const double r = rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(x_l, x_r, y_l, y_r, d_l, d_r, second_derivative_r);
    const double r_min = minimum_rational_cubic_control_parameter(d_l, d_r, (y_r - y_l) / (x_r - x_l), preferShapePreservationOverSmoothness);
    return std::max(r, r_min);
}