An Explicit Implied Volatility Formula
International Journal of Theoretical and Applied Finance, Vol. 20, no. 7, 2017
Dan Stefanica
Baruch College, City University of New York
Rados Radoicic
CUNY Baruch College
[CODE=cpp]// 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';
}[/CODE]
