时间序列主要特性的分析·进阶篇
（2/3）· 自己写统计脚本太耗时？这篇进阶篇把序列诊断做成免参数通用工具
从峰度到频谱：时序统计量的桶化与递归
这段代码把一段已排序的价格序列 TSort 先按峰度 Kurt 和样本数 NumTS 算直方图柱数：n 取 MathRound((Kurt+1.5)*NumTS^0.4/6.0)，若为偶数就减一，且下限锁在 5。柱数过少会让分布形状失真，过多则样本被稀释，5 柱是经验底限。 直方图横轴用 (v+(i+0.5)*delta-Mean)/StDev 做标准化，纵轴 YHist 先数落入每柱的样本数，再除以 NumTS*delta 乘 StDev 转成密度。开 MT5 把这段塞进脚本，换不同品种看 EURUSD 与 XAUUSD 的峰度差异，柱数可能从 7 跳到 11。 自相关部分用 cor[i-1]=c/a 算滞后 i 的归一化自相关，a 是序列中心化后的平方和。频谱则取 320 个 X 点，按 Spect[i]=2.0*(1+2*b) 用余弦窗累加 ACF，b 里 ((double)NLags-k)/(NLags+1.0) 是三角窗，能压住长滞后的噪声。 LevinsonRecursion 是核心：拿自相关 R[] 跑莱文逊-德宾，每阶 m 算偏相关 K[m]=km，同时更新 AR 系数 A[]，Em 按 (1-km^2)*Em 收缩。外汇与贵金属杠杆高、跳空频繁，AR 系数在重大数据前可能瞬间失效，验证时务必用分段样本。
n=(class="type">int)MathRound((Kurt+class="num">1.5)*MathPow(NumTS,class="num">0.4)/class="num">6.0); if((n&0x01)==class="num">0)n--; if(n<class="num">5)n=class="num">5; class=class="str">"cmt">// 柱数 ArrayResize(XHist,n); ArrayResize(YHist,n); ArrayInitialize(YHist,class="num">0.0); a=MathAbs(TSort[class="num">0]-Mean); b=MathAbs(TSort[NumTS-class="num">1]-Mean); if(a<b)a=b; v=Mean-a; delta=class="num">2.0*a/n; for(i=class="num">0;i<n;i++)XHist[i]=(v+(i+class="num">0.5)*delta-Mean)/StDev; class=class="str">"cmt">// 直方图. X-轴 for(i=class="num">0;i<NumTS;i++) { k=(class="type">int)((TS[i]-v)/delta); if(k>(n-class="num">1))k=n-class="num">1; YHist[k]++; } for(i=class="num">0;i<n;i++)YHist[i]=YHist[i]/NumTS/delta*StDev; class=class="str">"cmt">// 直方图. Y-轴 . . . . . . ArrayResize(cor,IP); a=class="num">0; for(i=class="num">0;i<NumTS;i++)a+=TSCenter[i]*TSCenter[i]; for(i=class="num">1;i<=IP;i++) { c=class="num">0; for(k=i;k<NumTS;k++)c+=TSCenter[k]*TSCenter[k-i]; cor[i-class="num">1]=c/a; class=class="str">"cmt">// 自动相关性 } . . . . . . n=class="num">320; class=class="str">"cmt">// X-点数量 ArrayResize(Spect,n); v=M_PI/n; for(i=class="num">0;i<n;i++) { a=i*v; b=class="num">0; for(k=class="num">0;k<NLags;k++)b+=((class="type">class="kw">double)NLags-k)/(NLags+class="num">1.0)*ACF[k]*MathCos(a*(k+class="num">1)); Spect[i]=class="num">2.0*(class="num">1+class="num">2*b); class=class="str">"cmt">// Y-轴频谱 } . . . class=class="str">"cmt">//----------------------------------------------------------------------------------- class=class="str">"cmt">// 为自动相关序列 R[] 计算莱文逊-德宾递归 class=class="str">"cmt">// 并返回自动回归系数 A[] 和部分自动相关 class=class="str">"cmt">// 系数 K[] class=class="str">"cmt">//----------------------------------------------------------------------------------- class="type">void TSAnalysis::LevinsonRecursion(const class="type">class="kw">double &R[],class="type">class="kw">double &A[],class="type">class="kw">double &K[]) { class="type">int p,i,m; class="type">class="kw">double km,Em,Am1[],err; p=ArraySize(R); ArrayResize(Am1,p); ArrayInitialize(Am1,class="num">0); ArrayInitialize(A,class="num">0); ArrayInitialize(K,class="num">0); km=class="num">0; Em=class="num">1; for(m=class="num">0;m<p;m++) { err=class="num">0; for(i=class="num">0;i<m;i++)err+=Am1[i]*R[m-i-class="num">1]; km=(R[m]-err)/Em; K[m]=km; A[m]=km; for(i=class="num">0;i<m;i++)A[i]=(Am1[i]-km*Am1[m-i-class="num">1]); Em=(class="num">1-km*km)*Em; ArrayCopy(Am1,A); } class="kw">return; } class=class="str">"cmt">//----------------------------------------------------------------------------------- class=class="str">"cmt">// 按频率抽取 (DIF) 的快速哈特利变换 (FHT) 基2算法. class=class="str">"cmt">// Length 等于 N = class="num">2 ** ldn class=class="str">"cmt">//-----------------------------------------------------------------------------------
「快速哈特莱变换的位反转重排」
上面的 fht 函数是快速哈特莱变换（FHT）的核心实现，输入数组 f[] 长度必须是 2 的 ldn 次幂，例如 ldn=10 时处理 1024 个点。前半段按 2 的幂次做蝶形迭代，用 MathSin/MathCos 做旋转因子累加，后半段才做位反转置换。 位反转那段容易被忽略：当 n>2 时，用 do-while 循环生成逆序索引 r，若 r>i 就交换 f[i] 与 f[r]。这一步不跑，频谱顺序就是乱的，直接画图会看到镜像错乱。 在 MT5 里把这段代码贴进自定义类，传一段 EURUSD 的 H1 收盘价数组（长度取 2 的幂），跑完打印 f[] 前 8 个值，能立刻验证重排是否生效。外汇与贵金属杠杆高，此类信号仅作概率参考，实盘须控仓。
class="type">void TSAnalysis::fht(class="type">class="kw">double &f[], class="type">ulong ldn) { const class="type">ulong n = ((class="type">ulong)class="num">1<<ldn); for (class="type">ulong ldm=ldn; ldm>=class="num">1; --ldm) { const class="type">ulong m = ((class="type">ulong)class="num">1<<ldm); const class="type">ulong mh = (m>>class="num">1); const class="type">ulong m4 = (mh>>class="num">1); const class="type">class="kw">double phi0 = M_PI / (class="type">class="kw">double)mh; for (class="type">ulong r=class="num">0; r<n; r+=m) { for (class="type">ulong j=class="num">0; j<mh; ++j) { class="type">ulong t1 = r+j; class="type">ulong t2 = t1+mh; class="type">class="kw">double u = f[t1]; class="type">class="kw">double v = f[t2]; f[t1] = u + v; f[t2] = u - v; } class="type">class="kw">double ph = class="num">0.0; for (class="type">ulong j=class="num">1; j<m4; ++j) { class="type">ulong k = mh-j; ph += phi0; class="type">class="kw">double s=MathSin(ph); class="type">class="kw">double c=MathCos(ph); class="type">ulong t1 = r+mh+j; class="type">ulong t2 = r+mh+k; class="type">class="kw">double pj = f[t1]; class="type">class="kw">double pk = f[t2]; f[t1] = pj * c + pk * s; f[t2] = pj * s - pk * c; } } } if(n>class="num">2) { class="type">ulong r = class="num">0; for (class="type">ulong i=class="num">1; i<n; i++) { class="type">ulong k = n; do {k = k>>class="num">1; r = r^k;} class="kw">while ((r & k)==class="num">0); if (r>i) {class="type">class="kw">double tmp = f[i]; f[i] = f[r]; f[r] = tmp;} } } }
◍ 把计算结果丢进浏览器看
TSAnalysis 类的可视化不走常规画图，而是把计算结果格式化成字符串，写进数据文件，再拉起系统浏览器加载同目录下的 TSA.htm 来渲染图形。这套路子在《HTML 中的图形与图表》里讲过，本质是用网页代替终端图表对象。 类里的显示方法会把所有要展示的结果一次性攒成一条长线，直接落盘到 TSDat.txt。文件由标准 MQL5 函数创建，实际落在 \MQL5\Files 目录下，随后靠外部系统函数挪到项目目录，再调浏览器打开 TSA.htm 读取该 txt 数据。 关键点在于：显示方法内部调用了系统 DLL 函数，所以终端必须勾选「允许使用外部 DLL」才能跑通 TSAnalysis。外汇与贵金属波动剧烈、杠杆高风险大，这类依赖外部调用的方案在实盘前务必先在策略测试器里验证路径与权限。
谱估值在合成信号与 EURUSD 上的实测表现
TSAexample.mq5 演示了 TSAnalysis 类的极简调用：准备好输入数组后丢给 Calc，跑完记得 delete 释放。脚本开头的 60 点数组只是冒烟测试，真正说明问题的是后面用生成序列做的几组对照。 先喂两条正弦曲线，频率响应图里峰位落在 0.0637 与 0.0712，和真实频率有偏差；单条正弦就不会出现这种乖离，可视为所选谱估计方法的固有效应。往里加单位方差正态随机量，正弦仍清楚；随机幅度翻五倍后明显被糊住。 把序列长度从 400 拉到 800，自回归阶数由 130 升到 145，峰值显著锐化、分辨力提高——这说明样本量直接决定模型阶与谱清晰度。 拿 2009–2010 年 EURUSD 的 D1 收盘做估值，输入 519 个点、模型阶约 135，图上冒出一堆明显峰。但外汇贵金属属高风险品种，峰值可能来自高阶随机分量或序列非定常性，不能单凭一张谱图断定周期分量存在。 别把正态当圣经 加噪实验里单位方差正态随机量只是理想假设，实盘报价分布常有厚尾，直接套用会低估谱峰模糊程度。 拿另一段序列或上移时间框架交叉验证，再试试用差分序列替代原序列，才可能挤掉假周期。下面这段是示例脚本里构造测试信号的骨架代码。
class=class="str">"cmt">//----------------------------------------------------------------------------------- class=class="str">"cmt">// TSAexample.mqh class=class="str">"cmt">// class="num">2011, victorg class=class="str">"cmt">// [MQL5官方文档] class=class="str">"cmt">//----------------------------------------------------------------------------------- class="macro">#class="kw">property copyright "class="num">2011, victorg" class="macro">#class="kw">property link "[MQL5官方文档] class="macro">#include "TSAnalysis.mqh" class=class="str">"cmt">//----------------------------------------------------------------------------------- class=class="str">"cmt">// 脚本程序起始函数 class=class="str">"cmt">//----------------------------------------------------------------------------------- class="type">void OnStart() { class="type">class="kw">double bd[]={class="num">47,class="num">64,class="num">23,class="num">71,class="num">38,class="num">64,class="num">55,class="num">41,class="num">59,class="num">48,class="num">71,class="num">35,class="num">57,class="num">40,class="num">58,class="num">44,class="num">80,class="num">55,class="num">37,class="num">74,class="num">51,class="num">57,class="num">50, class="num">60,class="num">45,class="num">57,class="num">50,class="num">45,class="num">25,class="num">59,class="num">50,class="num">71,class="num">56,class="num">74,class="num">50,class="num">58,class="num">45,class="num">54,class="num">36,class="num">54,class="num">48,class="num">55,class="num">45,class="num">57,class="num">50,class="num">62,class="num">44,class="num">64,class="num">43,class="num">52, class="num">38,class="num">59,class="num">55,class="num">41,class="num">53,class="num">49,class="num">34,class="num">35,class="num">54,class="num">45,class="num">68,class="num">38,class="num">50,class="num">60,class="num">39,class="num">59,class="num">40,class="num">57,class="num">54,class="num">23}; TSAnalysis *tsa=new TSAnalysis; tsa.Calc(bd); class="kw">delete tsa; } class="type">int i,n; class="type">class="kw">double a,x[]; n=class="num">400; ArrayResize(x,n); a=class="num">2*M_PI;