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
|
# File 'lib/pure_greeks/implied_volatility/brent_solver.rb', line 12
def find_root(lower:, upper:, tolerance: 1e-8, &f)
a = lower.to_f
b = upper.to_f
fa = f.call(a)
fb = f.call(b)
raise IVConvergenceError, "root not bracketed: f(#{a})=#{fa}, f(#{b})=#{fb}" if (fa * fb).positive?
if fa.abs < fb.abs
a, b = b, a
fa, fb = fb, fa
end
c = a
fc = fa
mflag = true
d = nil
MAX_ITERATIONS.times do
return b if fb.abs < tolerance || (b - a).abs < tolerance
s =
if fa != fc && fb != fc
(a * fb * fc / ((fa - fb) * (fa - fc))) +
(b * fa * fc / ((fb - fa) * (fb - fc))) +
(c * fa * fb / ((fc - fa) * (fc - fb)))
else
b - (fb * (b - a) / (fb - fa))
end
condition1 = !s.between?([(3 * a + b) / 4, b].min, [(3 * a + b) / 4, b].max)
condition2 = mflag && (s - b).abs >= (b - c).abs / 2
condition3 = !mflag && (s - b).abs >= (c - d).abs / 2
condition4 = mflag && (b - c).abs < tolerance
condition5 = !mflag && d && (c - d).abs < tolerance
if condition1 || condition2 || condition3 || condition4 || condition5
s = (a + b) / 2.0
mflag = true
else
mflag = false
end
fs = f.call(s)
d = c
c = b
fc = fb
if (fa * fs).negative?
b = s
fb = fs
else
a = s
fa = fs
end
if fa.abs < fb.abs
a, b = b, a
fa, fb = fb, fa
end
end
raise IVConvergenceError, "exceeded #{MAX_ITERATIONS} iterations"
end
|