ayantrader/catch22
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
catch22/Include/Catch22/Catch22.mqh

1288 lines
46 KiB
MQL5
Raw Permalink Normal View History

initial commit
2026-07-12 10:15:50 +00:00
//+------------------------------------------------------------------+
//| Catch22.mqh |
//| MMQ — Muhammad Minhas Qamar |
//| www.mql5.com/en/articles/23488 |
//+------------------------------------------------------------------+
#property copyright "MMQ — Muhammad Minhas Qamar"
#property link "https://www.mql5.com/en/articles/23488"
#property version "1.00"
#property strict
#include <Math\Alglib\alglib.mqh>
//+------------------------------------------------------------------+
//| catch22 — the canonical 22-feature time-series set (Lubba et al. |
//| 2019, DMKD). Each feature is one representative of a cluster of |
//| high-performing hctsa features, chosen to be minimally |
//| redundant. Unless noted, every feature runs on the z-scored |
//| window (mean 0, sd 1), exactly as the reference implementation. |
//| |
//| This single class exposes Compute() to fill a 22-element vector, |
//| plus named getters. Heavy math (FFT) uses ALGLIB; the rest is |
//| hand-rolled so the estimators match the published definitions. |
//+------------------------------------------------------------------+
//--- canonical feature order (matches the reference featureList.txt).
//--- Indices into the output vector so the caller can name a column
//--- without memorizing positions.
enum ENUM_CATCH22
{
C22_DN_HistogramMode_5=0, // mode of distribution, 5 bins
C22_DN_HistogramMode_10, // mode of distribution, 10 bins
C22_CO_f1ecac, // first 1/e crossing of the ACF
C22_CO_FirstMin_ac, // lag of first minimum of the ACF
C22_CO_HistogramAMI_even_2_5, // automutual information, lag 2, 5 bins
C22_CO_trev_1_num, // time-reversal asymmetry statistic
C22_MD_hrv_classic_pnn40, // fraction of |diffs| > 0.04 (pNN40)
C22_SB_BinaryStats_mean_longstretch1, // longest run above the mean
C22_SB_TransitionMatrix_3ac_sumdiagcov, // symbol transition-matrix structure
C22_PD_PeriodicityWang_th0_01, // Wang periodicity measure
C22_CO_Embed2_Dist_tau_d_expfit_meandiff, // embedding trajectory outlierness
C22_IN_AutoMutualInfoStats_40_gaussian_fmmi, // first min of Gaussian AMI
C22_FC_LocalSimple_mean1_tauresrat, // ACF-zero ratio of residuals
C22_DN_OutlierInclude_p_001_mdrmd, // timing of positive outliers
C22_DN_OutlierInclude_n_001_mdrmd, // timing of negative outliers
C22_SP_Summaries_welch_rect_area_5_1, // power in lowest 5th of spectrum
C22_SB_BinaryStats_diff_longstretch0, // longest run of consecutive falls
C22_SB_MotifThree_quantile_hh, // entropy of 3-letter motifs
C22_SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1,// rescaled-range (R/S) scaling
C22_SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1, // DFA scaling exponent
C22_SP_Summaries_welch_rect_centroid, // spectral centroid (Welch)
C22_FC_LocalSimple_mean3_stderr // local mean-3 forecast error
};
#define CATCH22_N 22 // number of features in the canonical set
//+------------------------------------------------------------------+
//| CCatch22 - the public engine. |
//| Call Compute(window, out) with a raw price/return window; it |
//| z-scores the window once, precomputes the shared machinery |
//| (autocorrelation, spectrum, symbolizations) and fills a |
//| 22-element feature vector in canonical order. Named getters read|
//| the last computed vector. |
//+------------------------------------------------------------------+
class CCatch22
{
private:
//--- last computed result
double m_feat[CATCH22_N]; // filled by Compute()
bool m_valid; // true if last Compute() succeeded
//--- shared working buffers (per Compute() call)
double m_x[]; // raw window copy
double m_z[]; // z-scored window (mean 0, sd 1)
int m_n; // window length
double m_mean; // raw window mean
double m_sd; // raw window population sd
double m_acf[]; // autocorrelation, lag 0..m_n-1 (of m_z)
int m_acfN; // number of valid ACF lags computed
//--- shared machinery -------------------------------------------
bool Prepare(const double &window[],int n);
void ComputeACF(void);
double Mean(const double &a[],int n) const;
double Sd(const double &a[],int n,double mean) const;
double Median(const double &a[],int n) const;
int HistogramBinCounts(const double &a[],int n,int nbins,
double &counts[],double &lo,double &binw) const;
void Diff(const double &a[],int n,double &d[]) const;
//--- feature estimators (one per canonical feature) -------------
double HistogramMode(int nbins) const;
double F1ecac(void) const;
double FirstMinAC(void) const;
double HistogramAMI(int lag,int nbins) const;
double Trev(int lag) const;
double HrvPnn40(void) const;
double BinaryStatsMeanLongstretch1(void) const;
double TransitionMatrix3ac(void) const;
double PeriodicityWang(void) const;
double Embed2DistExpfit(void) const;
double AutoMutualInfoFmmi(void) const;
double LocalSimpleMean1Tauresrat(void) const;
double OutlierInclude(int sign) const;
double WelchArea51(void) const;
double BinaryStatsDiffLongstretch0(void) const;
double MotifThreeHH(void) const;
double FluctAnal(int method) const; // 0=R/S, 1=DFA
double SegmentSSR(const double &x[],const double &y[],int a,int b) const;
double WelchCentroid(void) const;
double LocalSimpleMean3Stderr(void) const;
//--- helpers shared by several estimators
int FirstZeroACF(void) const;
void WelchSpectrum(double &psd[],int &npsd) const;
public:
CCatch22(void);
//--- main entry: fills out[CATCH22_N] from a raw window of length n.
bool Compute(const double &window[],int n,double &out[]);
//--- last-result access
bool IsValid(void) const { return m_valid; }
double Feature(int idx) const;
double HistogramMode5(void) const { return Feature(C22_DN_HistogramMode_5); }
double HistogramMode10(void)const { return Feature(C22_DN_HistogramMode_10); }
double DFA(void) const { return Feature(C22_SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1); }
double RSrange(void) const { return Feature(C22_SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1); }
//--- column name for a given index (diagnostics / CSV headers)
static string Name(int idx);
};
//+------------------------------------------------------------------+
//| Construct an empty, not-yet-computed engine. |
//+------------------------------------------------------------------+
CCatch22::CCatch22(void)
{
m_valid=false;
m_n=0;
m_acfN=0;
m_mean=0.0;
m_sd=0.0;
for(int i=0;i<CATCH22_N;i++)
m_feat[i]=0.0;
}
//+------------------------------------------------------------------+
//| Arithmetic mean of the first n elements of a[]. |
//+------------------------------------------------------------------+
double CCatch22::Mean(const double &a[],int n) const
{
if(n<=0)
return 0.0;
double s=0.0;
for(int i=0;i<n;i++)
s+=a[i];
return s/n;
}
//+------------------------------------------------------------------+
//| Population standard deviation of a[] given its mean. |
//+------------------------------------------------------------------+
double CCatch22::Sd(const double &a[],int n,double mean) const
{
if(n<=0)
return 0.0;
double s=0.0;
for(int i=0;i<n;i++)
{
double d=a[i]-mean;
s+=d*d;
}
return MathSqrt(s/n);
}
//+------------------------------------------------------------------+
//| Median of the first n elements of a[] (does not modify a[]). |
//+------------------------------------------------------------------+
double CCatch22::Median(const double &a[],int n) const
{
if(n<=0)
return 0.0;
double t[];
ArrayResize(t,n);
for(int i=0;i<n;i++)
t[i]=a[i];
ArraySort(t);
if((n&1)==1)
return t[n/2];
return 0.5*(t[n/2-1]+t[n/2]);
}
//+------------------------------------------------------------------+
//| First differences of a[]: d[i] = a[i+1]-a[i], length n-1. |
//+------------------------------------------------------------------+
void CCatch22::Diff(const double &a[],int n,double &d[]) const
{
int m=(n>1)?n-1:0;
ArrayResize(d,m);
for(int i=0;i<m;i++)
d[i]=a[i+1]-a[i];
}
//+------------------------------------------------------------------+
//| Copy the window, compute mean/sd, and z-score into m_z. Returns |
//| false on a degenerate (constant or too-short) window: catch22 is |
//| undefined without variance. The reference z-scores internally |
//| for every feature except mean/sd (not in the set), so we do it |
//| once here and all estimators read m_z. |
//+------------------------------------------------------------------+
bool CCatch22::Prepare(const double &window[],int n)
{
if(n<8) // too short for ACF/DFA scales to exist
return false;
m_n=n;
ArrayResize(m_x,n);
ArrayResize(m_z,n);
for(int i=0;i<n;i++)
m_x[i]=window[i];
m_mean=Mean(m_x,n);
m_sd =Sd(m_x,n,m_mean);
if(m_sd<=0.0 || !MathIsValidNumber(m_sd))
return false; // constant window -> undefined
for(int i=0;i<n;i++)
m_z[i]=(m_x[i]-m_mean)/m_sd;
return true;
}
//+------------------------------------------------------------------+
//| Biased autocorrelation of the z-scored window at lags 0..n-1. |
//| r[k] = (1/n) * sum_{i} z[i]*z[i+k] (divide by n, not n-k, to |
//| match the reference CO_ definitions). Stored in m_acf. |
//+------------------------------------------------------------------+
void CCatch22::ComputeACF(void)
{
int n=m_n;
ArrayResize(m_acf,n);
m_acfN=n;
//--- r[0] normalizes to 1 for a z-scored series (variance = 1)
for(int k=0;k<n;k++)
{
double s=0.0;
for(int i=0;i+k<n;i++)
s+=m_z[i]*m_z[i+k];
m_acf[k]=s/n;
}
}
//+------------------------------------------------------------------+
//| Equal-width histogram of a[] into nbins over [min,max]. Fills |
//| counts[nbins]; returns the bin holding the most mass (mode bin). |
//| lo and binw describe the binning for callers that need centers. |
//+------------------------------------------------------------------+
int CCatch22::HistogramBinCounts(const double &a[],int n,int nbins,
double &counts[],double &lo,double &binw) const
{
ArrayResize(counts,nbins);
ArrayInitialize(counts,0.0);
double mn=a[0],mx=a[0];
for(int i=1;i<n;i++)
{
if(a[i]<mn)
mn=a[i];
if(a[i]>mx)
mx=a[i];
}
lo=mn;
double range=mx-mn;
binw=(range>0.0)?range/nbins:1.0;
if(range<=0.0)
return 0;
for(int i=0;i<n;i++)
{
int b=(int)((a[i]-mn)/binw);
if(b<0)
b=0;
if(b>=nbins)
b=nbins-1; // include the right edge in the last bin
counts[b]+=1.0;
}
//--- locate the fullest bin (ties -> first, as in the reference)
int best=0;
double bestc=counts[0];
for(int b=1;b<nbins;b++)
if(counts[b]>bestc)
{
bestc=counts[b];
best=b;
}
return best;
}
//+------------------------------------------------------------------+
//| DN_HistogramMode_5 / _10 |
//| Center of the most populated histogram bin of the z-scored |
//| window, using 'nbins' equal-width bins. A robust location proxy |
//| for the bulk of the distribution. |
//+------------------------------------------------------------------+
double CCatch22::HistogramMode(int nbins) const
{
double counts[],lo,binw;
int best=HistogramBinCounts(m_z,m_n,nbins,counts,lo,binw);
if(binw<=0.0)
return 0.0;
return lo+(best+0.5)*binw; // center of the modal bin
}
//+------------------------------------------------------------------+
//| CO_f1ecac |
//| First lag at which the ACF drops below 1/e. Linear interpolation|
//| between the straddling lags gives a continuous decorrelation |
//| timescale rather than an integer lag. |
//+------------------------------------------------------------------+
double CCatch22::F1ecac(void) const
{
double thresh=1.0/M_E;
for(int k=0;k<m_acfN-1;k++)
{
if(m_acf[k]>thresh && m_acf[k+1]<=thresh)
{
//--- interpolate the crossing between lag k and k+1
double denom=m_acf[k]-m_acf[k+1];
double frac=(denom!=0.0)?(m_acf[k]-thresh)/denom:0.0;
return (double)k+frac;
}
}
return (double)m_n; // never crosses -> saturate at window length
}
//+------------------------------------------------------------------+
//| CO_FirstMin_ac |
//| Lag of the first local minimum of the ACF. Captures the |
//| dominant oscillation half-period when one exists. |
//+------------------------------------------------------------------+
double CCatch22::FirstMinAC(void) const
{
for(int k=1;k<m_acfN-1;k++)
if(m_acf[k]<m_acf[k-1] && m_acf[k]<m_acf[k+1])
return (double)k;
return (double)m_n;
}
//+------------------------------------------------------------------+
//| First lag at which the ACF becomes <= 0 (first zero crossing). |
//| Used by the local-forecast features to bound the horizon. |
//+------------------------------------------------------------------+
int CCatch22::FirstZeroACF(void) const
{
for(int k=1;k<m_acfN;k++)
if(m_acf[k]<=0.0)
return k;
return m_acfN;
}
//+------------------------------------------------------------------+
//| CO_trev_1_num |
//| Numerator of the time-reversal asymmetry statistic at lag 1: |
//| mean( (z[i+1]-z[i])^3 ). |
//| Non-zero values indicate irreversibility (nonlinearity). |
//+------------------------------------------------------------------+
double CCatch22::Trev(int lag) const
{
int m=m_n-lag;
if(m<=0)
return 0.0;
double s=0.0;
for(int i=0;i<m;i++)
{
double d=m_z[i+lag]-m_z[i];
s+=d*d*d;
}
return s/m;
}
//+------------------------------------------------------------------+
//| MD_hrv_classic_pnn40 |
//| pNN40 from heart-rate-variability: proportion of successive |
//| absolute differences exceeding 0.04. Applied to the z-scored |
//| series (the reference divides by 1000 internally on raw data; |
//| on z-scored input the 0.04 threshold is used directly). |
//+------------------------------------------------------------------+
double CCatch22::HrvPnn40(void) const
{
double d[];
Diff(m_z,m_n,d);
int m=ArraySize(d);
if(m<=0)
return 0.0;
int cnt=0;
for(int i=0;i<m;i++)
if(MathAbs(d[i])>0.04)
cnt++;
return (double)cnt/m;
}
//+------------------------------------------------------------------+
//| SB_BinaryStats_mean_longstretch1 |
//| Symbolize z as 1 (above mean) / 0 (below); return the length of |
//| the longest consecutive run of 1s. The z-scored mean is 0, so |
//| "above mean" is simply z[i] > 0. |
//+------------------------------------------------------------------+
double CCatch22::BinaryStatsMeanLongstretch1(void) const
{
int longest=0,run=0;
for(int i=0;i<m_n;i++)
{
if(m_z[i]>0.0)
{
run++;
if(run>longest)
longest=run;
}
else
run=0;
}
//--- hctsa measures the stretch spanned (transitions), i.e. run+1
return (longest>0)?(double)(longest+1):0.0;
}
//+------------------------------------------------------------------+
//| SB_BinaryStats_diff_longstretch0 |
//| Symbolize the sign of successive differences; return the length |
//| of the longest consecutive run of decreases (diff < 0). |
//+------------------------------------------------------------------+
double CCatch22::BinaryStatsDiffLongstretch0(void) const
{
double d[];
Diff(m_z,m_n,d);
int m=ArraySize(d);
int longest=0,run=0;
for(int i=0;i<m;i++)
{
if(d[i]<0.0)
{
run++;
if(run>longest)
longest=run;
}
else
run=0;
}
//--- hctsa measures the stretch spanned (transitions), i.e. run+1
return (longest>0)?(double)(longest+1):0.0;
}
//+------------------------------------------------------------------+
//| CO_HistogramAMI_even_2_5 |
//| Automutual information between z[t] and z[t+lag] estimated from |
//| a 2-D equal-width histogram with 'nbins' bins per axis over the |
//| data range (the "even" binning of the reference): |
//| AMI = sum_ij p_ij * log( p_ij / (p_i * p_j) ). |
//| A nonlinear dependence measure; unlike the ACF it detects |
//| structure invisible to linear correlation. |
//+------------------------------------------------------------------+
double CCatch22::HistogramAMI(int lag,int nbins) const
{
int m=m_n-lag;
if(m<=1 || nbins<=1)
return 0.0;
//--- shared bin edges over the full z range (both axes identical)
double mn=m_z[0],mx=m_z[0];
for(int i=0;i<m_n;i++)
{
if(m_z[i]<mn)
mn=m_z[i];
if(m_z[i]>mx)
mx=m_z[i];
}
double range=mx-mn;
if(range<=0.0)
return 0.0;
double binw=range/nbins;
//--- joint and marginal counts
double pj[];
ArrayResize(pj,nbins*nbins);
ArrayInitialize(pj,0.0);
double pa[];
ArrayResize(pa,nbins);
ArrayInitialize(pa,0.0);
double pb[];
ArrayResize(pb,nbins);
ArrayInitialize(pb,0.0);
for(int i=0;i<m;i++)
{
int ba=(int)((m_z[i] -mn)/binw);
if(ba<0)
ba=0;
if(ba>=nbins)
ba=nbins-1;
int bb=(int)((m_z[i+lag]-mn)/binw);
if(bb<0)
bb=0;
if(bb>=nbins)
bb=nbins-1;
pj[ba*nbins+bb]+=1.0;
pa[ba]+=1.0;
pb[bb]+=1.0;
}
//--- normalize to probabilities and accumulate the MI sum
double ami=0.0;
for(int a=0;a<nbins;a++)
{
double pai=pa[a]/m;
if(pai<=0.0)
continue;
for(int b=0;b<nbins;b++)
{
double pij=pj[a*nbins+b]/m;
if(pij<=0.0)
continue;
double pbi=pb[b]/m;
ami+=pij*MathLog(pij/(pai*pbi));
}
}
return ami;
}
//+------------------------------------------------------------------+
//| IN_AutoMutualInfoStats_40_gaussian_fmmi |
//| First minimum of the Gaussian-estimated automutual information |
//| as a function of lag, scanned up to tau_max=40 (capped at n-1). |
//| Under a Gaussian assumption AMI(k) = -0.5*log(1 - r[k]^2) where |
//| r[k] is the autocorrelation; we return the first lag at which |
//| that curve turns up. |
//+------------------------------------------------------------------+
double CCatch22::AutoMutualInfoFmmi(void) const
{
int taumax=(int)MathMin(40,m_acfN-1);
if(taumax<2)
return (double)taumax;
double ami[];
ArrayResize(ami,taumax+1);
for(int k=1;k<=taumax;k++)
{
double r=m_acf[k];
double t=1.0-r*r;
ami[k]=(t>0.0)?-0.5*MathLog(t):50.0; // large but finite if |r|~1
}
//--- first minimum = first lag where the AMI curve stops decreasing,
//--- i.e. the earliest k>=1 with ami[k+1] > ami[k] (the reference
//--- reports the lag at the turn, indexed from 1).
for(int k=1;k<taumax;k++)
if(ami[k+1]>ami[k])
return (double)k;
return (double)taumax;
}
//+------------------------------------------------------------------+
//| CO_Embed2_Dist_tau_d_expfit_meandiff |
//| Embed the series in 2-D at lag tau (= first zero of the ACF): |
//| points (z[i], z[i+tau]). Take the Euclidean step lengths between|
//| successive embedded points, fit an exponential to their |
//| distribution, and return the mean absolute deviation of the |
//| empirical distribution from that exponential fit. A measure of |
//| how "outlier-prone" trajectories are in phase space. |
//+------------------------------------------------------------------+
double CCatch22::Embed2DistExpfit(void) const
{
int tau=FirstZeroACF();
if(tau<1)
tau=1;
int m=m_n-tau; // number of embedded points
if(m<3)
return 0.0;
//--- successive step lengths in the 2-D embedding
double dist[];
ArrayResize(dist,m-1);
for(int i=0;i<m-1;i++)
{
double dx=m_z[i+1]-m_z[i];
double dy=m_z[i+1+tau]-m_z[i+tau];
dist[i]=MathSqrt(dx*dx+dy*dy);
}
int nd=m-1;
//--- exponential MLE: lambda = 1/mean, pdf f(x)=lambda*exp(-lambda x)
double meand=Mean(dist,nd);
if(meand<=0.0)
return 0.0;
double lambda=1.0/meand;
//--- compare a histogram of the distances to the fitted exponential
//--- density and return the mean absolute difference across bins.
int nbins=(int)MathSqrt((double)nd);
if(nbins<1)
nbins=1;
double counts[],lo,binw;
HistogramBinCounts(dist,nd,nbins,counts,lo,binw);
if(binw<=0.0)
return 0.0;
double sumabs=0.0;
for(int b=0;b<nbins;b++)
{
double emp=counts[b]/(nd*binw); // empirical density
double center=lo+(b+0.5)*binw;
double fit=lambda*MathExp(-lambda*center); // exponential density
sumabs+=MathAbs(emp-fit);
}
return sumabs/nbins;
}
//+------------------------------------------------------------------+
//| FC_LocalSimple_mean1_tauresrat |
//| One-step-ahead forecast using the previous single value |
//| (local mean of window length 1). Compute the residuals, then |
//| return the ratio of the first ACF-zero of the residuals to the |
//| first ACF-zero of the original series. ~1 means the simple |
//| predictor did not change the correlation timescale. |
//+------------------------------------------------------------------+
double CCatch22::LocalSimpleMean1Tauresrat(void) const
{
int m=m_n-1;
if(m<3)
return 0.0;
//--- residuals of the naive "predict previous value" model
double res[];
ArrayResize(res,m);
for(int i=0;i<m;i++)
res[i]=m_z[i+1]-m_z[i];
//--- first ACF-zero of the residuals (biased ACF, as elsewhere)
double rmean=Mean(res,m);
double r0=0.0;
for(int i=0;i<m;i++)
{
double d=res[i]-rmean;
r0+=d*d;
}
r0/=m;
int tauRes=m;
if(r0>0.0)
{
for(int k=1;k<m;k++)
{
double s=0.0;
for(int i=0;i+k<m;i++)
s+=(res[i]-rmean)*(res[i+k]-rmean);
if((s/m)/r0<=0.0)
{
tauRes=k;
break;
}
}
}
int tauOrig=FirstZeroACF();
if(tauOrig<1)
tauOrig=1;
return (double)tauRes/(double)tauOrig;
}
//+------------------------------------------------------------------+
//| FC_LocalSimple_mean3_stderr |
//| Forecast each point as the mean of the previous 3 values; return|
//| the standard error (population sd) of the forecast residuals. |
//| Higher = the series is less predictable from a short local mean.|
//+------------------------------------------------------------------+
double CCatch22::LocalSimpleMean3Stderr(void) const
{
int train=3;
int m=m_n-train;
if(m<=1)
return 0.0;
double res[];
ArrayResize(res,m);
for(int i=0;i<m;i++)
{
double f=(m_z[i]+m_z[i+1]+m_z[i+2])/3.0; // mean of previous 3
res[i]=m_z[i+train]-f;
}
double rmean=Mean(res,m);
return Sd(res,m,rmean);
}
//+------------------------------------------------------------------+
//| DN_OutlierInclude_p_001 / n_001 (mdrmd) |
//| Progressively raise a threshold from 0 upward (sign=+1) or lower|
//| it from 0 downward (sign=-1) in small steps. At each level keep |
//| the timings (indices) of points beyond the threshold, and track |
//| the median of those timings normalized to [-1,1] about the |
//| series midpoint. The feature is the median across thresholds of |
//| that centered median-timing (mdrmd). Detects whether extremes |
//| cluster early or late in the window. |
//+------------------------------------------------------------------+
double CCatch22::OutlierInclude(int sign) const
{
if(m_n<3)
return 0.0;
double inc=0.01; // threshold step (reference: 0.01)
double maxabs=0.0;
for(int i=0;i<m_n;i++)
if(MathAbs(m_z[i])>maxabs)
maxabs=MathAbs(m_z[i]);
if(maxabs<=0.0)
return 0.0;
int nlevels=(int)(maxabs/inc);
if(nlevels<1)
return 0.0;
double mids[];
ArrayResize(mids,nlevels);
int nvalid=0;
double half=(m_n-1)/2.0;
for(int L=1;L<=nlevels;L++)
{
double thr=L*inc;
//--- collect indices beyond the signed threshold
double idx[];
int c=0;
ArrayResize(idx,m_n);
for(int i=0;i<m_n;i++)
{
double v=sign*m_z[i];
if(v>=thr)
idx[c++]=(double)i;
}
if(c==0)
continue;
//--- reference stops once fewer than ~2% of points remain
if((double)c/m_n<0.02)
break;
ArrayResize(idx,c);
double medIdx=Median(idx,c);
//--- center to [-1,1] about the window midpoint
mids[nvalid++]=(medIdx-half)/half;
}
if(nvalid==0)
return 0.0;
ArrayResize(mids,nvalid);
return Median(mids,nvalid);
}
//+------------------------------------------------------------------+
//| SB_MotifThree_quantile_hh |
//| Symbolize z into 3 letters by quantile (terciles: below the |
//| 33rd pct = 0, middle = 1, above the 67th pct = 2). Count the 9 |
//| ordered 2-letter transitions and return the Shannon entropy of |
//| that 3x3 transition distribution (natural log). Higher entropy =|
//| richer short-range symbolic dynamics. |
//+------------------------------------------------------------------+
double CCatch22::MotifThreeHH(void) const
{
if(m_n<3)
return 0.0;
//--- tercile edges from the sorted z-scored window
double s[];
ArrayResize(s,m_n);
for(int i=0;i<m_n;i++)
s[i]=m_z[i];
ArraySort(s);
double q1=s[(int)(0.3333*(m_n-1))];
double q2=s[(int)(0.6667*(m_n-1))];
//--- symbolize into {0,1,2}
int sym[];
ArrayResize(sym,m_n);
for(int i=0;i<m_n;i++)
{
if(m_z[i]<=q1)
sym[i]=0;
else
if(m_z[i]<=q2)
sym[i]=1;
else
sym[i]=2;
}
//--- 3x3 transition counts over consecutive pairs
double cnt[9];
ArrayInitialize(cnt,0.0);
int total=m_n-1;
for(int i=0;i<total;i++)
cnt[sym[i]*3+sym[i+1]]+=1.0;
//--- Shannon entropy of the pair distribution
double h=0.0;
for(int k=0;k<9;k++)
{
double p=cnt[k]/total;
if(p>0.0)
h-=p*MathLog(p);
}
return h;
}
//+------------------------------------------------------------------+
//| SB_TransitionMatrix_3ac_sumdiagcov |
//| Symbolize z into 3 equal-width states, then build the 3x3 |
//| row-normalized transition-probability matrix at a step equal to |
//| the ACF-based downsample tau (first ACF zero). The feature is |
//| the sum of the covariance-matrix diagonal (i.e. the total |
//| variance) of the transition matrix's columns — a compact |
//| summary of how structured the state transitions are. |
//+------------------------------------------------------------------+
double CCatch22::TransitionMatrix3ac(void) const
{
int tau=FirstZeroACF();
if(tau<1)
tau=1;
if(m_n-tau<3)
return 0.0;
//--- 3 equal-width states over the z range
double mn=m_z[0],mx=m_z[0];
for(int i=0;i<m_n;i++)
{
if(m_z[i]<mn)
mn=m_z[i];
if(m_z[i]>mx)
mx=m_z[i];
}
double range=mx-mn;
if(range<=0.0)
return 0.0;
double binw=range/3.0;
int sym[];
ArrayResize(sym,m_n);
for(int i=0;i<m_n;i++)
{
int b=(int)((m_z[i]-mn)/binw);
if(b<0)
b=0;
if(b>2)
b=2;
sym[i]=b;
}
//--- transition counts at step tau, then row-normalize to probabilities
double T[9];
ArrayInitialize(T,0.0);
double rowsum[3];
ArrayInitialize(rowsum,0.0);
int m=m_n-tau;
for(int i=0;i<m;i++)
{
int a=sym[i], b=sym[i+tau];
T[a*3+b]+=1.0;
rowsum[a]+=1.0;
}
for(int a=0;a<3;a++)
if(rowsum[a]>0.0)
for(int b=0;b<3;b++)
T[a*3+b]/=rowsum[a];
//--- sum of the diagonal of the column-covariance matrix of T.
//--- Treat the 3 columns as variables observed over 3 rows; the
//--- covariance diagonal is the per-column variance, summed.
double colmean[3];
for(int b=0;b<3;b++)
{
double s=0.0;
for(int a=0;a<3;a++)
s+=T[a*3+b];
colmean[b]=s/3.0;
}
double sumdiag=0.0;
for(int b=0;b<3;b++)
{
double v=0.0;
for(int a=0;a<3;a++)
{
double d=T[a*3+b]-colmean[b];
v+=d*d;
}
sumdiag+=v/3.0; // population variance of column b
}
return sumdiag;
}
//+------------------------------------------------------------------+
//| PD_PeriodicityWang_th0_01 |
//| Wang's periodicity measure. Spline-detrend is approximated by a |
//| light three-point smoother; then find, in the ACF, the first |
//| peak (local maximum preceded by a trough) whose height exceeds |
//| a small threshold (0.01). Returns that lag as the estimated |
//| fundamental period, or 0 if none qualifies. |
//+------------------------------------------------------------------+
double CCatch22::PeriodicityWang(void) const
{
double th=0.01;
//--- walk the ACF: require a preceding trough, an upturn, and a peak
int i=0;
//--- skip the initial monotone descent from lag 0
while(i<m_acfN-1 && m_acf[i+1]<m_acf[i])
i++;
//--- from here look for the first local maximum above the threshold
for(int k=i+1;k<m_acfN-1;k++)
{
if(m_acf[k]>m_acf[k-1] && m_acf[k]>=m_acf[k+1] && m_acf[k]>th)
return (double)k;
}
return 0.0;
}
//+------------------------------------------------------------------+
//| Welch-style power spectral density of the z-scored window via a |
//| single periodogram (rectangular window), using ALGLIB's real |
//| FFT. psd[j] = |F[j]|^2 for j=0..N/2. The SP_ features summarize |
//| this density. A single segment matches the "rect" variant used |
//| by the reference at these window lengths. |
//+------------------------------------------------------------------+
void CCatch22::WelchSpectrum(double &psd[],int &npsd) const
{
int n=m_n;
double a[];
ArrayResize(a,n);
for(int i=0;i<n;i++)
a[i]=m_z[i];
complex f[];
CAlglib::FFTR1D(a,n,f); // real FFT -> complex spectrum length n
int half=n/2;
npsd=half+1; // 0..Nyquist inclusive
ArrayResize(psd,npsd);
for(int j=0;j<npsd;j++)
{
double re=f[j].real;
double im=f[j].imag;
psd[j]=(re*re+im*im)/n; // periodogram, normalized by n
}
}
//+------------------------------------------------------------------+
//| SP_Summaries_welch_rect_area_5_1 |
//| Fraction of total spectral power contained in the lowest fifth |
//| of the frequency axis. High values indicate low-frequency |
//| (trend-like) dominance. |
//+------------------------------------------------------------------+
double CCatch22::WelchArea51(void) const
{
double psd[];
int np;
WelchSpectrum(psd,np);
if(np<=1)
return 0.0;
double total=0.0;
for(int j=0;j<np;j++)
total+=psd[j];
if(total<=0.0)
return 0.0;
int cut=np/5; // lowest fifth of the frequency bins
if(cut<1)
cut=1;
double low=0.0;
for(int j=0;j<cut;j++)
low+=psd[j];
return low/total;
}
//+------------------------------------------------------------------+
//| SP_Summaries_welch_rect_centroid |
//| Spectral centroid: the power-weighted mean angular frequency of |
//| the periodogram. Frequencies run 0..pi across the np bins. |
//+------------------------------------------------------------------+
double CCatch22::WelchCentroid(void) const
{
double psd[];
int np;
WelchSpectrum(psd,np);
if(np<=1)
return 0.0;
double num=0.0,den=0.0;
for(int j=0;j<np;j++)
{
double w=M_PI*(double)j/(double)(np-1); // angular frequency 0..pi
num+=w*psd[j];
den+=psd[j];
}
return (den>0.0)?num/den:0.0;
}
//+------------------------------------------------------------------+
//| Sum of squared residuals of an OLS line fit to the points |
//| (x[a..b], y[a..b]) inclusive. Used by the piecewise fluctuation |
//| fit below. Returns 0 for a segment of fewer than 3 points (a |
//| line through <=2 points is exact). |
//+------------------------------------------------------------------+
double CCatch22::SegmentSSR(const double &x[],const double &y[],int a,int b) const
{
int m=b-a+1;
if(m<3)
return 0.0;
double sx=0,sy=0,sxx=0,sxy=0;
for(int i=a;i<=b;i++)
{
sx+=x[i];
sy+=y[i];
sxx+=x[i]*x[i];
sxy+=x[i]*y[i];
}
double denom=m*sxx-sx*sx;
if(denom==0.0)
return 0.0;
double slope=(m*sxy-sx*sy)/denom;
double icpt =(sy-slope*sx)/m;
double ssr=0.0;
for(int i=a;i<=b;i++)
{
double fit=icpt+slope*x[i];
double r=y[i]-fit;
ssr+=r*r;
}
return ssr;
}
//+------------------------------------------------------------------+
//| SC_FluctAnal_2_dfa / _rsrangefit (..._logi_prop_r1) |
//| hctsa fluctuation analysis of the integrated z-series. At each |
//| of a set of log-spaced scales s (tau_min..n/2) we measure a |
//| fluctuation F(s): |
//| method 0 (R/S) : mean rescaled range (range/std) of the |
//| profile over non-overlapping windows; |
//| method 1 (DFA) : rms of linear-detrend residuals over |
//| non-overlapping windows. |
//| On the log-log curve (log s, log F) we then locate the "first |
//| linear scaling region": for every interior split point we fit |
//| two separate lines (left and right segments) and pick the split |
//| that minimizes the combined residual sum of squares. The feature|
//| 'prop_r1' is the PROPORTION of scales in that first region, |
//| i.e. (split+1)/nScales — NOT the scaling exponent itself. |
//+------------------------------------------------------------------+
double CCatch22::FluctAnal(int method) const
{
int n=m_n;
//--- integrate the z-scored series into a profile (cumulative sum)
double prof[];
ArrayResize(prof,n);
double run=0.0;
for(int i=0;i<n;i++)
{
run+=m_z[i];
prof[i]=run;
}
//--- log-spaced (logi) scales from tau_min=5 to n/2, deduplicated
int smin=5, smax=n/2;
if(smax<=smin)
return 0.0;
int scales[];
double logS[],logF[];
ArrayResize(scales,64);
ArrayResize(logS,64);
ArrayResize(logF,64);
int npts=0;
double logmin=MathLog((double)smin);
double logmax=MathLog((double)smax);
int nsteps=50; // '50' in the feature name (max scale points)
int lastS=-1;
for(int k=0;k<nsteps && npts<64;k++)
{
double frac=(nsteps>1)?(double)k/(double)(nsteps-1):0.0;
int s=(int)MathRound(MathExp(logmin+frac*(logmax-logmin)));
if(s<=lastS) // enforce strictly increasing integer scales
s=lastS+1;
if(s>smax)
break;
lastS=s;
int nwin=n/s; // non-overlapping windows
if(nwin<1)
continue;
double fsum=0.0;
int used=0;
for(int w=0;w<nwin;w++)
{
int start=w*s;
if(method==1)
{
//--- DFA: linear detrend, rms of residuals
double sx=0,sy=0,sxx=0,sxy=0;
for(int i=0;i<s;i++)
{
double xx=(double)i;
double yy=prof[start+i];
sx+=xx;
sy+=yy;
sxx+=xx*xx;
sxy+=xx*yy;
}
double denom=s*sxx-sx*sx;
double slope=(denom!=0.0)?(s*sxy-sx*sy)/denom:0.0;
double icpt =(s>0)?(sy-slope*sx)/s:0.0;
double ss=0.0;
for(int i=0;i<s;i++)
{
double fit=icpt+slope*i;
double r=prof[start+i]-fit;
ss+=r*r;
}
fsum+=MathSqrt(ss/s);
used++;
}
else
{
//--- R/S: rescaled range of the profile within the window
double mn=prof[start],mx=prof[start],sm=0.0;
for(int i=0;i<s;i++)
{
double v=prof[start+i];
if(v<mn)
mn=v;
if(v>mx)
mx=v;
sm+=v;
}
double mean=sm/s;
double var=0.0;
for(int i=0;i<s;i++)
{
double d=prof[start+i]-mean;
var+=d*d;
}
double sd=MathSqrt(var/s);
if(sd>0.0)
{
fsum+=(mx-mn)/sd;
used++;
}
}
}
if(used>0)
{
double F=fsum/used;
if(F>0.0)
{
scales[npts]=s;
logS[npts]=MathLog((double)s);
logF[npts]=MathLog(F);
npts++;
}
}
}
if(npts<4)
return 0.0;
//--- find the split minimizing two-segment linear residuals. The first
//--- region is scales[0..split]; we require >=2 points per side.
int bestSplit=1;
double bestSSR=DBL_MAX;
for(int split=1;split<npts-2;split++)
{
double ssr=SegmentSSR(logS,logF,0,split)
+SegmentSSR(logS,logF,split+1,npts-1);
if(ssr<bestSSR)
{
bestSSR=ssr;
bestSplit=split;
}
}
//--- proportion of the scale range in the first linear region
return (double)(bestSplit+1)/(double)npts;
}
//+------------------------------------------------------------------+
//| Compute the full 22-feature vector from a raw window of length n.|
//| Z-scores once, precomputes the ACF, then fills every feature in |
//| canonical order. On a degenerate window returns false and leaves|
//| out[] zero-filled. out[] is resized to CATCH22_N. |
//+------------------------------------------------------------------+
bool CCatch22::Compute(const double &window[],int n,double &out[])
{
m_valid=false;
ArrayResize(out,CATCH22_N);
ArrayInitialize(out,0.0);
for(int i=0;i<CATCH22_N;i++)
m_feat[i]=0.0;
if(!Prepare(window,n))
{
for(int i=0;i<CATCH22_N;i++)
out[i]=0.0;
return false;
}
ComputeACF();
//--- fill in canonical order (indices from ENUM_CATCH22)
m_feat[C22_DN_HistogramMode_5] = HistogramMode(5);
m_feat[C22_DN_HistogramMode_10] = HistogramMode(10);
m_feat[C22_CO_f1ecac] = F1ecac();
m_feat[C22_CO_FirstMin_ac] = FirstMinAC();
m_feat[C22_CO_HistogramAMI_even_2_5] = HistogramAMI(2,5);
m_feat[C22_CO_trev_1_num] = Trev(1);
m_feat[C22_MD_hrv_classic_pnn40] = HrvPnn40();
m_feat[C22_SB_BinaryStats_mean_longstretch1] = BinaryStatsMeanLongstretch1();
m_feat[C22_SB_TransitionMatrix_3ac_sumdiagcov] = TransitionMatrix3ac();
m_feat[C22_PD_PeriodicityWang_th0_01] = PeriodicityWang();
m_feat[C22_CO_Embed2_Dist_tau_d_expfit_meandiff] = Embed2DistExpfit();
m_feat[C22_IN_AutoMutualInfoStats_40_gaussian_fmmi] = AutoMutualInfoFmmi();
m_feat[C22_FC_LocalSimple_mean1_tauresrat] = LocalSimpleMean1Tauresrat();
m_feat[C22_DN_OutlierInclude_p_001_mdrmd] = OutlierInclude(+1);
m_feat[C22_DN_OutlierInclude_n_001_mdrmd] = OutlierInclude(-1);
m_feat[C22_SP_Summaries_welch_rect_area_5_1] = WelchArea51();
m_feat[C22_SB_BinaryStats_diff_longstretch0] = BinaryStatsDiffLongstretch0();
m_feat[C22_SB_MotifThree_quantile_hh] = MotifThreeHH();
m_feat[C22_SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1]= FluctAnal(0);
m_feat[C22_SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1] = FluctAnal(1);
m_feat[C22_SP_Summaries_welch_rect_centroid] = WelchCentroid();
m_feat[C22_FC_LocalSimple_mean3_stderr] = LocalSimpleMean3Stderr();
//--- guard against NaN/Inf leaking into downstream ML
for(int i=0;i<CATCH22_N;i++)
{
if(!MathIsValidNumber(m_feat[i]))
m_feat[i]=0.0;
out[i]=m_feat[i];
}
m_valid=true;
return true;
}
//+------------------------------------------------------------------+
//| Read one feature from the last Compute() by canonical index. |
//+------------------------------------------------------------------+
double CCatch22::Feature(int idx) const
{
if(idx<0 || idx>=CATCH22_N)
return 0.0;
return m_feat[idx];
}
//+------------------------------------------------------------------+
//| Canonical feature name for a column index (CSV headers / logs). |
//+------------------------------------------------------------------+
static string CCatch22::Name(int idx)
{
switch(idx)
{
case C22_DN_HistogramMode_5:
return "DN_HistogramMode_5";
case C22_DN_HistogramMode_10:
return "DN_HistogramMode_10";
case C22_CO_f1ecac:
return "CO_f1ecac";
case C22_CO_FirstMin_ac:
return "CO_FirstMin_ac";
case C22_CO_HistogramAMI_even_2_5:
return "CO_HistogramAMI_even_2_5";
case C22_CO_trev_1_num:
return "CO_trev_1_num";
case C22_MD_hrv_classic_pnn40:
return "MD_hrv_classic_pnn40";
case C22_SB_BinaryStats_mean_longstretch1:
return "SB_BinaryStats_mean_longstretch1";
case C22_SB_TransitionMatrix_3ac_sumdiagcov:
return "SB_TransitionMatrix_3ac_sumdiagcov";
case C22_PD_PeriodicityWang_th0_01:
return "PD_PeriodicityWang_th0_01";
case C22_CO_Embed2_Dist_tau_d_expfit_meandiff:
return "CO_Embed2_Dist_tau_d_expfit_meandiff";
case C22_IN_AutoMutualInfoStats_40_gaussian_fmmi:
return "IN_AutoMutualInfoStats_40_gaussian_fmmi";
case C22_FC_LocalSimple_mean1_tauresrat:
return "FC_LocalSimple_mean1_tauresrat";
case C22_DN_OutlierInclude_p_001_mdrmd:
return "DN_OutlierInclude_p_001_mdrmd";
case C22_DN_OutlierInclude_n_001_mdrmd:
return "DN_OutlierInclude_n_001_mdrmd";
case C22_SP_Summaries_welch_rect_area_5_1:
return "SP_Summaries_welch_rect_area_5_1";
case C22_SB_BinaryStats_diff_longstretch0:
return "SB_BinaryStats_diff_longstretch0";
case C22_SB_MotifThree_quantile_hh:
return "SB_MotifThree_quantile_hh";
case C22_SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1:
return "SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1";
case C22_SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1:
return "SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1";
case C22_SP_Summaries_welch_rect_centroid:
return "SP_Summaries_welch_rect_centroid";
case C22_FC_LocalSimple_mean3_stderr:
return "FC_LocalSimple_mean3_stderr";
}
return "unknown";
}
//+------------------------------------------------------------------+
//+------------------------------------------------------------------+
Reference in a new issue Copy permalink