// The MIT License (MIT)
// © 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 with Gaussian 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 4×4 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 ---
        // Normal equations for degree-3 polynomial: M * a = rhs
        // M[i][j] = Σ x^(i+j),  rhs[i] = Σ x^i * y   for i,j = 0..3
        // x = 0 (newest) to p-1 (oldest), so a0 = fitted value at newest bar
        float s0 = 0.0,  s1 = 0.0,  s2 = 0.0,  s3 = 0.0
        float s4 = 0.0,  s5 = 0.0,  s6 = 0.0
        float r0 = 0.0,  r1 = 0.0,  r2 = 0.0,  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       // Σ x^0 = count
            s1 += x         // Σ x^1
            s2 += x2        // Σ x^2
            s3 += x3        // Σ x^3
            s4 += x2 * x2   // Σ x^4
            s5 += x2 * x3   // Σ x^5
            s6 += x3 * x3   // Σ x^6

            r0 += v          // Σ y
            r1 += x * v      // Σ x*y
            r2 += x2 * v     // Σ x²*y
            r3 += x3 * v     // Σ x³*y

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

        // --- Build 4×4 augmented matrix (row-major, 4 rows × 5 cols) ---
        var matrix<float> m = matrix.new<float>(4, 5, 0.0)

        // Row 0: [s0, s1, s2, s3 | r0]
        matrix.set(m, 0, 0, s0), matrix.set(m, 0, 1, s1), matrix.set(m, 0, 2, s2), matrix.set(m, 0, 3, s3), matrix.set(m, 0, 4, r0)
        // Row 1: [s1, s2, s3, s4 | r1]
        matrix.set(m, 1, 0, s1), matrix.set(m, 1, 1, s2), matrix.set(m, 1, 2, s3), matrix.set(m, 1, 3, s4), matrix.set(m, 1, 4, r1)
        // Row 2: [s2, s3, s4, s5 | r2]
        matrix.set(m, 2, 0, s2), matrix.set(m, 2, 1, s3), matrix.set(m, 2, 2, s4), matrix.set(m, 2, 3, s5), matrix.set(m, 2, 4, r2)
        // Row 3: [s3, s4, s5, s6 | r3]
        matrix.set(m, 3, 0, s3), matrix.set(m, 3, 1, s4), matrix.set(m, 3, 2, s5), matrix.set(m, 3, 3, s6), matrix.set(m, 3, 4, r3)

        // --- Gaussian elimination with partial pivoting ---
        bool singular = false
        for col = 0 to 3
            // Find pivot row
            int pivot_row = col
            float pivot_max = math.abs(matrix.get(m, col, col))
            for row = col + 1 to 3
                float absval = math.abs(matrix.get(m, row, col))
                if absval > pivot_max
                    pivot_max := absval
                    pivot_row := row

            if pivot_max < 1e-12
                singular := true
                break

            // Swap rows if needed
            if pivot_row != col
                for k = col to 4
                    float tmp = matrix.get(m, col, k)
                    matrix.set(m, col, k, matrix.get(m, pivot_row, k))
                    matrix.set(m, pivot_row, k, tmp)

            // Eliminate below
            float diag = matrix.get(m, col, col)
            for row = col + 1 to 3
                float factor = matrix.get(m, row, col) / diag
                for k = col to 4
                    matrix.set(m, row, k, matrix.get(m, row, k) - factor * matrix.get(m, col, k))

        float result = price
        if not singular
            // Back-substitution
            var array<float> a = array.new_float(4, 0.0)
            for row = 3 to 0
                float val = matrix.get(m, row, 4)
                for k = row + 1 to 3
                    val -= matrix.get(m, row, k) * array.get(a, k)
                array.set(a, row, val / matrix.get(m, row, row))

            // a[0] is the fitted value at x = 0 (most recent bar)
            result := array.get(a, 0)

        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)