Article-22258-Volatility-Mo.../Include/VolatilityModels/Arch/Univariate/acf.mqh

357 lines
10 KiB
MQL5
Raw Permalink Normal View History

2026-06-03 20:03:04 +02:00
//+------------------------------------------------------------------+
//| acf.mqh |
//| Copyright 2025, MetaQuotes Ltd. |
//| https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2025, MetaQuotes Ltd."
#property link "https://www.mql5.com"
#include<VolatilityModels\np.mqh>
#include<Math\Stat\ChiSquare.mqh>
//+------------------------------------------------------------------+
//| autocorrelation function result struct |
//+------------------------------------------------------------------+
struct ACFResult
{
vector acf;
vector qstats;
vector pvalues;
matrix conf_intervals;
ACFResult(void)
{
acf = qstats = pvalues = vector::Zeros(0);
conf_intervals = matrix::Zeros(0,0);
}
ACFResult(vector& _acf, vector& _qstats, vector& _pvalues, matrix& _confintervals)
{
acf = _acf;
qstats = _qstats;
pvalues = _pvalues;
conf_intervals = _confintervals;
}
ACFResult(ACFResult& other)
{
acf = other.acf;
qstats = other.qstats;
pvalues = other.pvalues;
conf_intervals = other.conf_intervals;
}
void operator=(ACFResult& other)
{
acf = other.acf;
qstats = other.qstats;
pvalues = other.pvalues;
conf_intervals = other.conf_intervals;
}
};
//+------------------------------------------------------------------+
//| ACF function |
//+------------------------------------------------------------------+
ACFResult acf(vector& x, ulong nlags = 0, double alpha = 0.05, bool bartlett_conf=false, bool demean= true, bool adjusted = false)
{
ACFResult out;
ulong nobs = x.Size();
if(!nobs)
{
Print(__FUNCTION__, " input container is empty ");
return out;
}
if(!nlags)
nlags = MathMin(ulong(10*log10(nobs)), nobs-1);
vector avf = acovf(x,0,demean,adjusted);
if(!avf.Size())
return out;
out.acf = np::sliceVector(avf,0,long(nlags+1));
out.acf /=avf[0];
int ec = 0;
if(bartlett_conf)
{
vector varacf = vector::Ones(out.acf.Size())/double(nobs);
varacf[0] = 0.;
varacf[1] = 1./double(nobs);
vector temp = np::sliceVector(out.acf,1,out.acf.Size()-1);
vector tempsq = pow(temp,2.);
temp = np::sliceVector(varacf,2);
temp*=1.0+2.*tempsq.CumSum();
np::vectorCopy(varacf,temp,2);
vector interval = CAlglib::InvNormalCDF(1. - alpha/2.)*sqrt(varacf);
out.conf_intervals = matrix::Zeros(2,varacf.Size());
out.conf_intervals.Row(out.acf-interval,0);
out.conf_intervals.Row(out.acf+interval,1);
}
else
{
double varacf = 1./double(x.Size());
double nv = CAlglib::InvNormalCDF(1. - alpha/2.);
double interval = nv*sqrt(varacf);
out.conf_intervals = matrix::Zeros(2,out.acf.Size());
out.conf_intervals.Row(out.acf-interval,0);
out.conf_intervals.Row(out.acf+interval,1);
}
vector temp = np::sliceVector(out.acf,1);
q_stat(temp,nobs,out.qstats,out.pvalues);
return out;
}
/*//+------------------------------------------------------------------+
//|Find the next regular number greater than or equal to target. |
//+------------------------------------------------------------------+
long next_target(long target)
{
if(target<=6)
return target;
if(!(target&(target-1)))
return target;
double match = AL_POSINF;
long p5 = 1;
long p35,quotient,N,p2;
while(p5<target)
{
p35 = p5;
while(p35<target)
{
quotient =-(-target/p35);
p2=(long)pow(2,(long)bit_length(quotient-1));
N = p2*p35;
if(N == target)
return N;
else
if(double(N)<match)
match = double(N);
p35*=3;
if(p35==target)
return p35;
}
p5*=5;
if(p5==target)
return p5;
}
if(double(p5)<match)
match = double(p5);
return long(match);
}
//+------------------------------------------------------------------+
//| get the number bits for an integer |
//+------------------------------------------------------------------+
ulong bit_length(long x)
{
long var = MathAbs(x);
ulong count = 1;
while(var>=2)
{
var = var/2;
++count;
}
return count;
}*/
//+------------------------------------------------------------------+
//|acovf |
//+------------------------------------------------------------------+
vector acovf(vector &in, ulong lags = 0,bool demean= true, bool adjusted = false)
{
vector xo,acov;
if(demean)
xo = in - in.Mean();
else
xo = in;
ulong n = in.Size();
ulong laglen = lags;
if(!laglen)
laglen = n-1;
else
if(lags>n-1)
{
Print(__FUNCTION__, " nlag must be smaller than ",n - 1);
return vector::Zeros(0);
}
acov = vector::Zeros(laglen+1);
acov[0] = xo.Dot(xo);
if(lags)
{
vector a,b;
for(ulong i = 0; i<laglen; ++i)
{
a = np::sliceVector(xo,long(i+1));
b = np::sliceVector(xo,0,long(xo.Size()-(i+1)));
acov[i+1] = a.Dot(b);
}
if(adjusted)
acov/=(double(n) - np::arange(laglen+1));
else
acov/=double(n);
return acov;
}
vector xi,temp,d;
if(adjusted)
{
xi = np::arange(n,1.);
d = vector::Zeros((xi.Size()*2)-1);
np::vectorCopy(d,xi,0,xi.Size());
temp = np::sliceVector(xi,0,xi.Size()-1);
np::reverseVector(temp);
np::vectorCopy(d,temp,xi.Size());
}
else
d = double(n)*vector::Ones(2*n-1);
/*if(fft)
{
ulong nobs = xo.Size();
n = ulong(next_target(long(2*nobs+1)));
CRowComplex Frf;
CRowDouble inx = xo;
CFastFourierTransform::FFTR1D(inx,int(nobs),Frf);
vectorc frf = Frf.ToVector();
frf = frf*frf.Conjugate();
Frf = frf;
CFastFourierTransform::FFTR1DInv(Frf,Frf.Size(),inx);
temp = inx.ToVector();
temp = np::sliceVector(temp,0,nobs);
acov = temp/np::sliceVector(d,long(nobs-1));
}
else
{*/
acov = xo.Correlate(xo,VECTOR_CONVOLVE_FULL);
acov = np::sliceVector(acov,long(n-1))/np::sliceVector(d,long(n-1));
//}
if(lags)
return np::sliceVector(acov,0,long(laglen+1));
return acov;
}
//+------------------------------------------------------------------+
//|LjungBox q statistic |
//+------------------------------------------------------------------+
bool q_stat(vector& x, ulong nvars, vector& out_stats, vector& out_pvals)
{
vector r = np::arange(x.Size(),1.) ;
r = (1.0/(double(nvars) - r))*pow(x,2.);
out_stats = double(nvars)*double(nvars+2)*r.CumSum();
out_pvals = out_stats;
r = np::arange(x.Size(),1.);
int error_note = 0;
for(ulong i = 0; i<out_pvals.Size(); ++i)
out_pvals[i] = 1. - MathCumulativeDistributionChiSquare(out_stats[i],r[i],error_note);
if(error_note)
return false;
return true;
}
//+------------------------------------------------------------------+
//| plot the ACF |
//+------------------------------------------------------------------+
void plot_acf(ACFResult& acf,long chart_id = 0,int x_coordinate = 0, int y_coordinate = 0, int subwindow = 0,int width = 600, int height = 400,string plot_name = NULL, int fontsize = 15, long viewtimeseconds = 30, bool showchart=true )
{
if(!acf.acf.Size())
{
Print(__FUNCTION__, " Results object is empty ");
return;
}
if(StringLen(plot_name)<1)
plot_name = "Autocorrelation Function";
CGraphic* graph = new CGraphic();
//---
if(!chart_id)
chart_id = ChartID();
//---
ChartRedraw(chart_id);
//---
ChartSetInteger(chart_id, CHART_SHOW, showchart);
//---
if(!graph.Create(chart_id,plot_name,subwindow,x_coordinate,y_coordinate,width,height))
{
delete graph;
Print(__FUNCTION__," Failed to Create graphical object on the Main chart Err ", GetLastError());
return;;
}
//---
double x[], y[];
ArrayResize(x,(int)acf.acf.Size());
ArrayResize(y,(int)acf.acf.Size());
//---
int ncurves = 4;
//---
ENUM_CURVE_TYPE ctype = WRONG_VALUE;
//---
for(int i=0; i<ncurves; ++i)
{
for(uint j = 0; j<x.Size(); ++j)
{
if(i<2)
{
x[j] = double(j);
y[j] = acf.acf[j];
ctype = (i==0)?CURVE_HISTOGRAM:CURVE_POINTS;
}
else
{
ctype = CURVE_LINES;
if(i == 2)
y[j] = (acf.conf_intervals[1,j]-acf.conf_intervals[0,j])/2.;
else
y[j] = (acf.conf_intervals[0,j]-acf.conf_intervals[1,j])/2.;
}
}
//---
CCurve* newcurve = graph.CurveAdd(x, y,ctype,NULL);
//---
switch(ctype)
{
case CURVE_HISTOGRAM:
newcurve.Color(ColorToARGB(clrBlue));
newcurve.HistogramWidth(2);
break;
case CURVE_POINTS:
newcurve.PointsFill(true);
newcurve.PointsColor(ColorToARGB(clrBlue));
newcurve.Color(ColorToARGB(clrBlue));
newcurve.PointsSize(8);
break;
default:
newcurve.Color(ColorToARGB(clrRed));
newcurve.LinesStyle(STYLE_DOT);
break;
}
}
//---
graph.BackgroundMain(plot_name);
graph.BackgroundMainSize(fontsize+5);
graph.BackgroundMainColor(ColorToARGB(clrBlack));
//---
graph.HistoryNameWidth(0);
graph.HistorySymbolSize(0);
graph.HistoryNameSize(0);
//---
graph.XAxis().Name("Lags");
graph.XAxis().NameSize(fontsize);
//---
graph.YAxis().Name("Correlation Coefficient");
graph.YAxis().NameSize(fontsize);
//---
graph.CurvePlotAll();
graph.Update();
//---
Sleep(int(viewtimeseconds)*1000);
graph.Destroy();
delete graph;
ChartRedraw(chart_id);
//---
}