Implied Volatility Estimator using Black Scholes [Loxx]
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © loxx

//@version=5
indicator("Implied Volatility Estimator using Black Scholes [Loxx]", 
     shorttitle ="IVEBS",
     overlay = true, 
     max_lines_count = 500)

if not timeframe.isdaily
    runtime.error("Error: Invald timeframe. Indicator only works on daily timeframe.") 

//imports
import loxx/loxxexpandedsourcetypes/4

string timtoolbar= "Time Now = Current time in UNIX format. It is the number of milliseconds that have elapsed since 00:00:00 UTC, 1 January 1970."
string timtoolnow = "Time Bar = The time function returns the UNIX time of the current bar for the specified timeframe and session or NaN if the time point is out of session."
string timetooltrade = "Trading Day = The beginning time of the trading day the current bar belongs to, in UNIX format (the number of milliseconds that have elapsed since 00:00:00 UTC, 1 January 1970)."

color darkGreenColor = #1B7E02 

string rogersatch = "Roger-Satchell"
string c2c = "Close-to-Close"
string parkinson = "Parkinson"
string gkvol = "Garman-Klass"
string gkzhvol = "Garman-Klass-Yang-Zhang"
string ewmavolstr = "Exponential Weighted Moving Average"

f_tickFormat() =>
    _s = str.tostring(syminfo.mintick)
    _s := str.replace_all(_s, '25', '00')
    _s := str.replace_all(_s, '5', '0')
    _s := str.replace_all(_s, '1', '0')
    _s

// N(0,1) density
f(float x)=> 
    math.exp(-x * x * 0.5) / math.sqrt(2 * math.pi)

// Boole's Rule
Boole(float StartPoint, float EndPoint, int n)=>
    float[] X = array.new<float>(n + 1 , 0)
    float[] Y = array.new<float>(n + 1 , 0)
    float delta_x = (EndPoint - StartPoint) / n
    for i = 0 to n
        array.set(X, i, StartPoint + i * delta_x)
        array.set(Y, i, f(array.get(X, i)))
    float sum = 0
    for t = 0 to (n - 1) / 4 
        int ind = 4 * t
        sum += (1 / 45.0) * 
             (14 * array.get(Y, ind) 
             + 64 * array.get(Y, ind + 1) 
             + 24 * array.get(Y, ind + 2) 
             + 64 * array.get(Y, ind + 3) 
             + 14 * array.get(Y, ind + 4)) 
             * delta_x
    sum

// N(0,1) cdf by Boole's Rule
N(float x)=> 
    Boole(-10.0, x, 240)

// Black-Scholes Call Price
BSPrice(float S, float K, float r, float T, float q, float v, string PutCall)=>
    float d = (math.log(S / K) + T * (r - q + 0.5 * v * v)) / (v * math.sqrt(T))
    float call = S * math.exp(-q * T) * N(d) - math.exp(-r * T) * K * N(d - v * math.sqrt(T))
    float out = 0
    if (PutCall == "Call")
        out := call
    else
        out := call - S * math.exp(-q * T) + K * math.exp(-r * T)
    out 

// Bisection Algorithm
BisecBSV(float S, float K, float r, float T, float q, float a, float b, float MktPrice, string PutCall)=>
    int MaxIter = 50000
    float Tol = 0.0000001
    float midP = 0
    float midCdif = 0
    float out = 0
    float lowCdif = MktPrice - BSPrice(S, K, r, T, q, a, PutCall)
    float highCdif = MktPrice - BSPrice(S, K, r, T, q, b, PutCall)
    float a1 = a
    float b1 = b
    if lowCdif * highCdif > 0
        out := -1
    else
        for i = 0 to MaxIter 
            midP := (a1 + b1) / 2.0
            midCdif := MktPrice - BSPrice(S, K, r, T, q, midP, PutCall)
            if (math.abs(midCdif) < Tol) 
                break
            else 
                if midCdif > 0
                    a1 := midP
                else 
                    b1 := midP
    midP

ewmavol(float src, int per) => 
    float lambda = (per - 1) / (per + 1)
    float temp = na
    temp := lambda * nz(temp[1], math.pow(src, 2)) + (1.0 - lambda) * math.pow(src, 2)
    out = math.sqrt(temp)
    out
    
rogerssatchel(int per) =>
    float sum = math.sum(math.log(high/ close) * math.log(high / open)
         + math.log(low / close) * math.log(low / open), per) / per
    float out = math.sqrt(sum) 
    out

closetoclose(float src, int per) => 
    float avg = ta.sma(src, per)
    float[] sarr = array.new_float(per, 0)
    for i = 0 to per - 1 
        array.set(sarr, i, math.pow(nz(src[i]) - avg, 2))
    float out = math.sqrt(array.sum(sarr) / (per - 1))
    out 

parkinsonvol(int per)=>
    float volConst = 1.0 / (4.0 * per * math.log(2))
    float sum = volConst * math.sum(math.pow(math.log(high / low), 2), per)
    float out = math.sqrt(sum) 
    out

garmanKlass(int per)=>
    float hllog = math.log(high / low)
    float oplog = math.log(close / open)
    float garmult = (2 * math.log(2) - 1)
    float parkinsonsum = 1 / (2 * per) * math.sum(math.pow(hllog, 2), per)
    float garmansum = garmult / per * math.sum(math.pow(oplog, 2), per)
    float sum = parkinsonsum - garmansum
    float devpercent = math.sqrt(sum)  
    devpercent

gkyzvol(int per)=>
    float gzkylog = math.log(open / nz(close[1]))
    float pklog = math.log(high / low)
    float gklog = math.log(close / open)
    float garmult = (2 * math.log(2) - 1)
    float gkyzsum = 1 / per * math.sum(math.pow(gzkylog, 2), per)
    float parkinsonsum = 1 / (2 * per) * math.sum(math.pow(pklog, 2), per)
    float garmansum = garmult / per * math.sum(math.pow(gklog, 2), per)
    float sum = gkyzsum + parkinsonsum - garmansum
    float devpercent = math.sqrt(sum)  
    devpercent

smthtype = input.string("Kaufman", "Heikin-Ashi Better Caculation Type", options = ["AMA", "T3", "Kaufman"], group = "Source Settings")

srcin = input.string("Close", "Spot Price", group= "Basic Settings", 
     options = 
     ["Close", "Open", "High", "Low", "Median", "Typical", "Weighted", "Average", "Average Median Body", "Trend Biased", "Trend Biased (Extreme)", 
     "HA Close", "HA Open", "HA High", "HA Low", "HA Median", "HA Typical", "HA Weighted", "HA Average", "HA Average Median Body", "HA Trend Biased", "HA Trend Biased (Extreme)",
     "HAB Close", "HAB Open", "HAB High", "HAB Low", "HAB Median", "HAB Typical", "HAB Weighted", "HAB Average", "HAB Average Median Body", "HAB Trend Biased", "HAB Trend Biased (Extreme)"])

float K = input.float(160, "Strike Price", group = "Basic Settings")   
float MktPrice = input.float(194.82, "Option Market Price", group = "Basic Settings")   

string PutCall = input.string("Call", "Option type", options = ["Call", "Put"], group = "Basic Settings") 

string rfrtype = input.string("USD", "Option Base Currency", options = ['USD', 'GBP', 'JPY', 'CAD', 'CNH', 'SGD', 'INR', 'AUD', 'SEK', 'NOK', 'DKK'], group = "Risk-free Rate Settings", tooltip = "Automatically pulls 10-year bond yield from corresponding currency")
float rfrman = input.float(3.97, "% Manual Risk-free Rate", group = "Risk-free Rate Settings") / 100
bool usdrsrman = input.bool(false, "Use manual input for Risk-free Rate?", group = "Risk-free Rate Settings")

float divsman = input.float(7.5, "% Manual Yearly Dividend Yield", group = "Dividend Settings") / 100
bool usediv = input.bool(true, "Adjust for Dividends?", tooltip = "Only works if divdends exist for the current ticker", group = "Dividend Settings")
bool autodiv = input.bool(true, "Automatically Calculate Yearly Dividend Yield?", tooltip = "Only works if divdends exist for the current ticker", group = "Dividend Settings")

int histvolper = input.int(22, "Historical Volatility Period", group = "Volatility Settings", tooltip =  "This is here to use for volatility input.")
string hvoltype = input.string(c2c, "Historical Volatility Type", options = [c2c, gkvol, gkzhvol, rogersatch, ewmavolstr, parkinson], group = "Volatility Settings")

string timein = input.string("Time Now", title = "Time Now Type", options = ["Time Now", "Time Bar", "Trading Day"], group = "Time Intrevals", tooltip = timtoolnow + "; " + timtoolbar + "; " + timetooltrade)
int daysinyear = input.int(252, title = "Days in Year", minval = 1, maxval = 365, group = "Time Intrevals", tooltip = "Typically 252 or 365")
float hoursinday = input.float(24, title = "Hours Per Day", minval = 1, maxval = 24, group = "Time Intrevals", tooltip = "Typically 6.5, 8, or 24")

float a = input.float(0.00000000001, "Bisection Algo starting value for lower volatility", group = "Bisection Algo Settings")
float b = input.float(7.0, "Bisection Algo starting value for higher volatility", group = "Bisection Algo Settings")

int thruMonth = input.int(3, title = "Expiry Month", minval = 1, maxval = 12, group = "Expiry Date/Time")
int thruDay = input.int(31, title = "Expiry Day", minval = 1, maxval = 31, group = "Expiry Date/Time")
int thruYear = input.int(2023, title = "Expiry Year", minval = 1970, group = "Expiry Date/Time")
int mins = input.int(0, title = "Expiry Minute", minval = 0, maxval = 60, group = "Expiry Date/Time")
int hours = input.int(9, title = "Expiry Hour", minval = 0, maxval = 24, group = "Expiry Date/Time")
int secs = input.int(0, title = "Expiry Second", minval = 0, maxval = 60, group = "Expiry Date/Time")

// seconds per year given inputs above
int spyr = math.round(daysinyear * hoursinday * 60 * 60)

// precision calculation miliseconds in time intreval from time equals now
start = timein == "Time Now" ? timenow : timein == "Time Bar" ? time : time_tradingday
finish = timestamp(thruYear, thruMonth, thruDay, hours, mins, secs) 
temp = (finish - start) 
float T = (finish - start) / spyr / 1000

float q = usediv ? (autodiv ? request.dividends(syminfo.tickerid) / close * 4 : divsman) : 0

string byield = switch rfrtype
    "USD"=> 'US10Y'
    "GBP"=> 'GB10Y'
    "JPY"=> 'US10Y'
    "CAD"=> 'CA10Y'
    "CNH"=> 'CN10Y'
    "SGD"=> 'SG10Y'
    "INR"=> 'IN10Y'
    "AUD"=> 'AU10Y'
    "USEKSD"=> 'SE10Y'
    "NOK"=> 'NO10Y'
    "DKK"=> 'DK10Y'
    => 'US10Y'

kfl=input.float(0.666, title="* Kaufman's Adaptive MA (KAMA) Only - Fast End", group = "Moving Average Inputs")
ksl=input.float(0.0645, title="* Kaufman's Adaptive MA (KAMA) Only - Slow End", group = "Moving Average Inputs")
amafl = input.int(2, title="* Adaptive Moving Average (AMA) Only - Fast", group = "Moving Average Inputs")
amasl = input.int(30, title="* Adaptive Moving Average (AMA) Only - Slow", group = "Moving Average Inputs")

haclose = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, close)
haopen = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, open)
hahigh = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, high)
halow = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, low)
hamedian = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, hl2)
hatypical = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, hlc3)
haweighted = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, hlcc4)
haaverage = request.security(ticker.heikinashi(syminfo.tickerid), timeframe.period, ohlc4)

float spot = switch srcin
	"Close" => loxxexpandedsourcetypes.rclose()
	"Open" => loxxexpandedsourcetypes.ropen()
	"High" => loxxexpandedsourcetypes.rhigh()
	"Low" => loxxexpandedsourcetypes.rlow()
	"Median" => loxxexpandedsourcetypes.rmedian()
	"Typical" => loxxexpandedsourcetypes.rtypical()
	"Weighted" => loxxexpandedsourcetypes.rweighted()
	"Average" => loxxexpandedsourcetypes.raverage()
    "Average Median Body" => loxxexpandedsourcetypes.ravemedbody()
	"Trend Biased" => loxxexpandedsourcetypes.rtrendb()
	"Trend Biased (Extreme)" => loxxexpandedsourcetypes.rtrendbext()
	"HA Close" => loxxexpandedsourcetypes.haclose(haclose)
	"HA Open" => loxxexpandedsourcetypes.haopen(haopen)
	"HA High" => loxxexpandedsourcetypes.hahigh(hahigh)
	"HA Low" => loxxexpandedsourcetypes.halow(halow)
	"HA Median" => loxxexpandedsourcetypes.hamedian(hamedian)
	"HA Typical" => loxxexpandedsourcetypes.hatypical(hatypical)
	"HA Weighted" => loxxexpandedsourcetypes.haweighted(haweighted)
	"HA Average" => loxxexpandedsourcetypes.haaverage(haaverage)
    "HA Average Median Body" => loxxexpandedsourcetypes.haavemedbody(haclose, haopen)
	"HA Trend Biased" => loxxexpandedsourcetypes.hatrendb(haclose, haopen, hahigh, halow)
	"HA Trend Biased (Extreme)" => loxxexpandedsourcetypes.hatrendb(haclose, haopen, hahigh, halow)
	"HAB Close" => loxxexpandedsourcetypes.habclose(smthtype, amafl, amasl, kfl, ksl)
	"HAB Open" => loxxexpandedsourcetypes.habopen(smthtype, amafl, amasl, kfl, ksl)
	"HAB High" => loxxexpandedsourcetypes.habhigh(smthtype, amafl, amasl, kfl, ksl)
	"HAB Low" => loxxexpandedsourcetypes.hablow(smthtype, amafl, amasl, kfl, ksl)
	"HAB Median" => loxxexpandedsourcetypes.habmedian(smthtype, amafl, amasl, kfl, ksl)
    "HAB Typical" => loxxexpandedsourcetypes.habtypical(smthtype, amafl, amasl, kfl, ksl)
	"HAB Weighted" => loxxexpandedsourcetypes.habweighted(smthtype, amafl, amasl, kfl, ksl)
	"HAB Average" => loxxexpandedsourcetypes.habaverage(smthtype, amafl, amasl, kfl, ksl)
    "HAB Average Median Body" => loxxexpandedsourcetypes.habavemedbody(smthtype, amafl, amasl, kfl, ksl)
	"HAB Trend Biased" => loxxexpandedsourcetypes.habtrendb(smthtype, amafl, amasl, kfl, ksl)
	"HAB Trend Biased (Extreme)" => loxxexpandedsourcetypes.habtrendbext(smthtype, amafl, amasl, kfl, ksl)
	=> haclose

float hvolout = switch hvoltype
    parkinson => parkinsonvol(histvolper)
    rogersatch => rogerssatchel(histvolper) 
    c2c => closetoclose(math.log(spot / nz(spot[1])), histvolper) 
    gkvol => garmanKlass(histvolper) 
    gkzhvol => gkyzvol(histvolper) 
    ewmavolstr => ewmavol(math.log(spot / nz(spot[1])), histvolper) 

float r = usdrsrman ? rfrman : request.security(byield, timeframe.period, spot) / 100

if barstate.islast
    IV = BisecBSV(spot, K, r, T, q, a, b, MktPrice, PutCall)
    var testTable = table.new(position = position.middle_right, columns = 1, rows = 16, bgcolor = color.yellow, border_width = 1)

    table.cell(table_id = testTable, column = 0, row = 0,  text = "  Inputs  ", bgcolor=color.yellow, text_color = color.black)
    table.cell(table_id = testTable, column = 0, row = 1,  text = "  Spot Price: " + str.tostring(spot, f_tickFormat()) , bgcolor=darkGreenColor, text_color = color.white)
    table.cell(table_id = testTable, column = 0, row = 2,  text = "  Strike Price: " + str.tostring(K, f_tickFormat()), bgcolor=darkGreenColor, text_color = color.white)
    table.cell(table_id = testTable, column = 0, row = 3,  text = "  Risk-free Rate Type: " + rfrtype , bgcolor=darkGreenColor, text_color = color.white)
    table.cell(table_id = testTable, column = 0, row = 4,  text = "  Risk-free Rate: " + str.tostring(r * 100, "##.##") + "%  ", bgcolor=darkGreenColor, text_color = color.white)
    table.cell(table_id = testTable, column = 0, row = 5,  text = "  Dividend Yield (annual): " + str.tostring(q * 100, "##.##") + "%  ", bgcolor=darkGreenColor, text_color = color.white)
    table.cell(table_id = testTable, column = 0, row = 6,  text = "  Option Type: " + PutCall, bgcolor=darkGreenColor, text_color = color.white)
    table.cell(table_id = testTable, column = 0, row = 7,  text = "  Market Price: " + str.tostring(MktPrice, f_tickFormat()), bgcolor=darkGreenColor, text_color = color.white)

    table.cell(table_id = testTable, column = 0, row = 8,  text = "  Time Now: " + str.format("{0,date,MMMM dd, yyyy - HH:mm:ss}", timenow), bgcolor=darkGreenColor, text_color = color.white)
    table.cell(table_id = testTable, column = 0, row = 9,  text = "  Expiry Date: " + str.format("{0,date,MMMM dd, yyyy - HH:mm:ss}", finish), bgcolor=darkGreenColor, text_color = color.white)

    table.cell(table_id = testTable, column = 0, row = 10,  text = "  Output  ", bgcolor=color.yellow, text_color = color.black)
    table.cell(table_id = testTable, column = 0, row = 11, text = "  Implied Volatility: " + str.tostring(IV * 100, "##.###") + "%", bgcolor=darkGreenColor, text_color = color.white)

    table.cell(table_id = testTable, column = 0, row = 12, text = "  Calculated Values  ", bgcolor=color.yellow, text_color = color.black)
    table.cell(table_id = testTable, column = 0, row = 13, text = "  Hist. Volatility Type: " + hvoltype, bgcolor=darkGreenColor, text_color = color.white)
    table.cell(table_id = testTable, column = 0, row = 14, text = "  Hist. Daily Volatility: " + str.tostring(hvolout * 100, "##.##") + "%  ", bgcolor=darkGreenColor, text_color = color.white)
    table.cell(table_id = testTable, column = 0, row = 15, text = "  Hist. Annualized Volatility: " + str.tostring(hvolout * math.sqrt(daysinyear) * 100, "##.##") + "%  ", bgcolor = darkGreenColor, text_color = color.white)
