// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
//
// Moving Regression
// © tbiktag
//
// MR is a generalization of moving average and polynomial regression.
// The procedure approximates a specified number of prior data points with
// a polynomial function of a user-defined degree. Polynomial interpolation
// of the last data point is used to construct a MR time series.
// The color of the MR curve will be green if the local polynomial is predicted
// to move up in the next time step, and red otherwise.
//
//@version=4
//
study("Moving Regression", shorttitle = "MR", 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)
_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
/// --- main ---
src = input(title="Source", defval=close)
degree = input(title="Local Polynomial Degree", type = input.integer, defval=2, minval = 0)
window = input(title="Length (must be larger than degree)", type = input.integer, defval=18, minval = 2)
isforecast = input(title="Predict Trend Direction", type = input.bool, defval=true)
// 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 ]
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:
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)
var plt_color = color.navy
if Y_f > Y and isforecast
plt_color := #006548
else if Y_f < Y and isforecast
plt_color := #CD5C5C
plot(Y,title='Smoothed',color=plt_color,linewidth=2)
.
