Machine Learning: Anchored Gaussian Process Regression [LuxAlgo]
// This work is licensed under a Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0) https://creativecommons.org/licenses/by-nc-sa/4.0/
// © LuxAlgo

//@version=5
indicator("Machine Learning: Anchored Gaussian Process Regression [LuxAlgo]", "LuxAlgo - Machine Learning: Anchored Gaussian Process Regression", overlay = true)
//------------------------------------------------------------------------------
//Settings
//-----------------------------------------------------------------------------{
iTim       = input.time  (timestamp("20 Oct 2023"),   'Anchor 1'  , confirm = true     , inline =  'window1')
iTim_      = input.time  (timestamp("20 Oct 2023"),   'Anchor 2'  , confirm = true     , inline =  'window2')

fitCss     = input.color (#2962ff  , 'Training Window Colours                        ' , inline =  'window3')
fitCssF    = input.color (#2962ff13,           ''                                      , inline =  'window3')

forecast   = input.int   (20       ,   'Forecasting Length   '    ,  minval = 0        , inline = 'forecast')
fcastCss   = input.color (#f23645  ,           ''                                      , inline = 'forecast')
fcastCssF  = input.color (#f2364513,           ''                                      , inline = 'forecast')

length     = input.float (20.      ,   'Smooth                   ',  minval = 1                             )
sigma      = input.float (0.1      ,           step = 0.1         ,  minval = 0                             )

iTime      = math.min(iTim, iTim_)
iTime_     = math.max(iTim, iTim_)

n          =      bar_index
Tdiff_bars = (iTime_ - iTime) / (timeframe.in_seconds(timeframe.period) * 1000)

//-----------------------------------------------------------------------------}
//Functions
//-----------------------------------------------------------------------------{
//Radial basis function
rbf(x1, x2, l)=> math.exp(-math.pow(x1 - x2, 2) / (2.0 * math.pow(l, 2)))

//Kernel matrix
kernel_matrix(X1, X2, l)=>
    km = matrix.new<float>(X1.size(), X2.size())

    i = 0
    for x1 in X1
        j = 0
        for x2 in X2
            rbf = rbf(x1, x2, l)
            km.set(i, j, rbf)
            j += 1
        i += 1
    
    km

//-----------------------------------------------------------------------------}
//Lines
//-----------------------------------------------------------------------------{
line.new(iTime , close, iTime , close + syminfo.mintick
       , color   = color.new(color.silver, 50)
       , xloc    = xloc.bar_time 
       , extend  = extend.both
      )
line.new(iTime_, close, iTime_, close + syminfo.mintick
       , color   = color.new(color.yellow, 50)
       , xloc    = xloc.bar_time 
       , extend  = extend.both
      )

//-----------------------------------------------------------------------------}
//Kernel matrix
//-----------------------------------------------------------------------------{
var int window = 1
var matrix <int>  identity = matrix.new<int>(1, 1, 0)
var matrix<float> K_source = na

if time == iTime 
    window   := math.max(1, math.min(325 - forecast, Tdiff_bars))

    identity := matrix.new<int>(window, window, 0)
    xtrain    =  array.new<int>(0)
    xtest     =  array.new<int>(0)

    //Build identity matrix
    for i = 0 to window -1
        for j = 0 to window -1
            identity.set(i, j, i == j ? 1 : 0)
        xtrain.push(i)
    
    for i = 0 to window + forecast -1
        xtest .push(i)
    
    //Compute kernel matrices
    s         = identity.mult(sigma * sigma)
    Ktrain    = kernel_matrix(xtrain, xtrain, length).sum(s)

    //matrix.pinv() -> The function is calculated using a Moore–Penrose inverse formula based on singular-value decomposition of a matrix. 
    //For non-singular square matrices this function returns the result of matrix.inv.
    K_inv     = Ktrain.pinv() 
    
    K_star    = kernel_matrix(xtrain, xtest , length)
    K_source := K_star.transpose().mult(K_inv)

//-----------------------------------------------------------------------------}
//Set forecast
//-----------------------------------------------------------------------------{
var polyline poly = polyline.new( na ), points = array.new<chart.point>()
var polyline paly = polyline.new( na ), paints = array.new<chart.point>()
float        mean =       ta.sma(close, window                          )

if time == iTime_ 
    //Remove previous polylines
    poly.delete(), paly.delete()

    //Dataset
    ytrain = array.new<float>(0)
    
    for i = 0 to window -1
        ytrain.unshift(close[i] - mean)

    //Estimate
    mu       = K_source.mult(ytrain)

    //Set forecast
    float y1 = na
    k        = -window +2

    for element in mu
        switch 
            k == 1 => points.push(chart.point.from_index(n+k, element + mean)), 
                      paints.push(chart.point.from_index(n+k, element + mean))
            k <  1 => points.push(chart.point.from_index(n+k, element + mean))
            =>        paints.push(chart.point.from_index(n+k, element + mean))

        y1 := element + mean
        k  += 1

    poly   := polyline.new(points
           , curved     = true
           , closed     = false
           , line_color = fitCss
           , line_width = 2
           , fill_color = fitCssF)
    paly   := polyline.new(paints
           , curved     = true
           , closed     = false
           , line_color = fcastCss
           , line_width = 2
           , fill_color = fcastCssF)
          

//-----------------------------------------------------------------------------}
//Plot Moving Average
//-----------------------------------------------------------------------------{
plot(mean, 'MA', color=color.new(color.silver, 90), linewidth=5)

//-----------------------------------------------------------------------------}