Price Level Stats (PLS)
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © peacefulLizard50262

//@version=5
indicator("Price Level Stats", "PLS", precision = 3, overlay = true)

import RicardoSantos/MathSpecialFunctionsGamma/1

// Types
type vec2
    float x
    float y

// Helper functions
gamma(float x) => MathSpecialFunctionsGamma.Gamma(x)

logReturn(float a, float b) => math.log(a / b)
percent(float a, float b) => (a - b) / b * 100
addLogReturn(float a, float b) => a * math.exp(b)
addPercentage(float a, float b) => a * (1 + b / 100)
subtractPercentage(float a, float b) => a * (1 - b / 100)
stepToInt(float stepSize) => int(stepSize * 100)

ema(float source = close, float length = 9)=>
    alpha = 2 / (length + 1)
    var float smoothed = na
    smoothed := alpha * source + (1 - alpha) * nz(smoothed[1])

t_dist_pdf(float x, int v) =>
    gamma((v + 1) / 2) / (math.sqrt(v * math.pi) * gamma(v / 2)) * math.pow(1 + x * x / v, -(v + 1) / 2)

rising_factorial(float c, int n) =>
    float product = 1.0
    if n > 0
        for i = 0 to n - 1
            product *= (c + i)
    product

hypergeometric(float a, float b, float c, float z) =>
    float sum = 0.0
    for i = 0 to 2
        sum += rising_factorial(a, i) * rising_factorial(b, i) * math.pow(z, i) / (rising_factorial(c, i) * rising_factorial(i, i))
    sum

t_dist_cdf(float x, int v) =>
    1/2.0 + x * gamma((v+1)/2.0) * (hypergeometric(1/2.0, (v+1)/2.0, 3.0/2, x*x/v) / (math.sqrt(v * math.pi)* gamma(v/2)))

// Price Movement Probability (T-Dist) function
priceMovementProbabilityTDist(float inputPrice, float mean, series int meanLength) =>
    float p = (1 - 0.99) / 2
    float lower = -100
    float upper = 100
    float mid = (lower + upper) / 2
    while math.abs(t_dist_cdf(mid, meanLength - 1) - p) > 0.000001
        if t_dist_cdf(mid, meanLength - 1) > p
            upper := mid
        else
            lower := mid
        mid := (lower + upper) / 2
    float t_value = mid
    if t_value < 0
        t_value := -t_value
    float squared_diff = math.pow(open-mean, 2)
    float summed_squares = math.sum(squared_diff, meanLength)
    float mean_squared_error = summed_squares / (meanLength-1)
    float rmse = math.sqrt(mean_squared_error)
    float root_standard_error = math.sqrt(1 + 1/meanLength + squared_diff/summed_squares)
    float offset = t_value * rmse * root_standard_error
    float zscore = (inputPrice-mean)/offset
    float probability = t_dist_pdf(zscore, meanLength - 1)
    probability

drawProbability(vec2 source, bool logReturns) =>
    price = array.new<label>(1)
    lines = array.new<line>(1)
    if barstate.islast
        array.push(price, label.new(bar_index + 2, source.x, text = str.tostring(math.round(source.y * (logReturns ? 1 : 100), 2))+"%", style=label.style_label_down, textcolor = source.x > open ? color.green : color.red, color=color.new(color.white, 100)))
        array.push(lines, line.new(bar_index, source.x, bar_index + 1, source.x, color = source.x > open ? color.green : color.red))
    if barstate.isconfirmed
        for i = 0 to array.size(price) - 1
            label.delete(array.get(price, i))
            line.delete(array.get(lines, i))

source = input.source(open, "Source")
stepSize = stepToInt(input.float(0.2, "Step Size", 0.001, 100, 0.001))
steps = input.int(1, "Steps", 1, 100, 1)
meanLength = input.int(20, "Mean/Deviation Length", 5)
style = input.string("Bar Estimate", "Style", ["Bar Estimate", "Log Bar Estimate", "Student's T"])
useLogReturns = style == "Log Bar Estimate" ? true : false

var movement_size_green = array.new<float>(100)
var count_green = array.new<int>(100)
var movement_size_red = array.new<float>(100)
var count_red = array.new<int>(100)
change = useLogReturns ? math.round(logReturn(close, close[1]), 4) : math.round(percent(close, close[1]), 4)

float mean = ema(open, meanLength)

if close > open
    if array.includes(movement_size_green, change)
        index = array.indexof(movement_size_green, change)
        array.set(count_green, index, array.get(count_green, index) + 1)
    else
        array.push(movement_size_green, change)
        array.push(count_green, 1)

if close < open
    if array.includes(movement_size_red, change)
        index = array.indexof(movement_size_red, change)
        array.set(count_red, index, array.get(count_red, index) + 1)
    else
        array.push(movement_size_red, change)
        array.push(count_red, 1)

probability_green = array.new<float>(array.size(count_green))
probability_red   = array.new<float>(array.size(count_red))

count_green_sum = array.sum(count_green)
count_red_sum = array.sum(count_red)

for i = 0 to array.size(count_green) - 1
    probability = array.get(count_green, i) / count_green_sum * 100
    array.set(probability_green, i, probability)

for i = 0 to array.size(count_red) - 1
    probability = array.get(count_red, i) / count_red_sum * 100
    array.set(probability_red, i, probability)

if style != "Student's T"
    for i = stepSize to stepSize * steps by stepSize
        inputPriceBullish = addPercentage(source, i/100.0)
        inputPcireBearish = subtractPercentage(source, i/100.0)
        arb_price_bullish = useLogReturns ? math.round(logReturn(inputPriceBullish, open), 4) : math.round(percent(inputPriceBullish, open), 4)
        arb_price_bearish = useLogReturns ? math.round(logReturn(inputPcireBearish, open), 4) : math.round(percent(inputPcireBearish, open), 4)
        green_index_up    = array.indexof(movement_size_green, arb_price_bullish)
        red_index_up      = array.indexof(movement_size_red, arb_price_bearish)
        green_index_down  = array.indexof(movement_size_green, arb_price_bullish)
        red_index_down    = array.indexof(movement_size_red, arb_price_bearish)
        probability_up    = inputPriceBullish > 0 and green_index_up != -1 ? array.get(probability_green, green_index_up) : inputPriceBullish < 0 and red_index_up != -1 ? array.get(probability_red, red_index_up) : 0

        probability_down = inputPcireBearish > 0 and green_index_down != -1 ? array.get(probability_green, green_index_down) : inputPcireBearish < 0 and red_index_down != -1 ? array.get(probability_red, red_index_down) : 0

        drawProbability(vec2.new(inputPriceBullish, probability_up), useLogReturns)
        drawProbability(vec2.new(inputPcireBearish, probability_down), useLogReturns)
else
    for i = stepSize to stepSize * steps by stepSize
        inputPriceBullish = addPercentage(source, i/100.0)
        inputPcireBearish = subtractPercentage(source, i/100.0)
        probability_up    = priceMovementProbabilityTDist(inputPriceBullish, mean, meanLength)
        probability_down  = priceMovementProbabilityTDist(inputPcireBearish, mean, meanLength)

        drawProbability(vec2.new(inputPriceBullish, probability_up), false)
        drawProbability(vec2.new(inputPcireBearish, probability_down), false)