QQE Student's T-Distribution Bollinger Bands Screener
// 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("QQETSBB Screener", overlay = true, precision = 0, max_bars_back = 1000)

// The number 2 * sqrt(e / pi) number:
var float math_TwoSqrtEOverPi = 1.8603827342052657173362492472666631120594218414085755

// The order of the Log Gamma aproximation:
var int GammaN = 10

// Auxiliary variable when evaluating the Log Gamma function:
var float GammaR = 10.900511

// Series of Polynomial coefficients for the Gamma approximation:
var float[] series_GammaDk = 
 array.from(2.48574089138753565546e-5, 
 1.05142378581721974210, 
 -3.45687097222016235469, 
 4.51227709466894823700, 
 -2.98285225323576655721, 
 1.05639711577126713077, 
 -1.95428773191645869583e-1, 
 1.70970543404441224307e-2, 
 -5.71926117404305781283e-4, 
 4.63399473359905636708e-6, 
 -2.71994908488607703910e-9)

gamma(input) =>
    // <needs>       : Precision.Constants
    //
    // <param> input : Any float input as argument of the Gamma Function.
    //
    // <returns>     : The Gamma Function.

    float z = input

    if z < 0.5

        float s = array.get(series_GammaDk, 0)
        for int i = 1 to GammaN by 1
            s += array.get(series_GammaDk, i) / (i - z)
            s

        return_Gamma = math.pi / math.sin(math.pi * z) * s * math_TwoSqrtEOverPi * math.pow((0.5 - z + GammaR) / math.e, 0.5 - z)
        return_Gamma
    else

        float s = array.get(series_GammaDk, 0)
        for int i = 1 to GammaN by 1
            s += array.get(series_GammaDk, i) / (z + i - 1.0)
            s

        return_Gamma = s * math_TwoSqrtEOverPi * math.pow((z - 0.5 + GammaR) / math.e, z - 0.5)
        return_Gamma

// EndOf: function_Gamma() }
ha_close(float Open = open, float High = high, float Low = low, float Close = close, bool enable = true) =>
    ha_close = (Open + High + Low + Close) / 4
    out      = enable == true ? ha_close : Close
 
ha_open(float Open = open, float High = high, float Low = low, float Close = close, bool enable = true) =>
    ha_open  = float(na)
    ha_close = ha_close(Open, High, Low, Close)
    ha_open := na(ha_open[1]) ? (Open + Close) / 2 : (nz(ha_open[1]) + nz(ha_close[1])) / 2 
    out      = enable == true ? ha_open : Open

diff(_src, _N) =>
    float derivative = 0.0
    if _N == 2
        derivative := (_src - _src[2]) / 2
        derivative
    else if _N == 3
        derivative := (_src + _src[1] - _src[2] - _src[3]) / 4
        derivative
    else if _N == 4
        derivative := (_src + 2 * _src[1] - 2 * _src[3] - _src[4]) / 8
        derivative
    else if _N == 5
        derivative := (_src + 3 * _src[1] + 2 * _src[2] - 2 * _src[3] - 3 * _src[4] - _src[5]) / 16
        derivative
    else if _N == 6
        derivative := (_src + 4 * _src[1] + 5 * _src[2] - 5 * _src[4] - 4 * _src[5] - _src[6]) / 32
        derivative
    else if _N == 7
        derivative := (_src + 5 * _src[1] + 9 * _src[2] + 5 * _src[3] - 5 * _src[4] - 9 * _src[5] - 5 * _src[6] - _src[7]) / 64
        derivative
    else if _N == 8
        derivative := (_src + 6 * _src[1] + 14 * _src[2] + 14 * _src[3] - 14 * _src[5] - 14 * _src[6] - 6 * _src[7] - _src[8]) / 128
        derivative
    else if _N == 9
        derivative := (_src + 7 * _src[1] + 20 * _src[2] + 28 * _src[3] + 14 * _src[4] - 14 * _src[5] - 28 * _src[6] - 20 * _src[7] - 7 * _src[8] - _src[9]) / 256
        derivative
    else if _N == 10
        derivative := (_src + 8 * _src[1] + 27 * _src[2] + 48 * _src[3] + 42 * _src[4] - 42 * _src[6] - 48 * _src[7] - 27 * _src[8] - 8 * _src[9] - _src[10]) / 512
        derivative
    derivative

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


style      = input.string("QQE", "Moving Average", ["QQE","SMA","EMA","WMA"], group="Deviation")
location = input.string("Middle Center", "Position", 
 ["Bottom Left","Bottom Middle","Bottom Right","Middle Left","Middle Center","Middle Right",
 "Top Left","Top Middle","Top Right"])
len        = input.int(200, title="Length", minval=2, tooltip="How long the rolling window for lookback calculations should be", group="Deviation")
itter      = input.int(0, "Cumulative Distribution Function Iterations", 0, group="Deviation")
geometric  = input.bool(true, 'Log-space', group="Features")
diff_n     = input.int(2, "Direction Look Back", 2, 10, group="Features")

smoothing  = input.int(title='QQE Smoothing (when QQE is used as signal line)', defval=5, minval=1, group="QQE")
qqe_factor = input.float(title='QQE Factor (when QQE is used as signal line)', defval=5, minval=1, group="QQE")

custom     = input.bool(false, "Custom Time Frame", inline = "Time", tooltip = "When using custom time frame use the minimum time on your chart.")
ctime      = input.timeframe("D", "", inline = "Time")
txt_col    = input.color(color.new(#c8d7ff, 0), "Title Color", inline = "bull", group = "Settings")
box_col    = input.color(color.new(#202236, 10), "", inline = "bull", group = "Settings")

location() => 
    switch location
        "Bottom Left" => position.bottom_left
        "Bottom Middle" => position.bottom_center
        "Bottom Right" => position.bottom_right
        "Middle Left" => position.middle_left
        "Middle Center" => position.middle_center
        "Middle Right" => position.middle_right
        "Top Left" => position.top_left
        "Top Middle" => position.top_center
        "Top Right" => position.top_right

hypergeometric(float a, float b, float c, float z) =>
    float sum = 0.0
    for i = 0 to itter
        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) =>
    // I hate myself for doing this
    1/2. + x * gamma((v+1)/2.) * (hypergeometric(1/2., (v+1)/2., 3./2, x*x/v) / (math.sqrt(v * math.pi)* gamma(v/2)))

ext1  = input.symbol('BTCUSD', 'Symbol 1', group = "Symbols")
ext2  = input.symbol('ETHUSD', 'Symbol 2', group = "Symbols")
ext3  = input.symbol('BNBUSD', 'Symbol 3', group = "Symbols")
ext4  = input.symbol('XRPUSD', 'Symbol 4', group = "Symbols")
ext5  = input.symbol('DOGEUSD', 'Symbol 5', group = "Symbols")
ext6  = input.symbol('ETCUSD', 'Symbol 6', group = "Symbols")
ext7  = input.symbol('LTCUSD', 'Symbol 7', group = "Symbols")
ext8  = input.symbol('XLMUSD', 'Symbol 8', group = "Symbols")
ext9  = input.symbol('XMRUSD', 'Symbol 9', group = "Symbols")
ext10 = input.symbol('TRXUSD', 'Symbol 10', group = "Symbols")
ext11 = input.symbol('AAPL', 'Symbol 11', group = "Symbols")
ext12 = input.symbol('MSFT', 'Symbol 12', group = "Symbols")
ext13 = input.symbol('GOOG', 'Symbol 13', group = "Symbols")
ext14 = input.symbol('AMZN', 'Symbol 14', group = "Symbols")
ext15 = input.symbol('TSLA', 'Symbol 15', group = "Symbols")
ext16 = input.symbol('UNH', 'Symbol 16', group = "Symbols")
ext17 = input.symbol('META', 'Symbol 17', group = "Symbols")
ext18 = input.symbol('NVDA', 'Symbol 18', group = "Symbols")
ext19 = input.symbol('V', 'Symbol 19', group = "Symbols")
ext20 = input.symbol('MA', 'Symbol 20', group = "Symbols")

timep = custom ? ctime : timeframe.period

var tbl = table.new(location(), 7, 21)
if barstate.islast
    table.cell(tbl, 0, 0, "Symbol", text_color = txt_col, bgcolor = box_col, tooltip = "Symbol for that line.")
    table.cell(tbl, 1, 0, "Price", text_color = txt_col, bgcolor = box_col, tooltip = "Price of Symbol")
    table.cell(tbl, 2, 0, "Band Spread", text_color = txt_col, bgcolor = box_col, tooltip = "The volitility range in %.")
    table.cell(tbl, 3, 0, "Segment Spread", text_color = txt_col, bgcolor = box_col, tooltip = "1/4 the volitility range in %.")
    table.cell(tbl, 4, 0, "Level", text_color = txt_col, bgcolor = box_col, tooltip = "What quarter of the band the price is currently at.")
    table.cell(tbl, 5, 0, "Sentiment", text_color = txt_col, bgcolor = box_col, tooltip = "A score bassed on all of the popular indicators. If the score is high, the popular indicators are showing over baught and visa versa. This is a combination of rsi, macd, stochastic, cci, mfi, and otheres.")
    table.cell(tbl, 6, 0, "Direction", text_color = txt_col, bgcolor = box_col, tooltip = "This is the direction the price is moving in.")

qqe(source) =>
    wilders_period = len * 2 - 1

    ma = ta.ema(source, smoothing)

    atr            = math.abs(nz(ma[1]) - ma)
    ma_atr         = ta.ema(atr, wilders_period)
    delta_fast_atr = ta.ema(ma_atr, wilders_period) * qqe_factor

    result = 0.0
    result := ma > nz(result[1]) ? ma - delta_fast_atr < result[1] ? result[1] : ma - delta_fast_atr : ma + delta_fast_atr > result[1] ? result[1] : ma + delta_fast_atr

    // Return
    result

ma(source) =>
    switch style
        "SMA" => ta.sma(source, len)
        "EMA" => ta.ema(source, len)
        "WMA" => ta.wma(source, len)
        "QQE" => qqe(source)

main(source, n) =>
    exn = str.tostring(source)
    ex  = request.security(source , timep, close)
    exh = ta.highest(ex, 1)
    exl = ta.lowest(ex, 1)
    // Calculate p-value given the confidence
    p = (1 - 0.9999) / 2
    // We need to find the value of x such that t_dist_cdf(x, length - 1) = p
    // We do this by using a binary search
    src = geometric ? math.log(ex) : ex
    float lower = -100
    float upper = 100
    float mid   = (lower + upper) / 2
    while math.abs(t_dist_cdf(mid, len - 1) - p) > 0.000001
        if t_dist_cdf(mid, len - 1) > p
            upper := mid
        else
            lower := mid
        mid := (lower + upper) / 2.
    float t_value = mid
    // t-value for the p-value
    if t_value < 0
        t_value := -t_value

    x_bar          = ma(src)
    squared_diff   = math.pow(src-x_bar, 2)
    summed_squares = math.sum(squared_diff, len)
    // length - 1 because we're using a sample rather than whole population
    mean_squared_error  = summed_squares / (len-1)
    rmse                = math.sqrt(mean_squared_error)
    root_standard_error = math.sqrt(1 + 1/len + squared_diff/summed_squares)
    offset              = t_value * rmse * root_standard_error

    // The midline
    mid_line = x_bar
    // Calculate the upper and lower bounds
    upper_bound = x_bar + offset
    lower_bound = x_bar - offset

    Bottom      = geometric ? math.exp(lower_bound) : lower_bound
    Top         = geometric ? math.exp(upper_bound) : upper_bound
    Center      = geometric ? math.exp(x_bar) : x_bar
    Signal      = geometric ? math.exp(src) : src
    Half_Bottom = Bottom - ((Bottom-Center)/2)
    Half_Top    = Top - ((Top-Center)/2)

    mean_diff  = diff(Center, diff_n)
    Diff_Score = mean_diff < 0 ? "↓" : mean_diff > 0 ? "↑" : "→"

    Score = Signal >= Top ? 5 : Signal > Half_Top ? 4 : Signal > Center ? 3 : Signal > Half_Bottom ? 2 : Signal > Bottom ? 1 : Signal <= Bottom ? 0 : -1

    Band    = (((Top - Bottom) / Bottom) - (((Top - Bottom) / Top) % 0.001)) * 100
    Segment = (((Half_Top - Center) / Center) - (((Half_Top - Center) / Center) % 0.001)) * 100

    oc_trend = ta.sma(ta.ema(ex, 20) - ta.ema(ex[1], 20), 2)
    oc_score = oc_trend > 0 ? 1 : 0

    fast_e    = ta.ema(ex, 9)
    slow_e    = ta.ema(ex, 21)
    ema_cross = fast_e > slow_e ? 1 : 0

    fast_s    = ta.sma(ex, 13)
    slow_s    = ta.sma(ex, 48)
    sma_cross = fast_s > slow_s ? 1 : 0

    rsi       = ta.rsi(ex, 14)
    rsi_score = rsi < 30 ? 0 : rsi > 30 and rsi < 70 ? 0.5 : rsi > 70 ? 1 : 0

    [macd1, macd2, hist] = ta.macd(ex, 12, 26, 9)
    macd_score           = macd1 > macd2 ? 0.5 : macd1 < macd2 ? 0 : 0
    hist_score           = (hist >=0 ? (hist[1] < hist ? 0.5 : 0.25) : (hist[1] < hist ? 0.25 : 0))

    stoch       = ta.sma(ta.stoch(ex, exh, exl, 14), 2)
    stoch_score = stoch < 20 ? 0 : stoch > 20 and stoch < 80 ? 0.5 : stoch > 80 ? 1 : 0

    ma200       = ta.sma(ex, 200)
    ma200_score = ex > ma200 ? 1 : ex <= ma200 ? 0 : 0

    mfi       = ta.mfi(ex, 14)
    mfi_score = mfi < 20 ? 0 : mfi > 20 and mfi < 80 ? 0.5 : mfi > 80 ? 1 : 0 

    cci       = ta.cci(ex, 20)
    cci_score = cci < -100 ? 0 : cci > -100 and cci < 100 ? 0.5 : cci > 100 ? 1 : 0

    ha_open  = ha_open(ex[1], exh, exl, ex)
    ha_close = ha_close(ex[1], exh, exl, ex)
    ha_score = ha_open < ha_close ? 1 : 0

    sent      = ema_cross + sma_cross + rsi_score + macd_score + hist_score + stoch_score + ma200_score + mfi_score + cci_score + oc_score + ha_score
    sentiment = sent

    if barstate.islast
        table.cell(tbl, 0, n, exn, text_color = txt_col, bgcolor = box_col)
        table.cell(tbl, 1, n, str.tostring(ex), text_color = txt_col, bgcolor = box_col)
        table.cell(tbl, 2, n, str.tostring(Band - (Band % 0.1)) + "%", text_color = txt_col, bgcolor = box_col)
        table.cell(tbl, 3, n, str.tostring(Segment - (Segment % 0.1)) + "%", text_color = txt_col, bgcolor = box_col)
        table.cell(tbl, 4, n, str.tostring(Score) + "/5", text_color = txt_col, bgcolor = box_col )
        table.cell(tbl, 5, n, str.tostring(sentiment) + "/10", text_color = txt_col, bgcolor = box_col )
        table.cell(tbl, 6, n, Diff_Score, text_color = txt_col, bgcolor = box_col)

main(ext1, 1)
main(ext2, 2)
main(ext3, 3)
main(ext4, 4)
main(ext5, 5)
main(ext6, 6)
main(ext7, 7)
main(ext8, 8)
main(ext9, 9)
main(ext10, 10)
main(ext11, 11)
main(ext12, 12)
main(ext13, 13)
main(ext14, 14)
main(ext15, 15)
main(ext16, 16)
main(ext17, 17)
main(ext18, 18)
main(ext19, 19)
main(ext20, 20)