449 lines
13 KiB
MQL5
449 lines
13 KiB
MQL5
//+------------------------------------------------------------------+
|
|
//| utils.mqh |
|
|
//| Copyright 2025, MetaQuotes Ltd. |
|
|
//| https://www.mql5.com |
|
|
//+------------------------------------------------------------------+
|
|
#property copyright "Copyright 2025, MetaQuotes Ltd."
|
|
#property link "https://www.mql5.com"
|
|
#include<Math\Alglib\solvers.mqh>
|
|
//+------------------------------------------------------------------+
|
|
//| Constant and trend used in regression model |
|
|
//+------------------------------------------------------------------+
|
|
|
|
enum ENUM_TREND
|
|
{
|
|
TREND_NONE=0,
|
|
TREND_CONST_ONLY,
|
|
TREND_LINEAR_ONLY,
|
|
TREND_LINEAR_CONST,
|
|
TREND_QUAD_LINEAR_CONST
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| Options for how to handle existing constants |
|
|
//+------------------------------------------------------------------+
|
|
|
|
enum ENUM_HAS_CONST
|
|
{
|
|
HAS_CONST_RAISE=0,
|
|
HAS_CONST_SKIP,
|
|
HAS_CONST_ADD
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| Options for trimming invalid observations |
|
|
//+------------------------------------------------------------------+
|
|
enum ENUM_TRIM
|
|
{
|
|
TRIM_NONE=0,
|
|
TRIM_FORWARD,
|
|
TRIM_BACKWARD,
|
|
TRIM_BOTH
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| options for how to handle original data set |
|
|
//+------------------------------------------------------------------+
|
|
|
|
enum ENUM_ORIGINAL
|
|
{
|
|
ORIGINAL_EX=0,
|
|
ORIGINAL_IN,
|
|
ORIGINAL_SEP
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| Yule walker calculation options |
|
|
//+------------------------------------------------------------------+
|
|
enum ENUM_YW_METHOD
|
|
{
|
|
YW_ADJUSTED=0,//adjusted
|
|
YW_MLE//mle
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| YW result struct |
|
|
//+------------------------------------------------------------------+
|
|
struct YWResult
|
|
{
|
|
double sigma;
|
|
vector rho;
|
|
matrix inv;
|
|
YWResult(void)
|
|
{
|
|
sigma = EMPTY_VALUE;
|
|
rho = vector::Zeros(0);
|
|
inv = matrix::Zeros(0,0);
|
|
}
|
|
~YWResult(void)
|
|
{
|
|
}
|
|
YWResult(YWResult& other)
|
|
{
|
|
sigma = other.sigma;
|
|
rho = other.rho;
|
|
inv = other.inv;
|
|
}
|
|
void operator=(YWResult& other)
|
|
{
|
|
sigma = other.sigma;
|
|
rho = other.rho;
|
|
inv = other.inv;
|
|
}
|
|
};
|
|
//+------------------------------------------------------------------+
|
|
//| sequence |
|
|
//+------------------------------------------------------------------+
|
|
//| The function generates a sequence of double numbers |
|
|
//| with given starting and ending values and step. |
|
|
//| |
|
|
//| Arguments: |
|
|
//| from : Starting value of the sequence |
|
|
//| to : Last value of the sequence |
|
|
//| step : Step of the sequence |
|
|
//| result[] : Array for the sequence values |
|
|
//| |
|
|
//| Return value: true if successful, otherwise false. |
|
|
//+------------------------------------------------------------------+
|
|
bool Sequence(const double from,const double to,const double step,double &result[])
|
|
{
|
|
//--- check NaN
|
|
if(!MathIsValidNumber(from) || !MathIsValidNumber(to) || !MathIsValidNumber(step))
|
|
return(false);
|
|
|
|
if(to<from)
|
|
return(false);
|
|
//--- calculate number of elements and prepare output array
|
|
int count=1+int((to-from)/step);
|
|
|
|
if(ArraySize(result)<count)
|
|
if(ArrayResize(result,count)!=count)
|
|
return(false);
|
|
//--- prepare sequence
|
|
for(int i=0; i<count; i++)
|
|
result[i]=from+i*step;
|
|
//---
|
|
return(true);
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//|helper function : adds selected trend and\or constant |
|
|
//+------------------------------------------------------------------+
|
|
bool addtrend(matrix &in,matrix &out,ENUM_TREND trend=TREND_CONST_ONLY, bool prepend=false, ENUM_HAS_CONST has_const=HAS_CONST_SKIP)
|
|
{
|
|
//---
|
|
ulong trendorder=0;
|
|
//---
|
|
if(trend==TREND_NONE)
|
|
return out.Copy(in);
|
|
//---
|
|
if(trend==TREND_CONST_ONLY)
|
|
trendorder = 0;
|
|
else
|
|
if(trend==TREND_LINEAR_CONST || trend==TREND_LINEAR_ONLY)
|
|
trendorder = 1;
|
|
else
|
|
if(trend==TREND_QUAD_LINEAR_CONST)
|
|
trendorder = 2;
|
|
//---
|
|
ulong nobs = in.Rows();
|
|
//---
|
|
vector tvector,temp;
|
|
//---
|
|
if(!tvector.Resize(nobs) || !temp.Resize(nobs))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
|
|
return false;
|
|
}
|
|
//---
|
|
double sequence[];
|
|
//---
|
|
if(!Sequence(1,nobs,1,sequence) || !tvector.Assign(sequence))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Error ", ::GetLastError());
|
|
return false;
|
|
}
|
|
//---
|
|
matrix trendmat;
|
|
//---
|
|
if(!trendmat.Resize(tvector.Size(),(trend==TREND_LINEAR_ONLY)?1:trendorder+1))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
|
|
return false;
|
|
}
|
|
//---
|
|
for(ulong i = 0; i<trendmat.Cols(); i++)
|
|
{
|
|
temp = MathPow(tvector,(trend==TREND_LINEAR_ONLY)?double(i+1):double(i));
|
|
if(!trendmat.Col(temp,i))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Matrix Assign Error ", ::GetLastError());
|
|
return false;
|
|
}
|
|
}
|
|
//---
|
|
vector ptp = in.Ptp(0);
|
|
//---
|
|
if(!ptp.Min())
|
|
{
|
|
if(has_const==HAS_CONST_RAISE)
|
|
{
|
|
Print("Input matrix contains one or more constant columns");
|
|
return false;
|
|
}
|
|
if(has_const==HAS_CONST_SKIP)
|
|
{
|
|
matrix vsplitted[];
|
|
|
|
ulong parts[] = {1,trendmat.Cols()-1};
|
|
|
|
trendmat.Vsplit(parts,vsplitted);
|
|
|
|
if(!trendmat.Copy(vsplitted[1]))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
|
|
return false;
|
|
}
|
|
}
|
|
}
|
|
//---
|
|
int order = (prepend)?1:-1;
|
|
//---
|
|
if(!out.Resize(trendmat.Rows(),trendmat.Cols()+in.Cols()))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
|
|
return false;
|
|
}
|
|
//---
|
|
ulong j = 0;
|
|
matrix mtemp;
|
|
//---
|
|
if(prepend)
|
|
mtemp.Copy(trendmat);
|
|
else
|
|
mtemp.Copy(in);
|
|
//---
|
|
for(j = 0; j<mtemp.Cols(); j++)
|
|
{
|
|
vector col = mtemp.Col(j);
|
|
|
|
if(!out.Col(col,j))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Matrix Assign Error ", ::GetLastError());
|
|
return false;
|
|
}
|
|
|
|
}
|
|
//---
|
|
ulong minus = j;
|
|
//---
|
|
if(prepend)
|
|
mtemp.Copy(in);
|
|
else
|
|
mtemp.Copy(trendmat);
|
|
//---
|
|
for(; j<out.Cols(); j++)
|
|
{
|
|
vector col = mtemp.Col(j-minus);
|
|
|
|
if(!out.Col(col,j))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Matrix Assign Error ", ::GetLastError());
|
|
return false;
|
|
}
|
|
|
|
}
|
|
//---
|
|
return true;
|
|
//---
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| helper function: transforms rhs matrix |
|
|
//+------------------------------------------------------------------+
|
|
bool lagmat(matrix &in,matrix &out[],ulong mlag,ENUM_TRIM trim=TRIM_BOTH,ENUM_ORIGINAL original=ORIGINAL_IN)
|
|
{
|
|
//---
|
|
ulong nobs=in.Rows();
|
|
ulong nvars =in.Cols();
|
|
ulong dropidx = 0;
|
|
//---
|
|
if(mlag>=nobs)
|
|
{
|
|
Print("mlag should be < number of rows of the input matrix");
|
|
return false;
|
|
}
|
|
//---
|
|
if(original==ORIGINAL_EX || original==ORIGINAL_SEP)
|
|
dropidx = nvars;
|
|
//---
|
|
matrix nn;
|
|
//---
|
|
if(!nn.Resize(nobs + mlag,nvars * (mlag + 1)))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
|
|
return false;
|
|
}
|
|
//---
|
|
nn.Fill(0.0);
|
|
//---
|
|
ulong maxlag=mlag;
|
|
//---
|
|
ulong row_end,row_start,col_end,col_start,j,z;
|
|
row_end=row_start=col_end=col_start=j=z=0;
|
|
//---
|
|
for(ulong k = 0; k<(maxlag+1); k++)
|
|
{
|
|
row_start = maxlag-k;
|
|
row_end = nobs+maxlag-k;
|
|
col_start = nvars*(maxlag-k);
|
|
col_end = nvars*(maxlag-k+1);
|
|
j = 0;
|
|
for(ulong irow = row_start; irow < row_end; irow++, j++)
|
|
{
|
|
z = 0;
|
|
for(ulong icol = col_start; icol < col_end; icol++, z++)
|
|
nn[irow][icol]=in[j][z];
|
|
}
|
|
}
|
|
//---
|
|
ulong startobs,stopobs;
|
|
if(trim==TRIM_NONE || trim==TRIM_FORWARD)
|
|
startobs=0;
|
|
else
|
|
startobs=maxlag;
|
|
//---
|
|
if(trim==TRIM_NONE || trim==TRIM_BACKWARD)
|
|
stopobs=nn.Rows();
|
|
else
|
|
stopobs=nobs;
|
|
//---
|
|
if(dropidx && original==ORIGINAL_SEP)
|
|
ArrayResize(out,2);
|
|
else
|
|
ArrayResize(out,1);
|
|
//---
|
|
if(!out[0].Resize(stopobs-startobs,nn.Cols()-dropidx))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Resize Error ", ::GetLastError());
|
|
return false;
|
|
}
|
|
//---
|
|
for(ulong irow = startobs; irow<stopobs; irow++)
|
|
for(ulong icol=dropidx; icol<nn.Cols(); icol++)
|
|
out[0][irow-startobs][icol-dropidx]=nn[irow][icol];
|
|
|
|
//---
|
|
if(out.Size()>1)
|
|
{
|
|
if(!out[1].Resize(stopobs-startobs,dropidx))
|
|
{
|
|
Print(__FUNCTION__," ",__LINE__," Resize Error ", GetLastError());
|
|
return false;
|
|
}
|
|
//---
|
|
for(ulong irow = startobs; irow<stopobs; irow++)
|
|
for(ulong icol=0; icol<out[1].Cols(); icol++)
|
|
out[1][irow-startobs][icol]=nn[irow][icol];
|
|
}
|
|
//---
|
|
return true;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| toeplitz matrix construction |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> toeplitz(vector<T>& x)
|
|
{
|
|
ulong n = x.Size();
|
|
matrix<T> out(n,n);
|
|
for(ulong i = 0; i<n; ++i)
|
|
for(ulong j = 0; j<n; ++j)
|
|
if(i<=j)
|
|
out[i,j] = x[j-i];
|
|
else
|
|
out[i,j] = x[i-j];
|
|
return out;
|
|
}
|
|
//+------------------------------------------------------------------+
|
|
//| toeplitz matrix construction |
|
|
//+------------------------------------------------------------------+
|
|
template<typename T>
|
|
matrix<T> toeplitz(vector<T>& x, vector<T>& y)
|
|
{
|
|
ulong nx = x.Size();
|
|
ulong ny = y.Size();
|
|
vector<T> z = vector<T>::Zeros(nx+ny-1);
|
|
ulong start = 1;
|
|
for(ulong i = 0; i<z.Size(); ++i)
|
|
{
|
|
if(i<nx)
|
|
z[i] = x[(nx-i)-1];
|
|
else
|
|
z[i] = y[start++];
|
|
}
|
|
|
|
matrix<T> out(nx,ny);
|
|
for(ulong i = 0; i<nx; ++i)
|
|
{
|
|
for(ulong j = 0; j<ny; ++j)
|
|
{
|
|
if(i<j)
|
|
out[i,j] = z[nx+(j-1)-i];
|
|
else
|
|
out[i,j] = z[nx-(i+1)+j];
|
|
}
|
|
}
|
|
return out;
|
|
}
|
|
//+---------------------------------------------------------------------------+
|
|
//| Estimate AR(p) parameters from a sequence using the Yule-Walker equations.|
|
|
//+---------------------------------------------------------------------------+
|
|
YWResult yule_walker(vector& x, long order=1, ENUM_YW_METHOD method =YW_ADJUSTED, long df = -1, bool inv = false, bool demean = true)
|
|
{
|
|
YWResult out;
|
|
vector in = x;
|
|
if(!in.Size())
|
|
{
|
|
Print(__FUNCTION__, " input sequence is empty ");
|
|
return out;
|
|
}
|
|
|
|
if(demean)
|
|
in = in - in.Mean();
|
|
long n = df>0?df:long(in.Size());
|
|
vector r = vector::Zeros(order+1);
|
|
r[0] = pow(in,2.0).Sum()/double(n);
|
|
double adj_needed = (double)(method == YW_ADJUSTED);
|
|
for(long k = 1; k<(order+1); ++k)
|
|
{
|
|
double sum = 0;
|
|
for(long i = 0; i<(long(in.Size()) - k); ++i)
|
|
sum+=(in[i]*in[k+i]);
|
|
r[k] = sum/double(n - k*adj_needed);
|
|
}
|
|
|
|
vector rslice = vector::Zeros(r.Size()-1);
|
|
for(ulong i = 0; i<rslice.Size(); rslice[i] = r[i], ++i);
|
|
matrix R = toeplitz(rslice);
|
|
for(ulong i = 0; i<rslice.Size(); rslice[i] = r[i+1], ++i);
|
|
|
|
double b[],sol[];
|
|
ArrayResize(b,(int)rslice.Size());
|
|
for(uint i = 0; i<b.Size(); b[i] = rslice[i], ++i);
|
|
|
|
CMatrixDouble Rr = R;
|
|
CDenseSolverLSReport ls;
|
|
int opinfo;
|
|
CDenseSolver::RMatrixSolveLS(Rr,Rr.Rows(),Rr.Cols(),b,0,opinfo,ls,sol);
|
|
if(opinfo<0)
|
|
{
|
|
Print(__FUNCTION__ " possible singular matrix. Using pinv ");
|
|
out.rho = (R.PInv()).MatMul(rslice);
|
|
}
|
|
else
|
|
out.rho.Assign(sol);
|
|
double sigmasq = r[0] - (rslice*out.rho).Sum();
|
|
if(MathClassify(sigmasq)==FP_NORMAL && sigmasq>0)
|
|
out.sigma = sqrt(sigmasq);
|
|
else
|
|
out.sigma = double("nan");
|
|
if(inv)
|
|
out.inv = R.Inv();
|
|
return out;
|
|
}
|
|
|
|
//+------------------------------------------------------------------+
|