Article-22714-Volatility-Mo.../Include/slsqp_article/Regression/utils.mqh

449 lines
13 KiB
MQL5
Raw Permalink Normal View History

2026-06-03 20:14:05 +02:00
//+------------------------------------------------------------------+
//| 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;
}
//+------------------------------------------------------------------+