Cryprich Moving Regression

// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © tbiktag
//
// Moving Regression Band Breakout strategy
//
// The strategy based on the Moving Regression Prediction Bands indicator:
//
// https://www.tradingview.com/script/zOaMXJ65-Moving-Regression-Prediction-Bands/
//
// Entry condition:
//    Long: the price crosses the Upper Band from below.
//    Short: the price crosses the Lower Band from above.
// Exit condition Long:
//    the price crosses the Lower Band (the Lower Band - ATR, or Central Line)
// Exit condition Short:
//    the price crosses the Upper Band (the Upper Band + ATR, or Central Line)
//
//
//@version=4
strategy("Cryprich Moving Regression ", shorttitle = "Cryp reg2", overlay=true)

matrix_get(_A,_i,_j,_nrows) =>
    // Get the value of the element of an implied 2d matrix
    //input:
    // _A :: array: pseudo 2d matrix _A = [[column_0],[column_1],...,[column_(n-1)]]
    // _i :: integer: row number
    // _j :: integer: column number
    // _nrows :: integer: number of rows in the implied 2d matrix
    array.get(_A,_i+_nrows*_j)

matrix_set(_A,_value,_i,_j,_nrows) =>
    // Set a value to the element of an implied 2d matrix
    //input:
    // _A :: array, changed on output: pseudo 2d matrix _A = [[column_0],[column_1],...,[column_(n-1)]]
    // _value :: float: the new value to be set
    // _i :: integer: row number
    // _j :: integer: column number
    // _nrows :: integer: number of rows in the implied 2d matrix
    array.set(_A,_i+_nrows*_j,_value)

transpose(_A,_nrows,_ncolumns) =>
    // Transpose an implied 2d matrix
    // input:
    // _A :: array: pseudo 2d matrix _A = [[column_0],[column_1],...,[column_(n-1)]]
    // _nrows :: integer: number of rows in _A
    // _ncolumns :: integer: number of columns in _A
    // output:
    // _AT :: array: pseudo 2d matrix with implied dimensions: _ncolums x _nrows
    var _AT = array.new_float(_nrows*_ncolumns,0)
    for i = 0 to _nrows-1
        for j = 0 to _ncolumns-1
            matrix_set(_AT, matrix_get(_A,i,j,_nrows),j,i,_ncolumns)
    _AT

multiply(_A,_B,_nrowsA,_ncolumnsA,_ncolumnsB) =>
    // Calculate scalar product of two matrices
    // input:
    // _A :: array: pseudo 2d matrix
    // _B :: array: pseudo 2d matrix
    // _nrowsA :: integer: number of rows in _A
    // _ncolumnsA :: integer: number of columns in _A
    // _ncolumnsB :: integer: number of columns in _B
    // output:
    // _C:: array: pseudo 2d matrix with implied dimensions _nrowsA x _ncolumnsB
    var _C = array.new_float(_nrowsA*_ncolumnsB,0)
    int _nrowsB = _ncolumnsA
    float elementC= 0.0
    for i = 0 to _nrowsA-1
        for j = 0 to _ncolumnsB-1
            elementC := 0
            for k = 0 to _ncolumnsA-1
                elementC := elementC + matrix_get(_A,i,k,_nrowsA)*matrix_get(_B,k,j,_nrowsB)
            matrix_set(_C,elementC,i,j,_nrowsA)
    _C

vnorm(_X,_n) =>
    //Square norm of vector _X with size _n
    float _norm = 0.0
    for i = 0 to _n-1
        _norm := _norm + pow(array.get(_X,i),2)
    sqrt(_norm)

qr_diag(_A,_nrows,_ncolumns) =>
    //QR Decomposition with Modified Gram-Schmidt Algorithm (Column-Oriented)
    // input:
    // _A :: array: pseudo 2d matrix _A = [[column_0],[column_1],...,[column_(n-1)]]
    // _nrows :: integer: number of rows in _A
    // _ncolumns :: integer: number of columns in _A
    // output:
    // _Q: unitary matrix, implied dimenstions _nrows x _ncolumns
    // _R: upper triangular matrix, implied dimansions _ncolumns x _ncolumns
    var _Q = array.new_float(_nrows*_ncolumns,0)
    var _R = array.new_float(_ncolumns*_ncolumns,0)
    var _a = array.new_float(_nrows,0)
    var _q = array.new_float(_nrows,0)
    float _r = 0.0
    float _aux = 0.0
    //get first column of _A and its norm:
    for i = 0 to _nrows-1
        array.set(_a,i,matrix_get(_A,i,0,_nrows))
    _r := vnorm(_a,_nrows)
    //assign first diagonal element of R and first column of Q
    matrix_set(_R,_r,0,0,_ncolumns)
    for i = 0 to _nrows-1
        matrix_set(_Q,array.get(_a,i)/_r,i,0,_nrows)
    if _ncolumns != 1
        //repeat for the rest of the columns
        for k = 1 to _ncolumns-1
            for i = 0 to _nrows-1
                array.set(_a,i,matrix_get(_A,i,k,_nrows))
            for j = 0 to k-1
                //get R_jk as scalar product of Q_j column and A_k column:
                _r := 0
                for i = 0 to _nrows-1
                    _r := _r + matrix_get(_Q,i,j,_nrows)*array.get(_a,i)
                matrix_set(_R,_r,j,k,_ncolumns)
                //update vector _a
                for i = 0 to _nrows-1
                    _aux := array.get(_a,i) - _r*matrix_get(_Q,i,j,_nrows)
                    array.set(_a,i,_aux)
            //get diagonal R_kk and Q_k column
            _r := vnorm(_a,_nrows)
            matrix_set(_R,_r,k,k,_ncolumns)
            for i = 0 to _nrows-1
                matrix_set(_Q,array.get(_a,i)/_r,i,k,_nrows)
    [_Q,_R]

pinv(_A,_nrows,_ncolumns) =>
//Pseudoinverse of matrix _A calculated using QR decomposition
// Input:
// _A:: array: implied as a (_nrows x _ncolumns) matrix _A = [[column_0],[column_1],...,[column_(_ncolumns-1)]]
// Output:
// _Ainv:: array implied as a (_ncolumns x _nrows) matrix _A = [[row_0],[row_1],...,[row_(_nrows-1)]]
// ----
// First find the QR factorization of A: A = QR,
// where R is upper triangular matrix.
// Then _Ainv = R^-1*Q^T.
// ----
    [_Q,_R] = qr_diag(_A,_nrows,_ncolumns)
    _QT = transpose(_Q,_nrows,_ncolumns)
    // Calculate Rinv:
    var _Rinv = array.new_float(_ncolumns*_ncolumns,0)
    float _r = 0.0
    matrix_set(_Rinv,1/matrix_get(_R,0,0,_ncolumns),0,0,_ncolumns)
    if _ncolumns != 1
        for j = 1 to _ncolumns-1
            for i = 0 to j-1
                _r := 0.0
                for k = i to j-1
                    _r := _r + matrix_get(_Rinv,i,k,_ncolumns)*matrix_get(_R,k,j,_ncolumns)
                matrix_set(_Rinv,_r,i,j,_ncolumns)
            for k = 0 to j-1
                matrix_set(_Rinv,-matrix_get(_Rinv,k,j,_ncolumns)/matrix_get(_R,j,j,_ncolumns),k,j,_ncolumns)
            matrix_set(_Rinv,1/matrix_get(_R,j,j,_ncolumns),j,j,_ncolumns)
    //
    _Ainv = multiply(_Rinv,_QT,_ncolumns,_ncolumns,_nrows)
    _Ainv

mae(_x, _xhat) =>
    // Mean Average Error
    // _x.   :: array float, original data
    // _xhat :: array float, model estimate
    // output
    // _nrmse:: float
    float _mae = 0.0
    if array.size(_x) != array.size(_xhat)
        _mae := na
    else
        _N = array.size(_x)
        for i = 0 to _N-1
            _mae := _mae + abs(array.get(_x,i) - array.get(_xhat,i))/_N
    _mae

mr(_src,_window,_degree) =>
    // Vandermonde matrix with implied dimensions (window x degree+1)
    // Linear form: J = [ [z]^0, [z]^1, ... [z]^degree], with z = [ (1-window)/2 to (window-1)/2 ]
    var _J = array.new_float(_window*(_degree+1),0)
    for i = 0 to _window-1
        for j = 0 to _degree
            matrix_set(_J,pow(i,j),i,j,_window)
    // Vector of raw datapoints:
    var _Y_raw = array.new_float(_window,na)
    for j = 0 to _window-1
        array.set(_Y_raw,j,_src[_window-1-j])
    // Calculate polynomial coefficients which minimize the loss function
    _C = pinv(_J,_window,_degree+1)
    _a_coef = multiply(_C,_Y_raw,_degree+1,_window,1)
    // For smoothing, approximate the last point (i.e. z=window-1) by a0
    float _Y = 0.0
    for i = 0 to _degree
        _Y := _Y + array.get(_a_coef,i)*pow(_window-1,i)
    // Trend Direction Forecast
    float _Y_f = 0.0
    for i = 0 to _degree
        _Y_f := _Y_f + array.get(_a_coef,i)*pow(_window,i)
    // Calculates data estimate (needed for rmse)
    _Y_hat = multiply(_J,_a_coef,_window,_degree+1,1)
    float _err = mae(_Y_raw,_Y_hat)
    [_Y,_Y_f,_err]

/// --- main ---
src = input(title="Source", defval=close, group = "Model Parameters:")
degree = input(title="Local Polynomial Degree", type = input.integer, defval=2, minval = 0,  group = "Model Parameters:")
// RSI
rsiGroup = "RSI"
lenRSI = input(14, title="Length", minval=1, group=rsiGroup)
// lessThan = input(50, title="Less than", minval=1 , maxval=100, group=rsiGroup)
RSI = rsi(src, lenRSI)

// MACD
fast_length = input(title="Fast Length", type=input.integer, defval=12, group="MACD")
slow_length = input(title="Slow Length", type=input.integer, defval=26, group="MACD")
signal_length = input(title="Signal Smoothing", type=input.integer, minval = 1, maxval = 50, defval = 9, group="MACD")
sma_source = input(title="Oscillator MA Type", type=input.string, defval="EMA", options=["SMA", "EMA"], group="MACD")
sma_signal = input(title="Signal Line MA Type", type=input.string, defval="EMA", options=["SMA", "EMA"], group="MACD")
fast_ma = sma_source == "SMA" ? sma(src, fast_length) : ema(src, fast_length)
slow_ma = sma_source == "SMA" ? sma(src, slow_length) : ema(src, slow_length)
macd = fast_ma - slow_ma
signal = sma_signal == "SMA" ? sma(macd, signal_length) : ema(macd, signal_length)

window = input(title="Length (must be larger than degree)", type = input.integer, defval=80, minval = 2,  group = "Model Parameters:")
mult = input(title="Multiplier", type=input.float,defval=2.0,minval=0.0,group = "Model Parameters:", tooltip = "Defines the Band Width.")
doLong = input(title="Allow Long Entries", type=input.bool,defval=true, inline = "linealwd", group = "Allowed Entries:")
doShort = input(title="Allow Short Entries", type=input.bool,defval=false, inline = "linealwd", group = "Allowed Entries:")
stoplong = input(title="Exit Long At", defval="Central Line",options = ["Lower Band","Upper Band - ATR","Central Line","Upper band"],inline = "lineexit", group = "Exit Conditions:")
stopshort = input(title="Exit Short At", defval="Central Line",options = ["Upper Band","Upper Band + ATR","Central Line", "Lower Band"],inline = "lineexit", group = "Exit Conditions:")
fixedstart =input(title="", group = "Fixed Backtest Period Start/End Dates:",inline = "linebac1", type = input.bool, defval = true)
backtest_start=input(title = "", type = input.time, inline = "linebac1",
              group = "Fixed Backtest Period Start/End Dates:",
              defval = timestamp("01 Jan 2017 13:30 +0000"),
              tooltip="If deactivated, backtest staring from the first available price bar.")
fixedend =  input(title="", group = "Fixed Backtest Period Start/End Dates:",
              inline = "linebac2", type = input.bool, defval = false)
backtest_end =input(title = "", type = input.time, inline = "linebac2",
              group = "Fixed Backtest Period Start/End Dates:",
              defval = timestamp("30 Dec 2080 23:30 +0000"),
              tooltip="If deactivated, backtesting ends at the last available price bar.")
istoploss = input(title="Show Stop-Loss Line", type=input.bool,defval=true, group="Additional Options")
slopefilter = input(title="Enter Only When MR Slope is Positive/Negative", type=input.bool,defval=false, group="Additional Options")
issignal = input(title="Show Breakout Signals Labels", type=input.bool,defval=false, group="Additional Options")
clinetype = input(title="Show Central Line As", defval="Previous-Period MR Prediction",
                  options = ["Ribbon","MR","Previous-Period MR Prediction"], group="Additional Options")

[MR,MR_f,div] = mr(src,window,degree)
div := div*mult



// plot bands
highband = MR_f[1]+div[1]
lowband = MR_f[1]-div[1]
centralline = MR_f[1]
hbplt = plot(highband,title='Upper Band',color=#3C94B8,linewidth=2, transp = 40)
lbplt = plot(lowband,title='Lower Band',color=#B83C94,linewidth=2, transp = 40)
fill(hbplt,lbplt,color=#3C94B8)

// plot central line
var plt_color = #3C94B8
if MR_f[1] < MR and clinetype=="Ribbon"
    plt_color := #94B83C
else if MR_f[1] > MR and clinetype=="Ribbon"
    plt_color := #B83C71
mrplt = plot(MR,title='MR',color=plt_color,linewidth=1,
             transp = clinetype=="Previous-Period MR Prediction"?100:20)
mrfplt = plot(centralline,title='MR Prediction',color=plt_color,linewidth=1,
             transp = clinetype=="MR"?100:20)
fill(mrplt,mrfplt,color=plt_color,transp = clinetype!="Ribbon"?100:20)


//  breakout signals
ATR = atr(14)
Breakout = crossover(src,highband)
Breakdown = crossunder(src,lowband)
ATRBreakout = crossover(src,highband + ATR)
ATRBreakdown = crossunder(src,lowband - ATR)
MRSlopesUp = MR>MR[1]
MRSlopesDw = MR<MR[1]
plotshape(Breakout and issignal ? high : na, location=location.abovebar,
          style=shape.triangleup, color=#88d4ca, size=size.tiny, title='Breakout')
plotshape(Breakdown and issignal? low : na, location=location.belowbar,
          style=shape.triangledown, color=#f67e7d, size=size.tiny,title='Breakdown')
goLong =  Breakdown and (slopefilter?MRSlopesUp:true)  and RSI < 30
goShort = Breakout and (slopefilter?MRSlopesDw:true) and RSI > 60
// Strategy Conditions

// trailing stop / stop loss
trailinglong = stoplong=="Central Line"?centralline:(stoplong=="Lower Band"?lowband:lowband-ATR)
trailingshort= stopshort=="Central Line"?centralline:(stopshort=="Upper Band"?highband:highband+ATR)
// if signal has been triggered at previous bar, set position open price and initial SL
var float opprice = 0.0
var float stoploss = 0.0
if (strategy.position_entry_name[1]=="short" or strategy.position_size[1]==0) and goLong[1]
    opprice := open
    stoploss := trailinglong
else if (strategy.position_entry_name[1]=="long" or strategy.position_size[1]==0) and goShort[1]
    opprice := open
    stoploss := trailingshort
// move SL to position opening price if lowband crosses the initial SL
if (strategy.position_entry_name=="long" and crossover(trailinglong, opprice)) or (strategy.position_entry_name=="short" and crossunder(trailingshort, opprice))
    stoploss := opprice
stoplosscolor = strategy.position_entry_name=="long"?#B83C94:(strategy.position_entry_name=="short"?#3C94B8:color.silver)
plot(istoploss?stoploss:na,title='sl',color=stoplosscolor,linewidth=1,transp=50)

endLong = crossunder(src,trailinglong)
endShort = crossover(src,trailingshort)

// backtest
isinrange = (fixedstart?time>=backtest_start:true) and (fixedend?time<=backtest_end:true)

if goLong and isinrange and doLong
    strategy.entry("long",   true)
    alert(syminfo.tickerid+" Long Signal Triggered",alert.freq_once_per_bar_close)
if goShort and isinrange and doShort
    strategy.entry("short",  false)
    alert(syminfo.tickerid+" Short Signal Triggered",alert.freq_once_per_bar_close)
if endLong or crossunder(src,stoploss)
    strategy.close("long")
    alert(syminfo.tickerid+" Exit Signal Triggered for Long Position",alert.freq_once_per_bar_close)
if endShort or crossover(src,stoploss)
    strategy.close("short")
    alert(syminfo.tickerid+" Exit Signal Triggered for Short Position",alert.freq_once_per_bar_close)
if (not isinrange)
    strategy.close_all()