Local Polynomial Regression
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © wzkkzw12345

//@version=5
indicator("Kernel Regression",overlay=true,max_bars_back=1000,max_lines_count=500,max_labels_count=500)
length = input.int(500,'Look Back Length',maxval=500,minval=0)
kernel_size = input.int(45,'Kernel Size',maxval=150,minval=0)
h      = input.float(8.5,'Bandwidth') 
src    = (close*2+high+low)/4//input.source(hlc3,'Estimation Source')
order = input.int(4, 'Order of Regression',maxval=10,minval=0)
x_scaling = input.int(10, 'Custom Scaling for Scalability (to compute (x/scaling)^i)',maxval=100,minval=1)
use_volume = input.bool(false, "To use volume as an additional weight (currently it seems impossible...)")
repaint = input.bool(true, "Repaint the curve when prices get known")


//var log_src = 0.
log_src = request.security(syminfo.tickerid, timeframe.period, math.log(src))

//max_bars_back(log_src, 1000)
up_col = input.color(#39ff14,'Colors',inline='col')
dn_col = input.color(#ff1100,'',inline='col')

//----
n = bar_index
var k = 2
var kernel = array.new_float(0) 
var ln = array.new_line(0) 
var labels = array.new_label(0)
var weight_matrix = array.new_float((kernel_size+1)*(order+1),0) 

  
lset(l,x1,y1,x2,y2,col)=>
    line.set_xy1(l,x1,y1)
    line.set_xy2(l,x2,y2)
    line.set_color(l,col)
    line.set_width(l,2)
labelset(l,x,y,text_,col,style,textcolor,textalign)=>
    label.set_xy(l,x,y)
    label.set_text(l,text_)
    label.set_color(l,col)
    label.set_style(l,style)
    label.set_textcolor(l,textcolor)
    label.set_textalign(l,textalign)

// ------------------------------------------------------------------------------------
//    Helper functions for simulating global variables using arrays
// ------------------------------------------------------------------------------------

new_string(value) =>
//{
    array.new_string(1, value)
//}

new_color(value) =>
//{
    array.new_color(1, value)
//}

get(global_variable_id) =>
//{
    if (array.size(global_variable_id) > 0)
        array.get(global_variable_id, 0)
//}

set(global_variable_id, value) =>
//{
    if (array.size(global_variable_id) > 0)
        array.set(global_variable_id, 0, value)
//}

//} GlobalVariables



// ====================================================================================================================================================
// 
//       LIBRARY: OutputStreams
// 
// ====================================================================================================================================================
//{

// ------------------------------------------------------------------------------------
//    Standard output
// ------------------------------------------------------------------------------------

var color NO_COLOR = color.new(#000000, 100)

var string[] output_stream = new_string("")
var color[] output_color = new_color(color.white)
var label output_label = label.new(0, 0, style = label.style_label_left, color = NO_COLOR, textcolor = array.get(output_color, 0))

clear_output() =>
//{
    array.clear(output_stream)
    array.push(output_stream, "")
    label.set_text(output_label, get(output_stream))
//}

set_output_visible(visible) =>
//{
    if (not visible)
    //{
        label.set_style(output_label, label.style_none)
        label.set_textcolor(output_label, NO_COLOR)
    //}
    else
    //{
        label.set_style(output_label, label.style_label_left)
        label.set_textcolor(output_label, get(output_color))
    //}
//}

set_output_color(clr) =>
//{
    set(output_color, clr)
//}

set_output_location(max_y, min_y) =>
//{
    label.set_y(output_label, (max_y + min_y) / 2)
//}

output(msg) =>
//{
    set(output_stream, get(output_stream) + msg)
    
    label.set_x(output_label, bar_index)
    label.set_text(output_label, get(output_stream))
//}

reset_output() =>
//{
    clear_output()
    set_output_color(color.white)
    set_output_location(0, 0)
//}

//} GlobalVariables
// ====================================================================================================================================================
// 
//       LIBRARY: ArrayLib
// 
// ====================================================================================================================================================
//{

// ------------------------------------------------------------------------------------
//    Create an array from a time series
// ------------------------------------------------------------------------------------

array_create_from_series(series, lookback) =>
//{
    float[] new_array = array.copy(array.new_float(lookback, 0))
    
    for i = 0 to lookback - 1
    //{
        if (i >= lookback)
            break
        
        array.set(new_array, i, nz(series[i]))
    //}
    
    new_array
//}



// ------------------------------------------------------------------------------------
//    Euclidean 2-norm of an array
// ------------------------------------------------------------------------------------

array_2_norm(array) =>
//{
    int size = array.size(array)
    float result = 0
    
    for i = 0 to size - 1
    //{
        if (i >= size)
            break
        
        result := result + math.pow(array.get(array, i), 2)
    //}
    
    math.sqrt(result)
//}



// ------------------------------------------------------------------------------------
//    Dot product of two arrays (vectors)
// ------------------------------------------------------------------------------------

array_dot_product(a_vector, b_vector) =>
//{
    int size_a = array.size(a_vector)
    int size_b = array.size(b_vector)
    
    float result = 0
    
    if (size_b >= size_a)
    //{
        for i = 0 to size_a - 1
        //{
            if (i >= size_a)
                break
            
            result := result + array.get(a_vector, i) * array.get(b_vector, i)
        //}
    //}
    
    result
//}

//} ArrayLib



// ====================================================================================================================================================
// 
//       LIBRARY: MatrixLib
// 
// ====================================================================================================================================================
//{

// ------------------------------------------------------------------------------------
//    Constants
// ------------------------------------------------------------------------------------
//{

var bool MM_END_IF = false
var float MM_RANK_CUTOFF = 10E-12

var int MM_META_SIZE = 2

var int MM_ROWS_INDEX = 0
var int MM_COLUMNS_INDEX = 1

var int MM_OP_ADD = 100
var int MM_OP_SUBTRACT = 101
var int MM_OP_MULTIPLY = 102
var int MM_OP_DIVIDE = 103
var int MM_OP_MODULUS = 104

var int MM_OP_SIN = 200
var int MM_OP_ASIN = 201
var int MM_OP_COS = 202
var int MM_OP_ACOS = 203
var int MM_OP_TAN = 204
var int MM_OP_ATAN = 205
var int MM_OP_TANH = 206
var int MM_OP_ATANH = 207
var int MM_OP_LOG = 208
var int MM_OP_LOG10 = 209
var int MM_OP_SIGMOID = 210
var int MM_OP_SILU = 211
var int MM_OP_GELU = 212
var int MM_OP_RELU = 213
var int MM_OP_LEAKY_RELU = 214
var int MM_OP_SIGN = 215

var int MM_MLR_ESTIMATE = 0
var int MM_MLR_SE_ESTIMATE = 1
var int MM_MLR_F_VALUE = 2
var int MM_MLR_CI = 3
var int MM_MLR_P_VALUE = 4
var int MM_MLR_R2 = 5
var int MM_MLR_ADJ_R2 = 6
var int MM_MLR_SSR = 7
var int MM_MLR_SSE = 8
var int MM_MLR_SST = 9
var int MM_MLR_MSR = 10
var int MM_MLR_MSE = 11
var int MM_MLR_MST = 12
var int MM_MLR_DFR = 13
var int MM_MLR_DFE = 14
var int MM_MLR_DFT = 15
var int MM_MLR_NUM_OF_OBSERVATIONS = 16
var int MM_MLR_NUM_OF_X_VARS = 17
var int MM_MLR_HAS_CONSTANT = 18

//} Constants



// ------------------------------------------------------------------------------------
//    Functions
// ------------------------------------------------------------------------------------
//{

// ------------------------------------------------------------------------------------
//    Create and return a copy of the specified matrix
// ------------------------------------------------------------------------------------

matrix_copy(matrix_array) =>
//{
    array.copy(matrix_array)
//}



// ------------------------------------------------------------------------------------
//    Create and return a new matrix with the same shape as the specified matrix, but
//    filled with 0's
// ------------------------------------------------------------------------------------

matrix_copy_shape(matrix_array) =>
//{
    int size = array.size(matrix_array)
    float[] result = array.copy(matrix_array)
    
    if (size > MM_META_SIZE)
    //{
        for i = MM_META_SIZE to size - 1
            array.set(result, i, 0)
    //}
    
    result
//}



// ------------------------------------------------------------------------------------
//    Clear a matrix, removing all rows and columns, resulting in an empty matrix
// ------------------------------------------------------------------------------------

matrix_clear(matrix_array) =>
//{
    array.clear(matrix_array)
    
    for i = 1 to MM_META_SIZE
        array.push(matrix_array, 0)
//}



// ------------------------------------------------------------------------------------
//    Create an empty matrix
// ------------------------------------------------------------------------------------

matrix_create_empty() =>
//{
    array.copy(array.new_float(MM_META_SIZE, 0))
//}



// ------------------------------------------------------------------------------------
//    Create a 'rows' by 'columns' matrix
// ------------------------------------------------------------------------------------

matrix_create(rows, columns) =>
//{
    float[] new_matrix = matrix_create_empty()
    
    if (rows > 0 and columns > 0)
    //{
        for i = 0 to rows * columns - 1
            array.push(new_matrix, 0)
        
        array.set(new_matrix, MM_ROWS_INDEX, rows)
        array.set(new_matrix, MM_COLUMNS_INDEX, columns)
    //}
    
    new_matrix
//}



// ------------------------------------------------------------------------------------
//    Create an identity matrix
// ------------------------------------------------------------------------------------

matrix_create_identity(size) =>
//{
    float[] new_matrix = matrix_create_empty()
    
    if (size > 0)
    //{
        for column = 0 to size - 1
            for row = 0 to size - 1
                array.push(new_matrix, column == row ? 1 : 0)
        
        array.set(new_matrix, MM_ROWS_INDEX, size)
        array.set(new_matrix, MM_COLUMNS_INDEX, size)
    //}
    
    new_matrix
//}



// ------------------------------------------------------------------------------------
//    Get the column count of a matrix
// ------------------------------------------------------------------------------------

matrix_get_column_count(matrix_array) =>
//{
    int(array.get(matrix_array, MM_COLUMNS_INDEX))
//}



// ------------------------------------------------------------------------------------
//    Get the row count of a matrix
// ------------------------------------------------------------------------------------

matrix_get_row_count(matrix_array) =>
//{
    int(array.get(matrix_array, MM_ROWS_INDEX))
//}



// ------------------------------------------------------------------------------------
//    Get the underlying array index of the start of the specified column
// ------------------------------------------------------------------------------------

matrix_get_column_index(matrix_array, column) =>
//{
    int(MM_META_SIZE + matrix_get_row_count(matrix_array) * column)
//}



// ------------------------------------------------------------------------------------
//    Get matrix order as a string in the format "ROWSxCOLUMNS"
// ------------------------------------------------------------------------------------

matrix_get_order(matrix_array) =>
//{
    str.tostring(matrix_get_row_count(matrix_array), "0") + "x" + str.tostring(matrix_get_column_count(matrix_array), "0")
//}



// ------------------------------------------------------------------------------------
//    Get the size of the matrix (row count multiplied by column count)
// ------------------------------------------------------------------------------------

matrix_get_size(matrix_array) =>
//{
    matrix_get_row_count(matrix_array) * matrix_get_column_count(matrix_array)
//}



// ------------------------------------------------------------------------------------
//    Returns true if the matrix is square (same number of rows as columns)
// ------------------------------------------------------------------------------------

matrix_is_square(matrix_array) =>
//{
    matrix_get_row_count(matrix_array) == matrix_get_column_count(matrix_array)
//}

// ------------------------------------------------------------------------------------
//    Get the value from a specific position in a matrix
// ------------------------------------------------------------------------------------

matrix_get(matrix_array, row, column) =>
//{
    int row_count = matrix_get_row_count(matrix_array)
    int column_count = matrix_get_column_count(matrix_array)
    int value_index = matrix_get_column_index(matrix_array, column) + row
    
    if (value_index <= array.size(matrix_array) - 1 and row < row_count and column < column_count and row >= 0 and column >= 0)
        array.get(matrix_array, value_index)
//}



// ------------------------------------------------------------------------------------
//    Set the value at a specific position in a matrix
// ------------------------------------------------------------------------------------

matrix_set(matrix_array, row, column, value) =>
//{
    int row_count = matrix_get_row_count(matrix_array)
    int column_count = matrix_get_column_count(matrix_array)
    int value_index = matrix_get_column_index(matrix_array, column) + row
    
    if (value_index <= array.size(matrix_array) - 1 and row < row_count and column < column_count and row >= 0 and column >= 0)
        array.set(matrix_array, value_index, value)
//}
// ------------------------------------------------------------------------------------
//    Transpose a matrix so that each row becomes a column and return the new matrix
// ------------------------------------------------------------------------------------

matrix_transpose(matrix_array) =>
//{
    int row_count = matrix_get_row_count(matrix_array)
    int column_count = matrix_get_column_count(matrix_array)
    
    // Create the matrix that will store the result of the operation
    float[] result_matrix = matrix_create_empty()
    
    // Set the metadata of the transposed matrix
    array.set(result_matrix, MM_ROWS_INDEX, column_count)
    array.set(result_matrix, MM_COLUMNS_INDEX, row_count)
    
    // Transpose the matrix
    int value_index = na
    
    if (row_count > 0 and column_count > 0)
    //{
        for row = 0 to row_count - 1
        //{
            for column = 0 to column_count - 1
            //{
                value_index := matrix_get_column_index(matrix_array, column) + row
                array.push(result_matrix, array.get(matrix_array, value_index))
            //}
        //}
    //}
    
    result_matrix
//}


// ------------------------------------------------------------------------------------
//    Perform a matrix multiplication and return the result in a new matrix
// ------------------------------------------------------------------------------------

matrix_multiply(matrix_a, matrix_b) =>
//{
    int row_count_a = matrix_get_row_count(matrix_a)
    int row_count_b = matrix_get_row_count(matrix_b)
    int column_count_a = matrix_get_column_count(matrix_a)
    int column_count_b = matrix_get_column_count(matrix_b)
    
    // Create a new matrix to store the result of the operation in
    float[] result_matrix = matrix_create(row_count_a, column_count_b)
    
    if (row_count_a > 0 and column_count_a > 0 and column_count_b > 0 and column_count_a == row_count_b)
    //{
        float value_a = na
        float value_b = na
        float value_result = na
        int result_index = na
        
        // Perform the multiplication matrix operation and store this in the result matrix
        for row_a = 0 to row_count_a - 1
        //{
            for column_b = 0 to column_count_b - 1
            //{
                for column_a = 0 to column_count_a - 1
                //{
                    float temp = 0
                    
                    for row_b = 0 to row_count_b - 1
                    //{
                        float element_a = matrix_get(matrix_a, row_a, row_b)
                        float element_b = matrix_get(matrix_b, row_b, column_b)
                        
                        temp := temp + element_a * element_b
                    //}
                    
                    matrix_set(result_matrix, row_a, column_b, temp)
                //}
            //}
        //}
    //}
    
    result_matrix
//}
// ------------------------------------------------------------------------------------
//    Decomposer the matrix into an orthogonal matrix and an upper triangular matrix,
//    using QR decomposition
//    
//    Original code taken from here (with consent from the author, tbiktag):
//    https://www.tradingview.com/script/0GhsW1KR-Moving-Regression/
// ------------------------------------------------------------------------------------

matrix_get_qr_decomposition(matrix_array) =>
//{
    int row_count = matrix_get_row_count(matrix_array)
    int column_count = matrix_get_column_count(matrix_array)
    
    float[] q = matrix_create(row_count, column_count)
    float[] r = matrix_create(column_count, column_count)
    
    if (row_count > 0 and column_count > 0)
    //{
        float norm = 0
        float value = 0
        
        float[] a = array.new_float(row_count)
        
        // Get the first column of 'matrix_array' and its 2-norm
        for row = 0 to row_count - 1
        //{
            value := matrix_get(matrix_array, row, 0)
            norm := norm + math.pow(value, 2)
            array.set(a, row, value)
        //}
        
        norm := math.sqrt(norm)
        
        // Assign first diagonal element of R and first column of Q
        matrix_set(r, 0, 0, norm)
        
        for row = 0 to row_count - 1
            matrix_set(q, row, 0, array.get(a, row) / norm)
        
        // Repeat for the rest of the columns
        if (column_count > 1)
        //{
            for k = 1 to column_count - 1
            //{
                for row = 0 to row_count - 1
                    array.set(a, row, matrix_get(matrix_array, row, k))
                
                for column = 0 to k - 1
                //{
                    if (column >= k)
                        break
                    
                    // Get R[column, k] as scalar product of Q[column] column an A[k] column
                    norm := 0
                    
                    for row = 0 to row_count - 1
                        norm := norm + matrix_get(q, row, column) * array.get(a, row)
                    
                    matrix_set(r, column, k, norm)
                    
                    // Update vector A
                    for row = 0 to row_count - 1
                    //{
                        value := array.get(a, row) - norm * matrix_get(q, row, column)
                        array.set(a, row, value)
                    //}
                //}
                
                // Get diagonal R[k, k] and Q[k] column
                norm := array_2_norm(a)
                matrix_set(r, k, k, norm)
                
                for row = 0 to row_count - 1
                    matrix_set(q, row, k, array.get(a, row) / norm)
            //}
        //}
        
        [q, r]
    //}
    else
        [matrix_create_empty(), matrix_create_empty()]
//}


// ------------------------------------------------------------------------------------
//    Get the pseudo-inverse from QR
//    
//    Original code taken from here (with consent from the author, tbiktag):
//    https://www.tradingview.com/script/0GhsW1KR-Moving-Regression/
// ------------------------------------------------------------------------------------

matrix_get_pseudo_inverse_from_qr(q, r) =>
//{
    int column_count = matrix_get_column_count(q)
    float[] q_transposed = matrix_transpose(q)
    var r_inv = matrix_create(column_count, column_count)
    float value = 0
    
    matrix_set(r_inv, 0, 0, 1 / matrix_get(r, 0, 0))
    
    if (column_count > 1)
    //{
        for column = 1 to column_count - 1
        //{
            for row = 0 to column - 1
            //{
                value := 0
                
                for k = row to column - 1
                //{
                    if (k >= column)
                        break
                    
                    value := value + matrix_get(r_inv, row, k) * matrix_get(r, k, column)
                //}
                
                matrix_set(r_inv, row, column, value)
            //}
            
            for k = 0 to column - 1
                matrix_set(r_inv, k, column, -matrix_get(r_inv, k, column) / matrix_get(r, column, column))
            
            matrix_set(r_inv, column, column, 1 / matrix_get(r, column, column))
        //}
    //}
    
    matrix_multiply(r_inv, q_transposed)
//}

matrix_get_pseudo_inverse(matrix_array) =>
//{
    [q, r] = matrix_get_qr_decomposition(matrix_array)
    matrix_get_pseudo_inverse_from_qr(q, r)
//}



set_output_location(ta.highest(close, 300), ta.lowest(close, 300))
var power_cache = array.new_float(0)
//var inverse_matrix = array.new_float(0)
//var x_data_matrix_to_solve = matrix_create(int(order)+1, int(order)+1)
// *********** Initialization **************
if barstate.isfirst
    for i = 0 to length-1 //  ***
        array.push(ln,line.new(na,na,na,na))
    for i = 0 to (length-1)/2
        array.push(labels,label.new(na,na,na))
    for dist = 0 to kernel_size
        array.push(kernel,math.exp(-(math.pow(dist/h,2)/2)))
    // prepare the inverse matrices for different points (from border to center) to solve for the constant factor in the polynomial regression
    if order !=0
        x_data_matrix_to_solve = matrix_create(int(order)+1, int(order)+1)
        // add the data from the left-hand-side, i.e. negative x values
        for o = 0 to 2*order
            sum_x_moment = 0.
            for x = - kernel_size to 0
                sum_x_moment += array.get(kernel,math.abs(x)) * math.pow(x/x_scaling,o)
            for index_on_row = 0 to order
                index_on_column = o - index_on_row
                if index_on_column >= 0 and index_on_column <= order
                    matrix_set(x_data_matrix_to_solve, index_on_row, index_on_column, sum_x_moment)
        inverse_matrix = matrix_get_pseudo_inverse(x_data_matrix_to_solve)
        for i = 0 to order
            array.set(weight_matrix, i, matrix_get(inverse_matrix, i, 0))
            //output( str.tostring(matrix_get(inverse_matrix, i, 0)))
        for x = 1 to kernel_size
            for o = 0 to 2*order
                additional_moment_for_data_on_the_right = array.get(kernel,math.abs(x)) * math.pow(x/x_scaling,o)
                for index_on_row = 0 to order
                    index_on_column = o - index_on_row
                    if index_on_column >= 0 and index_on_column <= order
                        matrix_set(x_data_matrix_to_solve, index_on_row, index_on_column, additional_moment_for_data_on_the_right + matrix_get(x_data_matrix_to_solve, index_on_row, index_on_column))
            inverse_matrix := matrix_get_pseudo_inverse(x_data_matrix_to_solve)
            for i = 0 to order
                array.set(weight_matrix, (order+1)*x+i, matrix_get(inverse_matrix, i, 0))
        for x = -kernel_size to kernel_size
            for o = 0 to order
                array.push(power_cache, math.pow(x/x_scaling, o))




//----
line l = na // ****
line up = na
line dn = na
//----
cross_up = 0.
cross_dn = 0.
if barstate.islast

    max_bars_back(src, 900)
    max_bars_back(log_src, 900)
    max_bars_back(time, 900)
    //log_src = array.new_float(length) 
    y = array.new_float(0)
    //for i = 0 to length-1
    //    array.set(log_src, i, math.log(src[i]))
    for i = 0 to length-1
        sum = 0.
        sumw = 0.
        prediction = 0.
        start_index = repaint ? math.max(0, i-kernel_size) : i
        if order == 0
            for j = start_index to (i+kernel_size)
                w = array.get(kernel, math.abs(i-j)) 
                sum +=  log_src[j]*w //math.log(src[j])*w//array.get(log_src, j)*w// 
                sumw += w
            prediction := math.exp(sum/sumw)
        else
            y_moments = array.new_float(order+1, 0) // including the zero moment
            volume_weight_total = 0.
            counts = 0
            for j = start_index to (i+kernel_size) // *********** Although in our mathematical treatment, we have defined data on the right as having x>0, here, larger indices refer to data on the left, and should be treated as negative
                weighted_y = log_src[j]*array.get(kernel, math.abs(i-j))
                power_cache_index = (i-j+kernel_size)*(order+1)
                for o = 0 to order
                    pow = array.get(power_cache, power_cache_index+o)
                    array.set(y_moments, o, array.get(y_moments, o) + weighted_y*pow) 
            weight_matrix_start_index = repaint ? math.min(i, kernel_size)*(order+1) : 0
            for o = 0 to order
                prediction += array.get(y_moments, o)*array.get(weight_matrix, weight_matrix_start_index+o)
            prediction := math.exp(prediction)
        array.push(y,prediction)
    
    i_label = 0
    trend_label = true
    for i = 0 to length-2
        y2 = array.get(y,i+1)
        y1 = array.get(y,i)
        if i <= length-1
            l := array.get(ln,i)
            lset(l ,n-i,y1 ,n-i-1,y2 ,y2 > y1 ? #ff1100 : #39ff14)
        
        label_ = array.get(labels, i_label)
        if i != length-2 and trend_label
            label_ := array.get(labels, i_label)
            y3 = array.get(y,i+2)
            if y1 > y2 and y2 < y3
                labelset(label_, n-i,y1,'▲',col=#00000000,style=label.style_label_up,textcolor=up_col,textalign=text.align_center) 
                i_label += 1
                trend_label:= false
            if y1 < y2 and y2 > y3
                labelset(label_, n-i,y1,'▼',col=#00000000,style=label.style_label_down,textcolor=dn_col,textalign=text.align_center) 
                i_label += 1
                trend_label:= false