Same as my answer.Just read this thread today and got

Volatility = 0.04

Opt_Price=9.01002030921698e-25

Implied_Volatilty =0.0400000003783122

using a biscetion with a smart initial bounds guess.

It converged using a tol 1e-30 after 24 steps

- Cuchulainn
**Posts:**19117**Joined:****Location:**89 19 79 15

Same as my answer.Just read this thread today and got

Volatility = 0.04

Opt_Price=9.01002030921698e-25

Implied_Volatilty =0.0400000003783122

using a biscetion with a smart initial bounds guess.

It converged using a tol 1e-30 after 24 steps

What's infeasible about it? Arbitrary large K is fine theoretically and has a strictly positive call value.With the same data above, SHGO algo fails to find a minimizer point when K > 290. Suggests infeasible region.

https://shgo.readthedocs.io/en/stable/docs/README.html

- Cuchulainn
**Posts:**19117**Joined:****Location:**89 19 79 15

I just started with SGHO as a black box. Here is the runnable Python code.What's infeasible about it? Arbitrary large K is fine theoretically and has a strictly positive call value.With the same data above, SHGO algo fails to find a minimizer point when K > 290. Suggests infeasible region.

https://shgo.readthedocs.io/en/stable/docs/README.html

What's strange is that iv == 0.02 for K > 290, no idea why, maybe because bounds = [0.02, 0.4] and it gets stuck at the lower limit?

Code: Select all

```
import numpy as np
from scipy import optimize
from scipy.stats import norm
# Option data, vol to be computed
# Stress test case
S = 100.0
K = 320 # K > 290 not OK
r = 0.0
b = r
x = 0.04 #sig/vol
T = 1
q = 0
# Cumulative normal N(0,1)
def N(x):
return norm.cdf(x)
# Black Scholes formula for a plain call option
def BS(x):
d1 = (np.log(S/K) + (r - q + x**2/2)*T) / (x*np.sqrt(T))
d2 = d1 - x* np.sqrt(T)
return S*np.exp(-q*T) * N(d1) - K * np.exp(-r*T)* N(d2)
bounds = [(0.02,0.4)]
marketValue = BS(x); print ("option price ", marketValue)
f = lambda x: (BS(x) - marketValue)*(BS(x) - marketValue)
def MySGHO(method, nSamples, nIter):
result = optimize.shgo(f, bounds, n = nSamples, iters = nIter, sampling_method =method)
print("roots ", result)
print("global minimum ", result.x)
print("function value: ", result.fun + marketValue)
# Testing
print('*******simplicial********')
MySGHO('simplicial', 100, 10)
print('*******sobol********')
#MySGHO('sobol', 128, 3)
print('*******halton********')
#MySGHO('halton', 128, 3)
```

- Cuchulainn
**Posts:**19117**Joined:****Location:**89 19 79 15

Forget SHGO.

It's very good, but was not built for iv.

It's very good, but was not built for iv.

Last edited by Cuchulainn on June 11th, 2023, 8:26 pm, edited 1 time in total.

- Cuchulainn
**Posts:**19117**Joined:****Location:**89 19 79 15

Now, this article by Dan Stefanica and Rados Radoicic of Baruch College, Lexington at 24th Street, NYU approximate N(x) by Polya's 1949 A(x) . I typed the code in C++ and and passes with flying colours till now.

https://papers.ssrn.com/sol3/papers.cfm ... id=2908494

const double S = 100;

const double K = 15000;

const double r = 0.00;

const double b = r;

const double sig = 0.004;

const double T = 10.0;

https://papers.ssrn.com/sol3/papers.cfm ... id=2908494

const double S = 100;

const double K = 15000;

const double r = 0.00;

const double b = r;

const double sig = 0.004;

const double T = 10.0;

I hope that my recent paper is helpful for the discussion: Choi et al. (2023), https://arxiv.org/abs/2302.08758 .

The paper finds quite many new IV bounds (both upper and lower) for a given price (this is main theme of the paper).

Then, it also formulates a new Newton-Raphon method on the log price to handle very small option premiums properly. The iteration formula is quite simple (see Eq. (30) in the paper), and it always converges as long as the initial guess is a lower bound. (The proof is a bit tedious. Just see the left side of Figure 1. The log price is a concave function of sigma.)

So we use a lower bound found in the earlier part of the paper (specifically L3 in Eq (23)).

I actually tested the famous example in this thread (S0=1, K=1.5, sigma=4%) and reported it in the paper (See Section 3.3 Numerical Example(. The new NR method reaches an IV value within 2e-11 just after three iterations.

The paper finds quite many new IV bounds (both upper and lower) for a given price (this is main theme of the paper).

Then, it also formulates a new Newton-Raphon method on the log price to handle very small option premiums properly. The iteration formula is quite simple (see Eq. (30) in the paper), and it always converges as long as the initial guess is a lower bound. (The proof is a bit tedious. Just see the left side of Figure 1. The log price is a concave function of sigma.)

So we use a lower bound found in the earlier part of the paper (specifically L3 in Eq (23)).

I actually tested the famous example in this thread (S0=1, K=1.5, sigma=4%) and reported it in the paper (See Section 3.3 Numerical Example(. The new NR method reaches an IV value within 2e-11 just after three iterations.

You joined this forum before Cuchulainn (with 18,934 posts) and you have bided your time until now to say something. I’m guessing that may be a record. Respect! Oh, and your paper looks nice!

- jherekhealy
**Posts:**20**Joined:**

Nice paper **jaesmine,** I also had proposed the idea of solving on the log price to speed up calculations a while ago in a blog post (and based on the nice guess of Dan Stefanica and Rados Radoicic).

https://quantsrus.github.io/post/implied_volatility_algorithms_go_julia/

https://quantsrus.github.io/post/implied_volatility_algorithms_go_julia/

I joined a while ago (guess I am old enough). I often checked this forum but did not post much. Recently, this thread got my attention as I was working on that paper.You joined this forum before Cuchulainn (with 18,934 posts) and you have bided your time until now to say something. I’m guessing that may be a record. Respect! Oh, and your paper looks nice!

Nice paperjaesmine,I also had proposed the idea of solving on the log price to speed up calculations a while ago in a blog post (and based on the nice guess of Dan Stefanica and Rados Radoicic).

https://quantsrus.github.io/post/implied_volatility_algorithms_go_julia/

- Cuchulainn
**Posts:**19117**Joined:****Location:**89 19 79 15

Here is the C++ "get it working" code I wrote based on the Dan Stefanica and Rados Radoicic paper.

As I mentioned, I tried many exotic and no so exotic methods but this one seems to the best.

As I mentioned, I tried many exotic and no so exotic methods but this one seems to the best.

Code: Select all

```
// DS_RR.cpp
//
/* An Explicit Implied Volatility Formula
International Journal of Theoreticaland Applied Finance, Vol. 20, no. 7, 2017
Dan Stefanica
Baruch College, City University of New York
Rados Radoicic
CUNY Baruch College
Date Written : January 30, 2017
*/
//
// https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2908494
//
// Programmed in C++ by Daniel J.Duffy 2023
// Datasim Education BV 2023
//
#include <cmath>
#include <iostream>
using value_type = double;
// Call parameters
/*
const value_type S = 2080;
const value_type K = 3050;
const value_type r = 0.02;
const value_type b = 0;
const value_type sig = 0.1097904;
const value_type T = 0.5;
*/
const value_type S = 100;
const value_type K = 15000;
const value_type r = 0.08;
const value_type b = r;
const value_type sig = 0.0004;
const value_type T = 1000.0;
// ATM + small vol OK
// OTM, ITM + small vol not OK
// r == vol == 0 is essentially a degenerate case??
/*
const value_type S = 60;
const value_type K = 65;
const value_type r = 0.08;
const value_type b = r;
const value_type sig = 0.302; // small vol has issues < 0.1 !!
const value_type T = 0.25;
*/
/*
const value_type S = 100;
const value_type K = 6.30957;
const value_type r = 0.0;
const value_type b = r;
const value_type sig = 0.3409436;
const value_type T = 1.0;
*/
// Polya's 1949 approximation to N(x)
value_type G(value_type x)
{
return 1 / (exp(-358 * x / 23 + 111 * atan(37 * x / 294) + 1));
value_type y = std::sqrt(1.0 - std::exp(-2.0 * x * x / 3.141592653589793));
if (x >= 0.0)
return 0.5 * (1.0 + y);
else
return 0.5 * (1.0 - y);
}
value_type N(value_type x)
{ // The approximation to the cumulative normal distribution
//return cndN(x);
return G(x);
}
value_type OptionPrice(value_type x)
{ // Call price, x == sig
value_type tmp = x * std::sqrt(T);
value_type d1 = (std::log(S / K) + (r + (x * x) * 0.5) * T) / tmp;
value_type d2 = d1 - tmp;
return (S * std::exp((b - r) * T) * N(d1)) - (K * std::exp(-r * T) * N(d2));
}
value_type Cm = OptionPrice(sig); // Market price
// Stefanica and Radoicic volatiliy estimator
value_type y = std::log(S / K);
value_type alpha = Cm / (K * std::exp(-r * T));
value_type R = 2 * alpha - std::exp(y) + 1.0;
value_type pi = 3.141592653589793;
value_type t1 = std::exp((1.0 - 2 / pi) * y);
value_type t2 = std::exp((-1.0 + 2 / pi) * y);
value_type A = std::pow(t1 - t2, 2);
value_type B1 = 4.0 * (std::exp(-y) + std::exp(-2 * y / pi));
value_type B2 = -2.0 * (std::exp(-y) * (t1 + t2) * (std::exp(2 * y) + 1.0 - R * R));
value_type B = B1 + B2;
value_type C = std::exp(-2.0 * y)* (R * R - std::pow(std::exp(y) + 1.0, 2) - R * R);
value_type beta = 2.0 * C / (B + std::sqrt(B * B + 4.0 * A*C));
value_type gamma = -0.5 * pi * std::log(beta);
value_type compute()
{
if (y >= 0.0)
{
value_type C0 = K * std::exp(-r * T) * (std::exp(y) * A * (std::sqrt(2 * y) - 0.5));
value_type sig;
if (Cm <= C0)
{
sig = (std::sqrt(gamma + y) - std::sqrt(gamma - y)) / std::sqrt(T);
}
else
{
sig = (std::sqrt(gamma + y) + std::sqrt(gamma - y)) / std::sqrt(T);
}
}
else
{
value_type C0 = K * std::exp(-r * T) * (std::exp(y)/2 - A * (std::sqrt(2 * y) - 0.5));
value_type sig;
if (Cm <= C0)
{
sig = (-std::sqrt(gamma + y) + std::sqrt(gamma - y)) / std::sqrt(T);
}
else
{
sig = (std::sqrt(gamma + y) - std::sqrt(gamma - y)) / std::sqrt(T);
}
}
return sig;
}
int main()
{
std::cout << "DS/RR: " << compute() << '\n';
}
```

- Cuchulainn
**Posts:**19117**Joined:****Location:**89 19 79 15

This looks like NR with some upfront work to find an initial seed (which is always the main challenge/bottleneck).I hope that my recent paper is helpful for the discussion: Choi et al. (2023), https://arxiv.org/abs/2302.08758 .

The paper finds quite many new IV bounds (both upper and lower) for a given price (this is main theme of the paper).

Then, it also formulates a new Newton-Raphon method on the log price to handle very small option premiums properly. The iteration formula is quite simple (see Eq. (30) in the paper), and it always converges as long as the initial guess is a lower bound. (The proof is a bit tedious. Just see the left side of Figure 1. The log price is a concave function of sigma.)

So we use a lower bound found in the earlier part of the paper (specifically L3 in Eq (23)).

I actually tested the famous example in this thread (S0=1, K=1.5, sigma=4%) and reported it in the paper (See Section 3.3 Numerical Example(. The new NR method reaches an IV value within 2e-11 just after three iterations.

Mention is made of iterating a number of times to get a desired accuracy which may affect performance.

If it takes more than 2-3 iterations for a simple BS model IV, then one is not bringing to bear everything that is already known about how to get there faster. Speaking of infeasible regions, it is always useful to first check using certain boundary conditions whether an IV calculation is even possible, I think Hull's texts used to point that out. Very, very good bang for the buck comes from finding a good seed calculation. I use a modified Brent method as per Numerical Recipes.

- Cuchulainn
**Posts:**19117**Joined:****Location:**89 19 79 15

Yeah, kind of begging the question "let's pretend we are near the solution and then we can fnd it". Kind of depressing. Spent time fixing silly problems MBS by BDT with initial yield 5% (because from a book....)

Improvements to NR are homotopy and continuation methods which converge with 'bad' seeds.

And fixed point always converges.

Improvements to NR are homotopy and continuation methods which converge with 'bad' seeds.

And fixed point always converges.

- Cuchulainn
**Posts:**19117**Joined:****Location:**89 19 79 15