// Licensed under the Apache License, Version 2.0
// © mihakralj
//@version=6
indicator("Cubic Regression Moving Average (CRMA)", "CRMA", overlay=true)

//@function Computes Cubic Regression Moving Average — fits a degree-3 polynomial
//          y = a0 + a1*x + a2*x² + a3*x³ to the most recent `period` bars via
//          normal equations solved by closed-form elimination, returns the fitted endpoint.
//@param source Series to analyze
//@param period Lookback window for the cubic regression
//@returns Fitted value at the most recent bar (x = 0)
//@reference Polynomial least-squares regression (degree 3), evaluated at endpoint
//@optimized O(period) per bar for accumulating sums; O(1) for solve
crma(series float source, simple int period) =>
    if period < 4
        runtime.error("Period must be at least 4 for cubic regression")

    float price = nz(source)

    // --- Circular buffer for rolling window ---
    var array<float> buffer = array.new_float(period, na)
    var int head = 0

    array.set(buffer, head, price)
    head := (head + 1) % period

    int p = math.min(bar_index + 1, period)
    if p < 4
        price
    else
        // --- Accumulate power sums and cross-products ---
        float s0 = 0.0
        float s1 = 0.0
        float s2 = 0.0
        float s3 = 0.0
        float s4 = 0.0
        float s5 = 0.0
        float s6 = 0.0
        float r0 = 0.0
        float r1 = 0.0
        float r2 = 0.0
        float r3 = 0.0

        int idx = (head - 1 + period) % period
        for i = 0 to p - 1
            float val = array.get(buffer, idx)
            float v = na(val) ? price : val
            float x  = float(i)
            float x2 = x * x
            float x3 = x2 * x

            s0 += 1.0
            s1 += x
            s2 += x2
            s3 += x3
            s4 += x2 * x2
            s5 += x2 * x3
            s6 += x3 * x3

            r0 += v
            r1 += x * v
            r2 += x2 * v
            r3 += x3 * v

            idx := (idx - 1 + period) % period

        // --- Solve 4x4 via scalar Gaussian elimination with partial pivoting ---
        // Row 0
        float a00 = s0
        float a01 = s1
        float a02 = s2
        float a03 = s3
        float b0  = r0
        // Row 1
        float a10 = s1
        float a11 = s2
        float a12 = s3
        float a13 = s4
        float b1  = r1
        // Row 2
        float a20 = s2
        float a21 = s3
        float a22 = s4
        float a23 = s5
        float b2  = r2
        // Row 3
        float a30 = s3
        float a31 = s4
        float a32 = s5
        float a33 = s6
        float b3  = r3

        // --- Column 0: partial pivoting among rows 0-3 ---
        float best0 = math.abs(a00)
        int piv0 = 0
        if math.abs(a10) > best0
            best0 := math.abs(a10)
            piv0 := 1
        if math.abs(a20) > best0
            best0 := math.abs(a20)
            piv0 := 2
        if math.abs(a30) > best0
            best0 := math.abs(a30)
            piv0 := 3

        // Swap pivot row with row 0
        if piv0 == 1
            float t0 = a00
            float t1 = a01
            float t2 = a02
            float t3 = a03
            float t4 = b0
            a00 := a10
            a01 := a11
            a02 := a12
            a03 := a13
            b0 := b1
            a10 := t0
            a11 := t1
            a12 := t2
            a13 := t3
            b1 := t4
        else if piv0 == 2
            float t0 = a00
            float t1 = a01
            float t2 = a02
            float t3 = a03
            float t4 = b0
            a00 := a20
            a01 := a21
            a02 := a22
            a03 := a23
            b0 := b2
            a20 := t0
            a21 := t1
            a22 := t2
            a23 := t3
            b2 := t4
        else if piv0 == 3
            float t0 = a00
            float t1 = a01
            float t2 = a02
            float t3 = a03
            float t4 = b0
            a00 := a30
            a01 := a31
            a02 := a32
            a03 := a33
            b0 := b3
            a30 := t0
            a31 := t1
            a32 := t2
            a33 := t3
            b3 := t4

        bool singular = best0 < 1e-12
        if not singular
            // Eliminate column 0 from rows 1,2,3
            float f1 = a10 / a00
            a11 -= f1 * a01
            a12 -= f1 * a02
            a13 -= f1 * a03
            b1  -= f1 * b0
            a10 := 0.0

            float f2 = a20 / a00
            a21 -= f2 * a01
            a22 -= f2 * a02
            a23 -= f2 * a03
            b2  -= f2 * b0
            a20 := 0.0

            float f3 = a30 / a00
            a31 -= f3 * a01
            a32 -= f3 * a02
            a33 -= f3 * a03
            b3  -= f3 * b0
            a30 := 0.0

            // --- Column 1: partial pivoting among rows 1-3 ---
            float best1 = math.abs(a11)
            int piv1 = 1
            if math.abs(a21) > best1
                best1 := math.abs(a21)
                piv1 := 2
            if math.abs(a31) > best1
                best1 := math.abs(a31)
                piv1 := 3

            if piv1 == 2
                float t1 = a11
                float t2 = a12
                float t3 = a13
                float t4 = b1
                a11 := a21
                a12 := a22
                a13 := a23
                b1 := b2
                a21 := t1
                a22 := t2
                a23 := t3
                b2 := t4
            else if piv1 == 3
                float t1 = a11
                float t2 = a12
                float t3 = a13
                float t4 = b1
                a11 := a31
                a12 := a32
                a13 := a33
                b1 := b3
                a31 := t1
                a32 := t2
                a33 := t3
                b3 := t4

            singular := best1 < 1e-12

        if not singular
            // Eliminate column 1 from rows 2,3
            float g2 = a21 / a11
            a22 -= g2 * a12
            a23 -= g2 * a13
            b2  -= g2 * b1
            a21 := 0.0

            float g3 = a31 / a11
            a32 -= g3 * a12
            a33 -= g3 * a13
            b3  -= g3 * b1
            a31 := 0.0

            // --- Column 2: partial pivoting among rows 2-3 ---
            if math.abs(a32) > math.abs(a22)
                float t2 = a22
                float t3 = a23
                float t4 = b2
                a22 := a32
                a23 := a33
                b2 := b3
                a32 := t2
                a33 := t3
                b3 := t4

            singular := math.abs(a22) < 1e-12

        if not singular
            // Eliminate column 2 from row 3
            float h3 = a32 / a22
            a33 -= h3 * a23
            b3  -= h3 * b2
            a32 := 0.0

            singular := math.abs(a33) < 1e-12

        float result = price
        if not singular
            // Back-substitution
            float c3 = b3 / a33
            float c2 = (b2 - a23 * c3) / a22
            float c1 = (b1 - a12 * c2 - a13 * c3) / a11
            float c0 = (b0 - a01 * c1 - a02 * c2 - a03 * c3) / a00
            result := c0

        result

// ---------- Main loop ----------

// Inputs
i_period = input.int(14, "Period", minval=4)
i_source = input.source(close, "Source")

// Calculation
crma_value = crma(i_source, i_period)

// Plot
plot(crma_value, "CRMA", color=color.yellow, linewidth=2)