FunctionBlackScholes
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © RicardoSantos

//@version=5

// @description Some methods for the Black Scholes Options Model, which demonstrates several approaches to the valuation of a European call.
library("FunctionBlackScholes")

// reference:
//      https://people.math.sc.edu/Burkardt/py_src/black_scholes/black_scholes.html
//      https://people.math.sc.edu/Burkardt/py_src/black_scholes/black_scholes.py

import RicardoSantos/Probability/1 as prob
import RicardoSantos/ArrayExtension/2 as ae
import RicardoSantos/ArrayGenerate/1 as ag
import RicardoSantos/DebugConsole/6 as console
[__T, __C] = console.init(25)



// @function Simulates the behavior of an asset price over time.
// @param s0 float, asset price at time 0.
// @param mu float, growth rate.
// @param sigma float, volatility.
// @param t1 float, time to expiry date.
// @param n int, time steps to expiry date.
// @returns option values at each equal timed step (0 -> t1)
export asset_path (float s0, float mu, float sigma, float t1, int n) => //{
    float _dt = t1 / n
    float[] _r = ag.normal_distribution(n, 0.0, 1.0)
    // float[] _arg = array.new_float(n)
    // float _a = math.pow(mu - sigma, 2.0) * _dt + sigma * math.sqrt(_dt)
    // float _cumprod = 1.0
    // float[] _s = array.new_float(n)
    // for _i = 0 to (n - 1)
    //     float _ri = array.get(_r, _i)
    //     float _ai = math.exp(_a * _ri)
    //     _cumprod *= _ai
    //     array.set(_arg, _i, _ai)
    //     array.set(_s, _i, s0 * _cumprod)
    // // insert the initial price on the start:
    // array.unshift(_s, s0)
    float[] _s = array.new_float(n+1, s0)
    float _p = s0
    for _i = 1 to n
        _p *= math.exp((mu - math.pow(sigma, 2.0)) * _dt + sigma * math.sqrt(_dt) * array.get(_r, _i-1))
        array.set(_s, _i, _p)
    _p
//{ usage
    _s
s0 = 2.0
mu = 0.1
sigma = 0.3
n = 100
t1 = 1.0
path = asset_path(s0, mu, sigma, t1, n)
str = '---<> asset_path() <>---'
str += str.format('\n  The asset price at time 0,     S0    = {0}', s0)
str += str.format('\n  The asset expected growth rate MU    = {0}', mu)
str += str.format('\n  The asset volatility           SIGMA = {0}', sigma)
str += str.format('\n  The expiry date                T1    = {0}', t1)
str += str.format('\n  The number of time steps       N     = {0}', n)
str += str.format('\n  {0}\n', str.tostring(path, '#.000'))
str += str.format('\n  {0}\n', array.get(path, 100))
console.queue_one(__C, str)
//}}

// @function Uses the binomial method for a European call.
// @param s0 float, asset price at time 0.
// @param e float, exercise price.
// @param r float, interest rate.
// @param sigma float, volatility.
// @param t1 float, time to expiry date.
// @param m int, time steps to expiry date.
// @returns option value at time 0.
export binomial (float s0, float e, float r, float sigma, float t1, int m) => //{
    float _dt = t1 / float(m)
    float _a = 0.5 * (math.exp(-r * _dt) + math.exp((r + math.pow(sigma, 2.0)) * _dt))
    float _sqrt = math.sqrt(math.pow(_a, 2.0) - 1.0)
    float _d = (_a - _sqrt)
    float _u = (_a + _sqrt)
    float _p = (math.exp(r * _dt) - _d) / (_u - _d)
    float[] _w = array.new_float(m + 1)
    for _i = 0 to m
        array.set(_w, _i, math.max(0.0, s0 * math.pow(_d, m - _i) * math.pow(_u, _i) - e))
    // Trace backwards to get the option value at time 0.
    for _n = m-2 to -1
        for _i = 0 to (_n + 1)
            float _wi = array.get(_w, _i)
            float _wi1 = array.get(_w, _i + 1)
            array.set(_w, _i, (1.0 - _p) * _wi + _p * _wi1)
    for _i = 0 to m-1
        float _wi = array.get(_w, _i)
        array.set(_w, _i, math.exp(- r * t1) * _wi)
    array.get(_w, 0)
//{
b_s0 = 2.0
b_e = 1.0
b_r = 0.05
b_sigma = 0.25
b_t1 = 3.0
b_m = 256
b_opt = binomial(b_s0, b_e, b_r, b_sigma, b_t1, b_m)
b_str = '---<> binomial() <>---'
b_str += str.format('\n  asset price at time 0 (S0) = {0}', b_s0)
b_str += str.format('\n  exercise price (E) = {0}', b_e)
b_str += str.format('\n  interest rate (R) = {0}', b_r)
b_str += str.format('\n  asset volatility (SIGMA) = {0}', b_sigma)
b_str += str.format('\n  expiry date (T1) = {0}', b_t1)
b_str += str.format('\n  number of time steps (M) = {0}', b_m)
b_str += str.format('\n  {0}\n', b_opt)
console.queue_one(__C, b_str)
//}

// @function Evaluates the Black-Scholes formula for a European call.
// @param s0 float, asset price at time 0.
// @param t0 float, time at which the price is known.
// @param e float, exercise price.
// @param r float, interest rate.
// @param sigma float, volatility.
// @param t1 float, time to expiry date.
// @returns option value at time 0.
export bsf (float s0, float t0, float e, float r, float sigma, float t1) => //{
    float _tau = t1 - t0
    if 0.0 < _tau
        float _sqrtau = math.sqrt(_tau)
        float _d1 = (math.log(s0 / e) + (r + 0.5 * math.pow(sigma, 2.0)) * _tau) / (sigma * _sqrtau)
        float _d2 = _d1 - sigma * _sqrtau
        float _sqrt2 = math.sqrt(2.0)
        float _n1 = 0.5 * (1.0 + prob.erf(_d1 / _sqrt2))
        float _n2 = 0.5 * (1.0 + prob.erf(_d2 / _sqrt2))
        s0 * _n1 - e * math.exp(-r * _tau) * _n2
    else
        math.max(0.0, s0 - e)
//{
bsf_s0 = 2.0
bsf_t0 = 0.0
bsf_e = 1.0
bsf_r = 0.05
bsf_sigma = 0.25
bsf_t1 = 3.0
bsf_opt = bsf(bsf_s0, bsf_t0, bsf_e, bsf_r, bsf_sigma, bsf_t1)
bsf_str = '---<> bsf() <>---'
bsf_str += str.format('\n  asset price at time 0 (S0) = {0}', bsf_s0)
bsf_str += str.format('\n  last known date (S0) = {0}', bsf_t0)
bsf_str += str.format('\n  exercise price (E) = {0}', bsf_e)
bsf_str += str.format('\n  interest rate (R) = {0}', bsf_r)
bsf_str += str.format('\n  asset volatility (SIGMA) = {0}', bsf_sigma)
bsf_str += str.format('\n  expiry date (T1) = {0}', bsf_t1)
bsf_str += str.format('\n  {0}\n', bsf_opt)
console.queue_one(__C, bsf_str)
//}

// @function Forward difference method to value a European call option.
// @param e float, exercise price.
// @param r float, interest rate.
// @param sigma float, volatility.
// @param t1 float, time to expiry date.
// @param nx int, number of space steps in interval (0, L).
// @param nt int, number of time steps.
// @param smax float, maximum value of S to consider.
// @returns option values for the european call, float array of size ((nx-1) * (nt+1)).
export forward (float e, float r, float sigma, float t1, int nx, int nt, float smax) => //{
    float _dt = t1 / float(nt)
    float _dx = smax / float(nx)
    float[] _a = array.new_float(nx)
    float[] _b = array.new_float(nx)
    float[] _c = array.new_float(nx)
    for _i = 0 to (nx - 1)
        array.set(_b, _i, 1.0 - r * _dt - _dt * math.pow(sigma * (_i + 1), 2.0))
    for _i = 0 to (nx - 2)
        array.set(_c, _i, 0.5 * _dt * math.pow(sigma * (_i + 1), 2.0) + 0.5 * _dt * r * (_i + 1.0))
    for _i = 1 to (nx - 1)
        array.set(_a, _i, 0.5 * _dt * math.pow(sigma * (_i + 1), 2.0) - 0.5 * _dt * r * (_i + 1))
    float[] _u = array.new_float((nx-1)*(nt+1))
    float _u0 = 0.0
    for _i = 0 to (nx-2)
        _u0 += _dx
        array.set(_u, ae.index_2d_to_1d(nx-1, nt+1, _i, 0), math.max(0.0, _u0 - e))
    for _j = 0 to (nt-1)
        float _t = float(_j) * t1 / float(nt)
        float _p = 0.5 * _dt * (nx - 1) * (math.pow(sigma, 2.0) * (nx - 1) + r) * (smax - e * math.exp(-r * _t))
        for _i = 0 to (nx-2)
            float _bi = array.get(_b, _i)
            float _uij = array.get(_u, ae.index_2d_to_1d(nx-1, nt+1, _i, _j))
            array.set(_u, ae.index_2d_to_1d(nx-1, nt+1, _i, _j + 1), _bi * _uij)
        for _i = 0 to (nx-3)
            float _ci = array.get(_c, _i)
            float _uij1 = array.get(_u, ae.index_2d_to_1d(nx-1, nt+1, _i, _j+1))
            float _ui1j = array.get(_u, ae.index_2d_to_1d(nx-1, nt+1, _i+1, _j))
            array.set(_u, ae.index_2d_to_1d(nx-1, nt+1, _i, _j + 1), _uij1 + _ci * _ui1j)
        for _i = 1 to (nx-2)
            float _ai = array.get(_a, _i)
            float _uij1 = array.get(_u, ae.index_2d_to_1d(nx-1, nt+1, _i, _j+1))
            float _u1ij = array.get(_u, ae.index_2d_to_1d(nx-1, nt+1, _i-1, _j))
            array.set(_u, ae.index_2d_to_1d(nx-1, nt+1, _i, _j + 1), _uij1 + _ai * _u1ij)
        float _unx2 = array.get(_u, ae.index_2d_to_1d(nx-1, nt+1, nx - 2, _j+1))
        array.set(_u, ae.index_2d_to_1d(nx-1, nt+1, nx - 2, _j+1), _unx2 + _p)
    _u
//{
for_e = 4.0
for_r = 0.03
for_sigma = 0.5
for_t1 = 1.0
for_nx = 11
for_nt = 29
for_smax = 10.0
for_opt = forward(for_e, for_r, for_sigma, for_t1, for_nx, for_nt, for_smax)
for_str = '---<> forward() <>---'
for_str += str.format('\n  exercise price (E) = {0}', for_e)
for_str += str.format('\n  interest rate (R) = {0}', for_r)
for_str += str.format('\n  asset volatility (SIGMA) = {0}', for_sigma)
for_str += str.format('\n  expiry date (T1) = {0}', for_t1)
for_str += str.format('\n  number of space steps (nx) = {0}', for_nx)
for_str += str.format('\n  number of time steps (nt) = {0}', for_nt)
for_str += str.format('\n  max value to consider (smax) = {0}', for_smax)
// for_str += str.format('\n  {0}\n', str.tostring(for_opt, '#.###'))
for_str += '\n [ '
for _i = 0 to for_nx-2
    for_str += str.tostring(array.get(for_opt, ae.index_2d_to_1d(for_nx-1, for_nt+1, _i, for_nt)), '#.000') + (_i != for_nx-2 ? ', ' : ' ]\n')
// for_str += str.format('\n  {0}\n', array.get(for_opt, array.size(for_opt)-1))
console.queue_one(__C, for_str)
//}

// @function Uses Monte Carlo valuation on a European call.
// @param s0 float, asset price at time 0.
// @param e float, exercise price.
// @param r float, interest rate.
// @param sigma float, volatility.
// @param t1 float, time to expiry date.
// @param m int, time steps to expiry date.
// @returns confidence interval for the estimated range of valuation.
export mc (float s0, float e, float r, float sigma, float t1, int m) => //{
    float[] _u = ag.normal_distribution(m, 0.0, 1.0)
    float[] _pvals = array.new_float(m, 0.0)
    for _i = 0 to (m - 1)
        float _ui = array.get(_u, _i)
        float _si = s0 * math.exp((r - 0.5 * math.pow(sigma, 2.0)) * t1 + sigma * math.sqrt(t1) * _ui)
        array.set(_pvals, _i, math.exp(-r * t1) * math.max(0.0, _si - e))
    float _pmean = array.avg(_pvals)
    float _std = array.stdev(_pvals)
    float _width = 1.96 * _std / math.sqrt(float(m))
    array.from(_pmean - _width, _pmean + _width) // confidence interval
//{
mc_s0 = 2.0
mc_e = 1.0
mc_r = 0.05
mc_sigma = 0.25
mc_t1 = 3.0
mc_m = 100
mc_opt = mc(mc_s0, mc_e, mc_r, mc_sigma, mc_t1, mc_m)
mc_str = '---<> mc() <>---'
mc_str += str.format('\n  asset price at time 0 (S0) = {0}', mc_s0)
mc_str += str.format('\n  exercise price (E) = {0}', mc_e)
mc_str += str.format('\n  interest rate (R) = {0}', mc_r)
mc_str += str.format('\n  asset volatility (SIGMA) = {0}', mc_sigma)
mc_str += str.format('\n  expiry date (T1) = {0}', mc_t1)
mc_str += str.format('\n  number of time steps (M) = {0}', mc_m)
mc_str += str.format('\n  {0}\n', str.tostring(mc_opt, '#.000'))
console.queue_one(__C, mc_str)
//}



console.update(__T, __C)