//+------------------------------------------------------------------+ //| 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 //+------------------------------------------------------------------+ //| 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;i1)?n-1:0; ArrayResize(d,m); for(int i=0;i undefined for(int i=0;imx) 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=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;bbestc) { 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;kthresh && 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;k0.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;i0.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;ilongest) 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;imx) 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=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;a0.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;kami[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;i0.0) { for(int k=1;kmaxabs) 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=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;i0.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;imx) 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;i2) 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;i0.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(im_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 complex spectrum length n int half=n/2; npsd=half+1; // 0..Nyquist inclusive ArrayResize(psd,npsd); for(int j=0;j0.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;i1)?(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;w0)?(sy-slope*sx)/s:0.0; double ss=0.0; for(int i=0;imx) mx=v; sm+=v; } double mean=sm/s; double var=0.0; for(int i=0;i0.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=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"; } //+------------------------------------------------------------------+ //+------------------------------------------------------------------+