Polynomial Regression Bands w/ Extrapolation of Price [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("Polynomial Regression Bands w/ Extrapolation of Price [Loxx]", 
     shorttitle ="PRBEP [Loxx]",
     overlay = true, 
     max_lines_count = 500)

import loxx/loxxexpandedsourcetypes/4
import loxx/loxxmas/1

greencolor = #2DD204
redcolor = #D2042D 

lightgreencolor = #96E881 
lightredcolor =  #DF4F6C

polyfit(float[] src, int deg, int len, int bar)=>
    sumout = 0.
    AX = matrix.new<float>(12, 12, 0.)
    BX = array.new<float>(12, 0.)
    ZX = array.new<float>(12, 0.)
    Pow = array.new<float>(12, 0.)
    Row = array.new<int>(12, 0)
    CX = array.new<float>(12, 0.)
    
    // FILL MATRIX FIRST
    if (len <= 1) 
        sumout := array.get(src, bar)
    else
        XK = 0.
        Prod = 0.
        for j = 1 to deg + 1
            array.set(BX, j, 0)

        for k = 1 to len 
            YK = array.get(src, len - k)
            XK := k
            Prod := 1
            for j = 1 to deg + 1
                array.set(BX, j, array.get(BX, j) + YK * Prod)
                Prod *= XK
    
        for j = 0 to 2 * deg
            array.set(Pow, j, 0)

        array.set(Pow, 0, len)
    
        for k = 1 to len  
            XK := k
            Prod := k
            for j = 1 to 2 * deg 
                array.set(Pow, j, array.get(Pow, j) + Prod)
                Prod *= XK
    
        for j = 1 to deg + 1 
            for l = 1 to deg + 1 
                matrix.set(AX, j, l, array.get(Pow, j + l - 2))
        
        //NOW SOLVE FOR COEFFICIENTS
        for j = 1 to deg + 1 
            array.set(Row, j, j)
    
        for i = 1 to deg 
            for k = i + 1 to deg + 1
                if math.abs(matrix.get(AX, array.get(Row, k), i)) > 
                     math.abs(matrix.get(AX, array.get(Row, i), i))
                     
                    temp = array.get(Row, i)
                    array.set(Row, i, array.get(Row, k))
                    array.set(Row, k, temp)
                    
            for k = i + 1 to  deg + 1  
                if matrix.get(AX, array.get(Row, i), i) != 0
                    matrix.set(AX, array.get(Row, k), i, 
                         matrix.get(AX, array.get(Row, k), i) / 
                         matrix.get(AX, array.get(Row, i), i))
                        
                for l = i + 1 to deg + 1  
                    matrix.set(AX, array.get(Row, k), l, 
                         matrix.get(AX, array.get(Row, k), l) - 
                         matrix.get(AX, array.get(Row, k), i) * 
                         matrix.get(AX, array.get(Row, i), l))
                         
        array.set(ZX, 1, array.get(BX, array.get(Row, 1)))
        for k = 2 to deg + 1 
            sum = 0.
            for l = 1 to k - 1 
                sum += matrix.get(AX, array.get(Row, k), l) * array.get(ZX, l)
            array.set(ZX, k, array.get(BX, array.get(Row, k)) - sum)
    
        array.set(CX, deg + 1, array.get(ZX, deg + 1) / matrix.get(AX, array.get(Row, deg + 1), deg + 1))
        for k = deg to 1
            sum = 0.
            for l = k + 1 to deg + 1  
                sum += matrix.get(AX, array.get(Row, k), l) * array.get(CX, l)
            array.set(CX, k, (array.get(ZX, k) - sum) / matrix.get(AX, array.get(Row, k), k))
            
        //NOW COMPUTE NEXT POINT IN SERIES AND RETURN
        sumout := array.get(CX, deg + 1)
        for k = deg to 1
            sumout := array.get(CX, k) + sumout * (len + bar)
    sumout

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

srcin = input.string("Close", "Source", group= "Source 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)"])

Length = input.int(30, "Period", group = "Basic Settings")
Degree = input.int(3, "Degree", group = "Basic Settings")
MA_Length = input.int(5, "Source Smoothing Period", group = "Basic Settings")

type = input.string("Exponential Moving Average - EMA", "Source Smoothing Type", options = ["ADXvma - Average Directional Volatility Moving Average",  "Ahrens Moving Average"
     , "Alexander Moving Average - ALXMA", "Double Exponential Moving Average - DEMA", "Double Smoothed Exponential Moving Average - DSEMA"
     , "Exponential Moving Average - EMA", "Fast Exponential Moving Average - FEMA", "Fractal Adaptive Moving Average - FRAMA"
     , "Hull Moving Average - HMA", "IE/2 - Early T3 by Tim Tilson", "Integral of Linear Regression Slope - ILRS"
     , "Instantaneous Trendline", "Laguerre Filter", "Leader Exponential Moving Average", "Linear Regression Value - LSMA (Least Squares Moving Average)"
     , "Linear Weighted Moving Average - LWMA", "McGinley Dynamic", "McNicholl EMA", "Non-Lag Moving Average", "Parabolic Weighted Moving Average"
     , "Recursive Moving Trendline", "Simple Moving Average - SMA", "Sine Weighted Moving Average", "Smoothed Moving Average - SMMA"
     , "Smoother", "Super Smoother", "Three-pole Ehlers Butterworth", "Three-pole Ehlers Smoother"
     , "Triangular Moving Average - TMA", "Triple Exponential Moving Average - TEMA", "Two-pole Ehlers Butterworth", "Two-pole Ehlers smoother"
     , "Volume Weighted EMA - VEMA", "Zero-Lag DEMA - Zero Lag Double Exponential Moving Average", "Zero-Lag Moving Average"
     , "Zero Lag TEMA - Zero Lag Triple Exponential Moving Average"],
     group = "Basic Settings")
     
ProjLength = input.int(10, maxval = 99, minval = 5, group = "Basic Settings")
K_Sigma = input.float(3, "Inner Bands Multiplier", group = "Basic Settings")
K_Sigma2 = input.float(6, "Outer Bands Multiplier", group = "Basic Settings")

FitMode = input.string("Moving", "Fit Mode", options =["Moving", "Fitting"], group = "Basic Settings")
barsbark = input.int(600, "Bars back", group = "Basic Settings", tooltip = "Number of bars painted on screen. This controls the computation speed")

colorbars = input.bool(true, "Color bars?", group = "UI Options")
showbands = input.bool(true, "Show bands?", group = "UI Options")
showmid = input.bool(true, "Show center line?", group = "UI Options")

frama_FC = input.int(defval=1, title="* Fractal Adjusted (FRAMA) Only - FC", group = "Moving Average Inputs")
frama_SC = input.int(defval=200, title="* Fractal Adjusted (FRAMA) Only - SC", group = "Moving Average Inputs")
instantaneous_alpha = input.float(defval=0.07, minval = 0, title="* Instantaneous Trendline (INSTANT) Only - Alpha", group = "Moving Average Inputs")
_laguerre_alpha = input.float(title="* Laguerre Filter (LF) Only - Alpha", minval=0, maxval=1, step=0.1, defval=0.7, group = "Moving Average Inputs")
lsma_offset = input.int(defval=0, title="* Least Squares Moving Average (LSMA) Only - Offset", group = "Moving Average Inputs")
_pwma_pwr = input.int(2, "* Parabolic Weighted Moving Average (PWMA) Only - Power", minval=0, group = "Moving Average Inputs")
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)

src = 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.hatrendbext(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
variant(type, src, len) =>
    sig = 0.0
    trig = 0.0
    special = false
    if type == "ADXvma - Average Directional Volatility Moving Average"
        [t, s, b] = loxxmas.adxvma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Ahrens Moving Average"
        [t, s, b] = loxxmas.ahrma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Alexander Moving Average - ALXMA"
        [t, s, b] = loxxmas.alxma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Double Exponential Moving Average - DEMA"
        [t, s, b] = loxxmas.dema(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Double Smoothed Exponential Moving Average - DSEMA"
        [t, s, b] = loxxmas.dsema(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Exponential Moving Average - EMA"
        [t, s, b] = loxxmas.ema(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Fast Exponential Moving Average - FEMA"
        [t, s, b] = loxxmas.fema(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Fractal Adaptive Moving Average - FRAMA"
        [t, s, b] = loxxmas.frama(src, len, frama_FC, frama_SC)
        sig := s
        trig := t
        special := b
    else if type == "Hull Moving Average - HMA"
        [t, s, b] = loxxmas.hma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "IE/2 - Early T3 by Tim Tilson"
        [t, s, b] = loxxmas.ie2(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Integral of Linear Regression Slope - ILRS"
        [t, s, b] = loxxmas.ilrs(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Instantaneous Trendline"
        [t, s, b] = loxxmas.instant(src, instantaneous_alpha)
        sig := s
        trig := t
        special := b
    else if type == "Laguerre Filter"
        [t, s, b] = loxxmas.laguerre(src, _laguerre_alpha)
        sig := s
        trig := t
        special := b
    else if type == "Leader Exponential Moving Average"
        [t, s, b] = loxxmas.leader(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Linear Regression Value - LSMA (Least Squares Moving Average)"
        [t, s, b] = loxxmas.lsma(src, len, lsma_offset)
        sig := s
        trig := t
        special := b
    else if type == "Linear Weighted Moving Average - LWMA"
        [t, s, b] = loxxmas.lwma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "McGinley Dynamic"
        [t, s, b] = loxxmas.mcginley(src, len)
        sig := s
        trig := t
        special := b
    else if type == "McNicholl EMA"
        [t, s, b] = loxxmas.mcNicholl(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Non-Lag Moving Average"
        [t, s, b] = loxxmas.nonlagma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Parabolic Weighted Moving Average"
        [t, s, b] = loxxmas.pwma(src, len, _pwma_pwr)
        sig := s
        trig := t
        special := b
    else if type == "Recursive Moving Trendline"
        [t, s, b] = loxxmas.rmta(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Simple Moving Average - SMA"
        [t, s, b] = loxxmas.sma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Sine Weighted Moving Average"
        [t, s, b] = loxxmas.swma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Smoothed Moving Average - SMMA"
        [t, s, b] = loxxmas.smma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Smoother"
        [t, s, b] = loxxmas.smoother(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Super Smoother"
        [t, s, b] = loxxmas.super(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Three-pole Ehlers Butterworth"
        [t, s, b] = loxxmas.threepolebuttfilt(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Three-pole Ehlers Smoother"
        [t, s, b] = loxxmas.threepolesss(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Triangular Moving Average - TMA"
        [t, s, b] = loxxmas.tma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Triple Exponential Moving Average - TEMA"
        [t, s, b] = loxxmas.tema(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Two-pole Ehlers Butterworth"
        [t, s, b] = loxxmas.twopolebutter(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Two-pole Ehlers smoother"
        [t, s, b] = loxxmas.twopoless(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Volume Weighted EMA - VEMA"
        [t, s, b] = loxxmas.vwema(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Zero-Lag DEMA - Zero Lag Double Exponential Moving Average"
        [t, s, b] = loxxmas.zlagdema(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Zero-Lag Moving Average"
        [t, s, b] = loxxmas.zlagma(src, len)
        sig := s
        trig := t
        special := b
    else if type == "Zero Lag TEMA - Zero Lag Triple Exponential Moving Average"
        [t, s, b] = loxxmas.zlagtema(src, len)
        sig := s
        trig := t
        special := b
    trig


int len = 0

if FitMode == "Moving"
    len := Length + MA_Length
else 
    len := Length + 1
        
Polynom = 0.
PolyFitMA = 0.
UpBand = 0.
DnBand = 0.
UpBand1 = 0.
DnBand1 = 0.

PriceArray = array.new<float>(Length, 0)

var Project = array.new_float(ProjLength + 1, 0)
var Projection = array.new_float(ProjLength + 1, 0)
var DnProject = array.new_float(ProjLength + 1, 0)
var UpProject = array.new_float(ProjLength + 1, 0)

var DnProject1 = array.new_float(ProjLength + 1, 0)
var UpProject1 = array.new_float(ProjLength + 1, 0)
poly = 0.

var uplines = array.new_line(0)
var uplines1 = array.new_line(0)

var fvlines = array.new_line(0)
var dnlines = array.new_line(0)
var dnlines1 = array.new_line(0)

if barstate.isfirst
    for i = 0 to 100 - 1
        array.push(uplines, line.new(na, na, na, na))
        array.push(fvlines, line.new(na, na, na, na))
        array.push(dnlines, line.new(na, na, na, na))
        array.push(uplines1, line.new(na, na, na, na))
        array.push(dnlines1, line.new(na, na, na, na))

if last_bar_index - bar_index < barsbark

    for j = 0 to Length  - 1
        array.set(PriceArray, j, nz(variant(type, src, MA_Length)[j]))
    
    Sum = 0.
    for j = Length - 1 to 0
        poly := polyfit(PriceArray, Degree, Length, -j)
        del = array.get(PriceArray, j) - poly
        Sum += del * del
        if FitMode == "Fitting"
            Polynom := poly
    
    if FitMode == "Moving"
        Polynom := poly
    
    StdDev = 0.
    if (Length - 1 > 0 and Sum > 0) 
        if Length < 32
            StdDev := math.sqrt(Sum / (Length - 1))
        else
            StdDev := math.sqrt(Sum / Length)
    
    UpBand := Polynom + K_Sigma * StdDev
    DnBand := Polynom - K_Sigma * StdDev
    
    UpBand1 := Polynom + K_Sigma2 * StdDev
    DnBand1 := Polynom - K_Sigma2 * StdDev
    
    if barstate.islast
        for j = 0 to ProjLength  
            array.set(Project, j, polyfit(PriceArray, Degree, Length, ProjLength - j))
    
        for j = 0 to ProjLength  
            array.set(Projection, j, array.get(Project, j))
            array.set(UpProject, j, array.get(Project, j) + K_Sigma * StdDev)
            array.set(DnProject, j, array.get(Project, j) - K_Sigma * StdDev)
            
            array.set(UpProject1, j, array.get(Project, j) + K_Sigma2 * StdDev)
            array.set(DnProject1, j, array.get(Project, j) - K_Sigma2 * StdDev)

        psize = array.size(Projection)

        skipperfv = psize >= 165 * 3 ? 8 : psize >= 165 * 2 ? 4 : psize >= 165 ? 2 : 1
        
        i = 0
        j = 0
        
        array.reverse(Projection) 
        array.reverse(UpProject) 
        array.reverse(DnProject) 
        array.reverse(UpProject1) 
        array.reverse(DnProject1) 

        outer = math.min(psize, ProjLength)
        while i < outer - skipperfv
            
            if j > array.size(fvlines) - 1
                break
                
            if showmid
                fvline = array.get(fvlines, j)
                line.set_xy1(fvline, bar_index + i + 1, array.get(Projection, i + skipperfv))
                line.set_xy2(fvline, bar_index + i + 1 - skipperfv, array.get(Projection, i))
                line.set_color(fvline, color.white)
                line.set_style(fvline, line.style_dotted)
                line.set_width(fvline, 3)
            
            if showbands
                upline = array.get(uplines, j)
                line.set_xy1(upline, bar_index + i + 1, array.get(UpProject, i + skipperfv))
                line.set_xy2(upline, bar_index + i + 1 - skipperfv, array.get(UpProject, i))
                line.set_color(upline, lightgreencolor)
                line.set_style(upline, line.style_dotted)
                line.set_width(upline, 2)
                
                upline1 = array.get(uplines1, j)
                line.set_xy1(upline1, bar_index + i + 1, array.get(UpProject1, i + skipperfv))
                line.set_xy2(upline1, bar_index + i + 1 - skipperfv, array.get(UpProject1, i))
                line.set_color(upline1, greencolor)
                line.set_style(upline1, line.style_dotted)
                line.set_width(upline1, 2)
                
                dnline = array.get(dnlines, j)
                line.set_xy1(dnline, bar_index + i + 1, array.get(DnProject, i + skipperfv))
                line.set_xy2(dnline, bar_index + i + 1 - skipperfv, array.get(DnProject, i))
                line.set_color(dnline, lightredcolor)
                line.set_style(dnline, line.style_dotted)
                line.set_width(dnline, 2)
    
                dnline1 = array.get(dnlines1, j)
                line.set_xy1(dnline1, bar_index + i + 1, array.get(DnProject1, i + skipperfv))
                line.set_xy2(dnline1, bar_index + i + 1 - skipperfv, array.get(DnProject1, i))
                line.set_color(dnline1, redcolor)
                line.set_style(dnline1, line.style_dotted)
                line.set_width(dnline1, 2)
            
            i += skipperfv
            j += 1
            

colorout = Polynom > Polynom[1] ? greencolor : redcolor

plot(last_bar_index - bar_index < barsbark and showmid ? Polynom : na, "Poly MA", color = colorout, linewidth = 4)

plot(last_bar_index - bar_index < barsbark and showbands ? UpBand : na, "Upper inner band", color = lightgreencolor)
plot(last_bar_index - bar_index < barsbark and showbands ? UpBand1 : na, "Upper outer band", color = greencolor)

plot(last_bar_index - bar_index < barsbark and showbands ? DnBand : na, "Lower inner band", color = lightredcolor)
plot(last_bar_index - bar_index < barsbark and showbands ? DnBand1 : na, "Lower outer band", color = redcolor)

barcolor(last_bar_index - bar_index < barsbark and colorbars ? colorout : na)



