Modified Covariance Autoregressive Estimator 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("Modified Covariance Autoregressive Estimator of Price [Loxx]", 
     shorttitle = "MCAEP [Loxx]",
     overlay = true, 
     max_lines_count = 500)

import loxx/loxxexpandedsourcetypes/4

greencolor = #2DD204
redcolor = #D2042D 

darkGreenColor =  #1B7E02 
darkRedColor = #93021F

//+-------------------------------------------------------------------------+
// Fast algorithm for the solution of the modified covariance least squares normal equations.
//
//    This implementation is based on Modified Covariance method from Marple's book
//
//    x :    - Array of complex data samples X(1) through X(N)
//    ip:    - Order of linear prediction model (integer)
//
//    variables
//    v   - Real linear prediction variance at order IP
//    a    - Array of complex linear prediction coefficients
//
//    Bool stop value at time of exit
//      false - for normal exit (no numerical ill-conditioning)
//      true - if v is not a positive value
//      true - if delta and gamma do not lie in the range 0 to 1
//      true - if v is not a positive value
//      true - if delta and gamma do not lie in the range 0 to 1
//    Error message
//      see above for bool stop
//+-------------------------------------------------------------------------+
_MCov(float[] x, int ip)=>
    string errormessage = ""
    bool stop = false
    int n = array.size(x)
    float[] c = array.new_float(ip + 1, 0.)
    float[] d = array.new_float(ip + 1, 0.)
    float[] r = array.new_float(ip, 0.)
    float[] a = array.new<float>(n, 0.)

    float v = 0.
    int mk = 0
    float r1 = 0.
    float r2 = 0.
    float r3 = 0.
    float r4 = 0.
    float r5 = 0.
    float delta = 0.
    float gamma = 0.
    float lambda = 0.
    float theta = 0.
    float psi = 0.
    float xi = 0.
    float save1 = 0.
    float save2 = 0.
    float save3 = 0.
    float save4 = 0.
    float c1 = 0.
    float c2 = 0.
    float c3 = 0.
    float c4 = 0.
    float ef = 0.
    float eb = 0.
    
    r1 := 0.0
    for k = 1 to n - 2
        r1 += 2.0 * math.pow(array.get(x, k), 2)
    r2 := math.pow(array.get(x, n - 1), 2)
    r3 := math.pow(array.get(x, 0), 2)
    r4 := 1.0 / (r1 + 2.0 * (r2 + r3))
    v := r1 + r2 + r3
    delta := 1.0 - r2 * r4
    gamma := 1.0 - r3 * r4
    lambda := array.get(x, 0) * array.get(x, n - 1) * r4
    array.set(c, 0, array.get(x, 0) * r4)
    array.set(c, 0, array.get(x, n - 1) * r4)

    int m = 0
    while true 
        save1 := 0.0
        for k = m + 1 to n - 1
            save1 += array.get(x, n - 1 - k) * array.get(x, n - k + m)
        save1 *= 2.0
        array.set(r, m, save1)
        theta := array.get(x, 0) * array.get(d, 0)
        psi := array.get(x, 0) * array.get(c, 0)
        xi := array.get(x, n - 1) * array.get(d, 0)
        if m > 0
            for k = 1 to m 
                theta += array.get(x, k) * array.get(d, k)
                psi += array.get(x, k) * array.get(c, k)
                xi += array.get(x, n - 1 - k) * array.get(d, k)
                array.set(r, k - 1, array.get(r, k - 1) - array.get(x, m) * array.get(x, m - k) + array.get(x, n - 1 - m) * array.get(x, n - 1 - m + k))
                save1 += array.get(r, k - 1) * array.get(a, m - k)

        //Order update of a vector
        c1 := -save1 / v
        array.set(a, m, c1)
        v *= (1.0 - c1 * c1)

        if m > 0
            for k = 0 to (m + 1) / 2  - 1 
                mk := m - k - 1
                save1 := array.get(a, k)
                array.set(a, k, save1 + c1 * array.get(a, mk))
                if k != mk
                    array.set(a, mk, array.get(a, mk) + c1 * save1)

        if m == ip - 1
            v *= (0.5 / (n - 1 - m))
            break
            
        //Time update of c,d vectors and gamma, delta, lambda scalars
        r1 := 1.0 / (delta * gamma - lambda * lambda)
        c1 := (theta * lambda + psi * delta) * r1
        c2 := (psi * lambda + theta * gamma) * r1
        c3 := (xi * lambda + theta * delta) * r1
        c4 := (theta * lambda + xi * gamma) * r1
        for k = 0 to m / 2 
            mk := m - k
            save1 := array.get(c, k)
            save2 := array.get(d, k) 
            save3 := array.get(c, mk)  
            save4 := array.get(d, mk) 
            array.set(c, k, array.get(c, k) + (c1 * save3 + c2 * save4))
            array.set(d, k, array.get(d, k) + (c3 * save3 + c4 * save4))
            if (k != mk) 
                array.set(c, k, array.get(c, k) + (c1 * save1 + c2 * save2))
                array.set(d, k, array.get(d, k) + (c3 * save1 + c4 * save2))
            
        r2 := psi * psi
        r3 := theta * theta
        r4 := xi * xi
        r5 := gamma - (r2 * delta + r3 * gamma + 2.0 * psi * lambda * theta) * r1
        r2 := delta - (r3 * delta + r4 * gamma + 2. * theta * lambda * xi) * r1
        gamma := r5
        delta := r2
        lambda += (c3 * psi + c4 * theta)
        if (v <= 0.0) 
            errormessage := "Error: Negative or zero value of the v variable"
            stop := true
            break
        
        if (delta <= 0.0 or delta > 1.0 or gamma <= 0.0 or gamma > 1.0) 

            errormessage := "Error: delta and gamma variables values out of the range (0,1)"
            stop := true
            break
        
        r1 := 1.0 / v
        r2 := 1.0 / (delta * gamma - lambda * lambda)
        ef := array.get(x, n - m - 2)
        eb := array.get(x, m + 1)
        for k = 0 to m  
            ef += array.get(a, k) * array.get(x, n - 1 - m + k)
            eb += array.get(a, k) * array.get(x, m - k)
        
        c1 := eb * r1
        c2 := ef * r1
        c3 := (eb * delta + ef * lambda) * r2
        c4 := (ef * gamma + eb * lambda) * r2
        for k = m to 0  
            save1 := array.get(a, k)
            array.set(a, k, save1 + c3 * array.get(c, k) + c4 * array.get(d, k))
            array.set(c, k + 1, array.get(c, k) + c1 * save1)
            array.set(d, k + 1, array.get(d, k) + c2 * save1)
        
        array.set(c, 0, c1)
        array.set(d, 0, c2)

        r3 := eb * eb
        r4 := ef * ef
        v -= (r3 * delta + r4 * gamma + 2.0 * ef * eb * lambda) * r2
        delta -= r4 * r1
        gamma -= r3 * r1
        lambda += ef * eb * r1
        if (v <= 0.0) 
            errormessage := "Error: Negative or zero value of the v variable"
            stop := true
            break
        
        if (delta <= 0.0 or delta > 1.0 or gamma <= 0.0 or gamma > 1.0) 

            errormessage := "Error: delta and gamma variables values out of the range (0,1)"
            stop := true
            break
        m += 1
    [a, stop, errormessage]
    

smthtype = input.string("Kaufman", "Heikin-Ashi Better Caculation Type", options = ["AMA", "T3", "Kaufman"], group = "Source Settings")
srcin = input.string("Open", "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)"])

LastBar = input.int(30, "Last Bar", group = "Basic Settings", tooltip = "Bar from where to start prediction")
PastBars = input.int(500, "Past Bars", group = "Basic Settings", maxval = 2000)

LPOrder = input.float(0.1, "Order of Linear Prediction", group = "Basic Settings", minval = 0, maxval = 0.67, step = 0.01, 
     tooltip = "Order must be less than or equal to 2/3 the input frame size")
     
FutBars = input.int(20, "Future Bars", group = "Basic Settings", maxval = 500)

colorbars = input.bool(true, "Mute bar colors?", group = "UI Options")

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
	
barcolor(colorbars ? color.gray : na)

int lb = LastBar
int np = PastBars
int no = math.ceil(LPOrder * PastBars)
int nf = np - no - 1

float[] x = array.new<float>(np, 0.)
float[] pv = array.new<float>(np, 0.)
float[] fv = array.new<float>(nf + 1, 0.)

//Find average of past values

var pvlines = array.new_line(0)
var fvlines = array.new_line(0)

cnp = np >= 250 ? 250 : np
cnf = nf >= 250 ? 250 : nf

if barstate.isfirst
    for i = 0 to 250 - 1
        array.push(pvlines, line.new(na, na, na, na))
        array.push(fvlines, line.new(na, na, na, na))

if barstate.islast
    //Prepare data
    
    float av = 0.
    avar = array.new<float>(np, 0.)
    for i = 0 to np - 1
        array.set(avar, i, nz(src[i + lb]))
    av := array.avg(avar)
    
    for i = 0 to np - 1
        array.set(x, np - 1 - i, nz(src[i + lb]) - av)
    
    //Use linear prediction _Geom
    [result, stop, errormessage] = _MCov(x, no)
    
    //grab test results
    strout = stop ? "🚨 " + errormessage : "✅ Passed"
    colort = stop ? darkRedColor : darkGreenColor
    var testTable = table.new(position = position.bottom_right, columns = 2, rows = 1, bgcolor = colort, border_width = 1)
    table.cell(table_id = testTable, column = 0, row = 0, text = strout, text_color = color.white, text_size = size.large)

    //correct data
    for i = no to 1 
        array.set(result, i, array.get(result, i - 1))
            
    //Calculate linear predictions

    for n = no to np + nf - 1
        float sum = 0.
        for i = 1 to no 
            if (n - i < np)
                sum -= array.get(result, i) * array.get(x, n - i)
            else
                sum -= array.get(result, i) * array.get(fv, n - i - np + 1)
        if (n < np)
            array.set(pv, np - 1 - n, sum)
        else
            array.set(fv, n - np + 1, sum)
    
    array.set(fv, 0, array.get(pv, 0))
    for i = 0 to np - no - 1
        array.set(pv, i, array.get(pv, i) + av)
        array.set(fv, i, array.get(fv, i) + av)
        
    //+------------------------------------------------------------------+
    //| Draw lines w/ skipping to stay within 500 line limit
    //+------------------------------------------------------------------+
    skipperpv = array.size(pv) >= 2000 ? 8 : array.size(pv) >= 1000 ? 4 : array.size(pv) >= 500 ? 2 : 1
    int i = 0
    int j = 0

    while i < np - no - 1 - skipperpv

        if j > array.size(pvlines) - 1
            break

        pvline = array.get(pvlines, j)
        line.set_xy1(pvline, bar_index - i - skipperpv - LastBar, array.get(pv, i + skipperpv))
        line.set_xy2(pvline, bar_index - i - LastBar, array.get(pv, i))
        line.set_color(pvline, greencolor)
        line.set_style(pvline, line.style_solid)
        line.set_width(pvline, 3)

        i += skipperpv
        j += 1
    
    skipperfv = array.size(fv) >= 2000 ? 8 : array.size(fv) >= 1000 ? 4 : array.size(fv) >= 500 ? 2 : 1
    
    i := 0
    j := 0
    outer = math.min(np - no - 1, FutBars)
    while i < outer - skipperfv
        
        if j > array.size(fvlines) - 1
            break
            
        fvline = array.get(fvlines, j)
        line.set_xy1(fvline, bar_index + i + 1 - LastBar, array.get(fv, i + skipperfv))
        line.set_xy2(fvline, bar_index + i + 1 - LastBar - skipperfv, array.get(fv, i))
        line.set_color(fvline, color.blue)
        line.set_style(fvline, line.style_solid)
        line.set_width(fvline, 2)
        
        i += skipperfv
        j += 1


