Endpointed SSA 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("Endpointed SSA of Price [Loxx]", overlay = true)

color greencolor = #2DD204
color redcolor = #D2042D 

int Maxncomp = 25
int MaxLag = 200
int MaxArrayLength = 1000

// Calculation of the function Sn, needed to calculate the eigenvalues
// Negative determinants are counted there 
gaussSn(matrix<float> A, float l, int n)=>
    array<float> w = array.new<float>(n, 0)
    matrix<float> B = matrix.copy(A)
    int count = 0
    int cp = 0
    float c = 0.
    float s1 = 0.
    float s2 = 0.
    for i = 0 to n - 1 
        matrix.set(B, i, i, matrix.get(B, i, i) - l)
    for k = 0 to n - 2 
        for i = k + 1 to n - 1 
            if matrix.get(B, k, k) == 0
                for i1 = 0 to n - 1
                    array.set(w, i1, matrix.get(B, i1, k)) 
                    matrix.set(B, i1, k, matrix.get(B, i1, k + 1))
                    matrix.set(B, i1, k + 1, array.get(w, i1))
                cp := cp + 1
            c := matrix.get(B, i, k) / matrix.get(B, k, k) 
            for j = 0 to n - 1 
                matrix.set(B, i, j, matrix.get(B, i, j) - matrix.get(B, k, j) * c)
    count := 0
    s1 := 1
    for i = 0 to n - 1 
        s2 := matrix.get(B, i, i) 
        if s2 < 0
            count := count + 1
    count

// Calculation of eigenvalues by the bisection method}
// The good thing is that as many eigenvalues are needed, so many will count,
// saves a lot of resources
gaussbisectionl(matrix<float> A, int k, int n)=>
    float e1 = 0.
    float maxnorm = 0.
    float cn = 0.
    float a1 = 0.
    float b1 = 0.
    float c = 0.
    for i = 0 to n - 1 
        cn := 0
        for j = 0 to n - 1 
            cn := cn + matrix.get(A, i, i) 
        if maxnorm < cn 
            maxnorm := cn
    a1 := 0
    b1 := 10 * maxnorm
    e1 := 1.0 * maxnorm / 10000000
    while math.abs(b1 - a1) > e1
        c := 1.0 * (a1 + b1) / 2
        if gaussSn(A, c, n) < k
            a1 := c
        else
            b1 := c
    float out = (a1 + b1) / 2.0
    out

// Calculates eigenvectors for already computed eigenvalues 
svector(matrix<float> A, float l, int n, array<float> V)=>
    int cp = 0
    matrix<float> B = matrix.copy(A)
    float c = 0
    array<float> w = array.new<float>(n, 0)
    for i = 0 to n - 1 
        matrix.set(B, i, i, matrix.get(B, i, i) - l)
    for k = 0 to n - 2 
        for i = k + 1 to n - 1
            if matrix.get(B, k, k) == 0
                for i1 = 0 to n - 1 
                    array.set(w, i1, matrix.get(B, i1, k))
                    matrix.set(B, i1, k, matrix.get(B, i1, k + 1))
                    matrix.set(B, i1, k + 1, array.get(w, i1))
                cp += 1
            
            c := 1.0 * matrix.get(B, i, k) / matrix.get(B, k, k) 
            for j = 0 to n - 1 
                matrix.set(B, i, j, matrix.get(B, i, j) - matrix.get(B, k, j) * c)
    array.set(V, n - 1, 1)
    c := 1
    for i = n - 2 to 0 
        array.set(V, i, 0) 
        for j = i to n - 1 
            array.set(V, i, array.get(V, i) - matrix.get(B, i, j) * array.get(V, j))
        array.set(V, i, array.get(V, i) / matrix.get(B, i, i))
        c += math.pow(array.get(V, i), 2)
    for i = 0 to n - 1 
        array.set(V, i, array.get(V, i) / math.sqrt(c))

// Fast Singular SSA - "Caterpillar" method
// X-vector of the original series
// n-length
// l-lag length
// s-number of eigencomponents
// (there the original series is divided into components, and then restored, here you set how many components you need)
// Y - the restored row (smoothed by the caterpillar) 
fastsingular(array<float> X, int n1, int l1, int s1)=>
    int n = math.min(MaxArrayLength, n1)
    int l = math.min(MaxLag, l1)
    int s = math.min(Maxncomp, s1)
    
    matrix<float> A = matrix.new<float>(l, l, 0.)
    matrix<float> B = matrix.new<float>(n, l, 0.)
    matrix<float> Bn = matrix.new<float>(l, n, 0.)
    matrix<float> V = matrix.new<float>(l, n, 0.)
    matrix<float> Yn = matrix.new<float>(l, n, 0.)
    
    var array<float> vtarr = array.new<float>(l, 0.)
    array<float> ls = array.new<float>(MaxLag, 0)
    array<float> Vtemp = array.new<float>(MaxLag, 0)
    array<float> Y = array.new<float>(n, 0)
    
    int k = n - l + 1
    
    // We form matrix A in the method that I downloaded from the site of the creators of this matrix S 
    for i = 0 to l - 1  
        for j = 0 to l - 1 
            matrix.set(A, i, j, 0)
            for m = 0 to k - 1 
                matrix.set(A, i, j, matrix.get(A, i, j) + array.get(X, i + m) * array.get(X, m + j))
                matrix.set(B, m, j, array.get(X, m + j))

    //Find the eigenvalues and vectors of the matrix A
    for i = 0 to s - 1  
        array.set(ls, i, gaussbisectionl(A, l - i, l))
        svector(A, array.get(ls, i), l, Vtemp)  
        for j = 0 to l - 1 
			matrix.set(V, i, j, array.get(Vtemp, j))

    // The restored matrix is formed
    for i1 = 0 to s - 1
        for i = 0 to k - 1  
            matrix.set(Yn, i1, i, 0) 
            for j = 0 to l - 1  
				matrix.set(Yn, i1, i, matrix.get(Yn, i1, i) + matrix.get(B, i, j) * matrix.get(V, i1, j))

        for i = 0 to l - 1 
            for j = 0 to k - 1 
				matrix.set(Bn, i, j, matrix.get(V, i1, i) * matrix.get(Yn, i1, j))
				
        //Diagonal averaging (series recovery)
        int kb = k
        int lb = l
        for i = 0 to n - 1  
            matrix.set(Yn, i1, i, 0)
            if i < lb - 1 
                for j = 0 to i
                    if l <= k
						matrix.set(Yn, i1, i, matrix.get(Yn, i1, i) + matrix.get(Bn, j, i - j))
                    if l > k
						matrix.set(Yn, i1, i, matrix.get(Yn, i1, i) + matrix.get(Bn, i - j, j))
                
                matrix.set(Yn, i1, i, matrix.get(Yn, i1, i) / (1.0 * (i + 1)))
            if (lb - 1 <= i) and (i < kb - 1)
                for j = 0 to lb - 1  
                    if l <= k
						matrix.set(Yn, i1, i, matrix.get(Yn, i1, i) + matrix.get(Bn, j, i - j))
                    if l > k
						matrix.set(Yn, i1, i, matrix.get(Yn, i1, i) + matrix.get(Bn, i - j, j))
                matrix.set(Yn, i1, i, matrix.get(Yn, i1, i) / (1.0 * lb))
            if kb - 1 <= i 
                for j = i - kb + 1 to n - kb  
                    if l <= k
						matrix.set(Yn, i1, i, matrix.get(Yn, i1, i) + matrix.get(Bn, j, i - j))
                    if l > k
						matrix.set(Yn, i1, i, matrix.get(Yn, i1, i) + matrix.get(Bn, i - j, j))
    
                matrix.set(Yn, i1, i, matrix.get(Yn, i1, i) / (1.0 * (n - i))) 

    // Here, if not summarized, then there will be separate decomposition components
    // process by own functions 
    for i = 0 to n - 1
        array.set(Y, i, 0) 
        for i1 = 0 to s - 1  
            array.set(Y, i, array.get(Y, i) + matrix.get(Yn, i1, i))
    Y

float src = input.source(close, "Source", group = "Basic Settings")
int lag = input.int(10, "Lag", group = "Basic Settings")
int ncomp = input.int(2, "Number of Computations", group = "Basic Settings")
int ssapernorm = input.int(30, "SSA Period Normalization", group = "Basic Settings")
int numbars = input.int(330, "Number of Bars", group = "Basic Settings")
int backbars = input.int(400, "Number of Bars to Render", group = "Basic Settings", tooltip ="How many bars to plot.The higher the number the slower the computation.")

bool colorbars = input.bool(true, "Color bars?", group = "UI Options")
bool showSigs = input.bool(false, "Show signals?", group = "UI Options")

float out = 0.
float sig = 0.
float[] srcVal = array.new_float(numbars + 1, 0)

if last_bar_index - bar_index < backbars
    for i = 0 to numbars - 1
        array.set(srcVal, i, nz(src[i]))

    float[] pv = fastsingular(srcVal, numbars, lag, ncomp) 
    out := array.get(pv, 0)
    sig := out[1]

color colorout = out > sig ? greencolor : redcolor

bool goLong = ta.crossover(out, sig)
bool goShort = ta.crossunder(out, sig)

plot(last_bar_index - bar_index < backbars ? out : na, "Endpointed SSA of Price", color = colorout, linewidth = 2)
barcolor(last_bar_index - bar_index < backbars and colorbars ? colorout : na)

plotshape(showSigs and goLong, title = "Long", color = color.yellow, textcolor = color.yellow, text = "L", style = shape.triangleup, location = location.belowbar, size = size.auto)
plotshape(showSigs and goShort, title = "Short", color = color.fuchsia, textcolor = color.fuchsia, text = "S", style = shape.triangledown, location = location.abovebar, size = size.auto)

alertcondition(goLong, title = "Long", message = "Endpointed SSA of Price [Loxx]: Long\nSymbol: {{ticker}}\nPrice: {{close}}")
alertcondition(goShort, title = "Short", message = "Endpointed SSA of Price [Loxx]: Short\nSymbol: {{ticker}}\nPrice: {{close}}")