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)

// 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';
}