基于隐马尔可夫模型的趋势跟踪波动率预测·进阶篇
(2/3)· 从数据获取到模型训练再到MT5集成,手把手拆开HMM波动率过滤链的落地缝隙
用EA把K线收盘数据落盘成CSV
做波动率状态分类(高波动=1,低波动=0)前,先得拿到干净的历史OHLC。MT5终端自带数据多偏实时tick,要跨年度训练样本,最好自己写个EA把收盘价和时间戳按根存下来。 下面这段EA逻辑很直接:每个tick进来时比对已统计K线数,一旦发现新K线成型,就把服务器交易时间和前一根收盘价写进二维数组;脚本退出时通过 CFileCSV 一次性写表头加数据行到CSV。 在策略测试器里跑,文件会落到 /Tester/Agent-sth000 目录。训练样本取2020.1.1–2024.1.1,样本外验证用2024.1.1–2025.1.1,四年窗口足够覆盖几种宏观波动 regime。 代码逐行看: #include <FileCSV.mqh> 引入CSV读写类 int barsTotal = 0 记录已处理K线总数,用于检测新K线 CFileCSV csvFile 声明文件对象 input string fileName = "Name.csv" 外部参数指定输出文件名 string headers[]={"time","close"} 定义CSV两列表头 string data[100000][2] 预开10万行缓存,存时间+收盘 data数组下标由 indexx 自增控制 input bool SaveData=true 总开关,false则不写盘 OnInit 仅返回成功,无初始化负载 OnDeinit 中若SaveData且文件可写,先WriteHeader再WriteLine写全部缓存,最后Close;失败则打印错误 OnTick 用 iBars 取当前K线数,与barsTotal不等即新K线:更新barsTotal,把 TimeTradeServer 转字符串存列0,iClose取前一根收盘转8位精度字符串存列1,indexx++。 实盘或回测前先把 fileName 改成你的品种名,否则多个品种会互盖文件。外汇和贵金属杠杆高、滑点跳空频繁,落盘数据也仅反映历史,不构成方向暗示。
class="macro">#include <FileCSV.mqh> class="type">int barsTotal = class="num">0; CFileCSV csvFile; input class="type">class="kw">string fileName = "Name.csv"; class="type">class="kw">string headers[] = { "time", "close" }; class="type">class="kw">string data[class="num">100000][class="num">2]; class="type">int indexx = class="num">0; vector xx; input class="type">bool SaveData = true; class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Expert initialization function | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">int OnInit() {class=class="str">"cmt">//Initialize model class="kw">return(INIT_SUCCEEDED); } class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Expert deinitialization function | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">void OnDeinit(const class="type">int reason) { if (!SaveData) class="kw">return; if(csvFile.Open(fileName, FILE_WRITE|FILE_ANSI)) { class=class="str">"cmt">//Write the header csvFile.WriteHeader(headers); class=class="str">"cmt">//Write data rows csvFile.WriteLine(data); class=class="str">"cmt">//Close the file csvFile.Close(); } else { Print("File opening error!"); } } class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Expert tick function | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">void OnTick() { class="type">int bars = iBars(_Symbol,PERIOD_CURRENT); if (barsTotal!= bars){ barsTotal = bars; data[indexx][class="num">0] =(class="type">class="kw">string)TimeTradeServer() ; data[indexx][class="num">1] = DoubleToString(iClose(_Symbol,PERIOD_CURRENT,class="num">1), class="num">8); indexx++; } }
「把波动率喂给隐马尔可夫模型」
原始 CSV 里时间戳和收盘价是分号混在一列里的字符串,先拆开再转类型,否则后面算波动率会直接报错。下面这段 Python 是最小可跑的预处理与建模流程,复制去本地环境改个文件路径就能用。 ADF 检验在样本上给出的统计量是 -13.12055,p-value 约 1.57e-24,观测数 23516,结论倾向序列平稳。虽然 HMM 本身带滚动更新不强制要平稳数据,但用平稳输入聚类效果大概率更干净。 波动率用 50 根 K 线的滚动标准差,再走 StandardScaler 压成 N(0,1)。注意 scaler 的 mean 和 scale 必须存下来——后面在 MT5 里做实时推断时要用同一组 μ 和 σ 反标准化,不然状态切分会偏。 模型用 GaussianHMM、n_components=2、n_iter=100,拟合后对每个点预测隐藏状态,聚类分布看起来合理,仅有微小误差。最后把转移矩阵、均值、方差和 scaler 参数格式化成 MQL5 全局变量常量,塞进 JSON 头文件,直接粘进 EA 代码当矩阵初值。 外汇和贵金属波动受杠杆与事件驱动影响大,这套训练结果只代表历史样本特征,实盘迁移须用小布盯盘回测验证,高风险。
class="kw">import pandas as pd data = pd.read_csv("XAU_test.csv",sep=";") data = data.dropna() data["close"] = data["close"].astype(class="type">class="kw">float) data[&class="macro">#x27;time&class="macro">#x27;] = pd.to_datetime(data[&class="macro">#x27;time&class="macro">#x27;]) data.set_index(&class="macro">#x27;time&class="macro">#x27;, inplace=True) data[&class="macro">#x27;volatility&class="macro">#x27;] = data[&class="macro">#x27;returns&class="macro">#x27;].rolling(window=class="num">50).std() Augmented Dickey-Fuller Test: Volatility ADF Statistic: -class="num">13.120552520156329 p-value: class="num">1.5664189630119278e-24 # Lags Used: class="num">46 Number of Observations Used: class="num">23516 => The series is likely stationary. from sklearn.preprocessing class="kw">import StandardScaler scaler = StandardScaler() scaled_volatility = scaler.fit_transform(data[[&class="macro">#x27;volatility&class="macro">#x27;]]) scaled_volatility = scaled_volatility.reshape(-class="num">1, class="num">1) scaler_mean = scaler.mean_[class="num">0] # Mean of the volatility feature scaler_std = scaler.scale_[class="num">0] # Standard deviation of the volatility feature from hmmlearn class="kw">import hmm class="kw">import numpy as np # Define the number of hidden states n_states = class="num">2 # Initialize the Gaussian HMM model = hmm.GaussianHMM(n_components=n_states, covariance_type="full", n_iter=class="num">100, random_state=class="num">42, verbose = True) # Fit the model to the scaled volatility data model.fit(scaled_volatility) # Predict the hidden states hidden_states = model.predict(scaled_volatility) # Add the hidden states to your dataframe data[&class="macro">#x27;hidden_state&class="macro">#x27;] = hidden_states plt.figure(figsize=(class="num">14, class="num">6)) for state in range(n_states): state_mask = data[&class="macro">#x27;hidden_state&class="macro">#x27;] == state plt.plot(data.index[state_mask], data[&class="macro">#x27;volatility&class="macro">#x27;][state_mask], &class="macro">#x27;o&class="macro">#x27;, label=f&class="macro">#x27;State {state}&class="macro">#x27;) plt.title(&class="macro">#x27;Hidden States and Rolling Volatility&class="macro">#x27;) plt.xlabel(&class="macro">#x27;Time&class="macro">#x27;) plt.ylabel(&class="macro">#x27;Volatility&class="macro">#x27;) plt.legend() plt.show() class="kw">import json # Your HMM model parameters transition_matrix = model.transmat_ means = model.means_ covars = model.covars_ # Construct the data in the required format data = { "A": [ [transition_matrix[class="num">0, class="num">0], transition_matrix[class="num">0, class="num">1]], [transition_matrix[class="num">1, class="num">0], transition_matrix[class="num">1, class="num">1]] ], "mu": [means[class="num">0, class="num">0], means[class="num">1, class="num">0]], "sigma_sq": [covars[class="num">0, class="num">0], covars[class="num">1, class="num">0]], "scaler_mean": scaler_mean, "scaler_std": scaler_std } # Create the output content in the desired format output_str = """ const class="type">class="kw">double A[class="num">2][class="num">2] = { {%.16f, %.16f}, {%.16f, %.16f} }; const class="type">class="kw">double mu[class="num">2] = {%.16f, %.16f}; const class="type">class="kw">double sigma_sq[class="num">2] = {%.16f, %.16f}; const class="type">class="kw">double scaler_mean = %.16f; const class="type">class="kw">double scaler_std = %.16f; """ % ( data["A"][class="num">0][class="num">0], data["A"][class="num">0][class="num">1], data["A"][class="num">1][class="num">0], data["A"][class="num">1][class="num">1],
◍ 把 Python 导出的隐马尔可夫参数落进 MT5 头文件
上面这段是 Python 侧把训练好的 HMM 参数写出到 model_parameters.h 的收尾动作,以及 MT5 端直接 include 的常量定义。Python 用 with open 把 mu、sigma_sq、scaler_mean、scaler_std 拼成字符串写盘,终端打印确认保存成功。 MT5 这边用 const 双精度数组接住转移矩阵 A[2][2],其中 A[0][0]=0.9941485184089348 表示状态 0 自转移概率极高,A[1][1]=0.9876122774141759 同理,说明两状态都偏粘滞。mu 与 sigma_sq 分别对应两状态的均值与方差,sigma_sq[1]=1.4515804806463273 明显大于状态 0 的 0.1073520489683213,意味着状态 1 的波动发散程度高一个数量级。 scaler_mean=0.0018685496675093、scaler_std=0.0008350190448735 是归一化用的平移与缩放量,回测前必须拿原始收益率减 mean 再除 std,否则状态判别会整体偏移。外汇与贵金属杠杆高、滑点随机,直接套用未缩放参数可能触发误判,开 MT5 把这份 .h 挂到 EA 里先跑历史 tick 验证状态切换频率。
const class="type">class="kw">double A[class="num">2][class="num">2] = { {class="num">0.9941485184089348, class="num">0.0058514815910651}, {class="num">0.0123877225858242, class="num">0.9876122774141759} }; const class="type">class="kw">double mu[class="num">2] = {-class="num">0.4677410059727503, class="num">0.9797900996225393}; const class="type">class="kw">double sigma_sq[class="num">2] = {class="num">0.1073520489683212, class="num">1.4515804806463273}; const class="type">class="kw">double scaler_mean = class="num">0.0018685496675093; const class="type">class="kw">double scaler_std = class="num">0.0008350190448735;
把波动率塞进维特比的状态机
在 MT5 代码编辑器里接着原策略往下写,先补一个滚动波动率模块。GetVolatility() 做的是把标准化后的价格百分比变动求标准差,再塞进长度为 50 的环形缓冲,供后面的隐马尔可夫推断使用;ComputeStdDev() 则单纯算指定数组的标准差,不碰全局状态。 真正的核心是两个推断函数。Viterbi() 用动态规划在给定观测序列 obs[] 时反推最可能的隐含状态链:T1[s][t] 存时刻 t 落在状态 s 的最优路径对数概率,T2[s][t] 存最优前驱;t=0 用先验 π[s] 乘首个发射概率初始化,t=1..49 每步把转移概率 A 转对数空间防下溢,最后从 t=49 最大概率状态沿 T2 回溯填进 states[]。 PredictCurrentState() 只是个薄封装:开 50 长 states[],调 Viterbi(obs) 算完返回 states[49] 当作当前状态估计。外汇与贵金属波动剧烈,这种状态判别只反映概率倾向,实盘前务必在策略测试器跑多周期验证。 下面这段是 GetVolatility 的主体,注意它拿最近两根收盘做百分比变动,缓冲满 50 笔后才算 stddev 并做 (x-mean)/std 的缩放。 别拿未填满的缓冲当信号 percent_change_filled 变 true 之前,scaled_stddev 根本没更新,obs[] 前段会是 0 值。直接拿去跑 Viterbi 会把早期状态推断带偏,建议加一段『缓冲未满则返回上一帧状态』的保护。
class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Get volatility Function | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">void GetVolatility(){ class=class="str">"cmt">// Step class="num">1: Get the last two close prices to compute the latest percent change class="type">class="kw">double close_prices[class="num">2]; class="type">int copied = CopyClose(_Symbol, PERIOD_CURRENT, class="num">1, class="num">2, close_prices); if(copied != class="num">2){ Print("Failed to copy close prices. Copied: ", copied); class="kw">return; } class=class="str">"cmt">// Step class="num">2: Compute the latest percent change class="type">class="kw">double latest_close = close_prices[class="num">0]; class="type">class="kw">double previous_close = close_prices[class="num">1]; class="type">class="kw">double percent_change = class="num">0.0; if(previous_close != class="num">0){ percent_change = (latest_close - previous_close) / previous_close; } else{ Print("Previous close price is zero. Percent change set to class="num">0."); } class=class="str">"cmt">// Step class="num">3: Update the percent_changes buffer percent_changes[percent_change_index] = percent_change; percent_change_index++; if(percent_change_index >= class="num">50){ percent_change_index = class="num">0; percent_change_filled = true; } class=class="str">"cmt">// Step class="num">4: Once the buffer is filled, compute the rolling std dev if(percent_change_filled){ class="type">class="kw">double current_stddev = ComputeStdDev(percent_changes, class="num">50); class=class="str">"cmt">// Step class="num">5: Scale the std dev class="type">class="kw">double scaled_stddev = (current_stddev - scaler_mean) / scaler_std; class=class="str">"cmt">// Step class="num">6: Update the volatility array(ring buffer for Viterbi) class=class="str">"cmt">// Shift the volatility array to make room for the new std dev for(class="type">int i = class="num">0; i < class="num">49; i++){ volatility[i] = volatility[i+class="num">1];
「波动率数组落点与标准差函数」
上一段把缩放后的标准差塞进长度为 50 的 volatility 数组第 49 位,相当于用最新一根棒的离散度刷新隐状态观测序列的末端。 下面这段 ComputeStdDev 是纯样本标准差实现:先判 size<=1 直接返 0,避免除零;随后用一次遍历累加 sum 与 sum_sq,按方差公式 (sum_sq - sum^2/size)/(size-1) 算无偏估计,最后 MathSqrt 开方。
class="type">class="kw">double ComputeStdDev(class="type">class="kw">double &data[], class="type">int size) { if(size <= class="num">1) class="kw">return class="num">0.0; class="type">class="kw">double sum = class="num">0.0; class="type">class="kw">double sum_sq = class="num">0.0; for(class="type">int i = class="num">0; i < size; i++) { sum += data[i]; sum_sq += data[i] * data[i]; } class="type">class="kw">double mean = sum / size; class="type">class="kw">double variance = (sum_sq - (sum * sum) / size) / (size - class="num">1); class="kw">return MathSqrt(variance); }
class="type">class="kw">double ComputeStdDev(class="type">class="kw">double &data[], class="type">int size) { if(size <= class="num">1) class="kw">return class="num">0.0; class="type">class="kw">double sum = class="num">0.0; class="type">class="kw">double sum_sq = class="num">0.0; for(class="type">int i = class="num">0; i < size; i++) { sum += data[i]; sum_sq += data[i] * data[i]; } class="type">class="kw">double mean = sum / size; class="type">class="kw">double variance = (sum_sq - (sum * sum) / size) / (size - class="num">1); class="kw">return MathSqrt(variance); }
◍ 维特比回溯与当前状态预测的实现落点
前一段把前向概率 T1 和回溯指针 T2 填完之后,真正的隐藏状态序列要靠反向追踪拿回来。代码里先在最后一列(索引 49,对应 50 根 K 线窗口)扫一遍 T1,找出 max_final_prob 对应的 last_state,这一步决定了整条路径的终点归属。 回溯循环从 t=48 倒推到 0,每一帧的 states[t] 直接取 T2[states[t+1]][t+1],也就是下一帧选中状态所指向的前驱。这样一条长度为 50 的状态链就还原了,数组用 ArrayResize(states, 50) 显式开好,避免 MT5 里越界报毒。 PredictCurrentState 只是把 Viterbi 包了一层:丢进 obs 数组,拿到 states 后直接返回 states[49]。也就是说,你只要维护好最近 50 个观测值(比如标准化后的收盘价变动),调这个函数就能拿到「当前最可能的隐状态」是 0 还是 1。 实盘里这俩状态常被映射成「震荡」与「趋势」倾向,但外汇和贵金属波动受事件驱动,隐状态只是概率归类,不代表后市必沿该结构走,杠杆品种高风险,参数窗口和均值方差要随品种重估。
if(prob > max_prob) { max_prob = prob; max_state = s_prev; } class=class="str">"cmt">// Calculate emission probability with epsilon class="type">class="kw">double emission_prob = (class="num">1.0 / MathSqrt(class="num">2 * M_PI * sigma_sq[s])) * MathExp(-MathPow(obs[t] - mu[s], class="num">2) / (class="num">2 * sigma_sq[s])) + class="num">1e-10; T1[s][t] = max_prob + MathLog(emission_prob); T2[s][t] = max_state; } } class=class="str">"cmt">// Backtrack to find the optimal state sequence class=class="str">"cmt">// Find the state with the highest probability in the last column class="type">class="kw">double max_final_prob = -DBL_MAX; class="type">int last_state = class="num">0; for(class="type">int s = class="num">0; s < class="num">2; s++) { if(T1[s][class="num">49] > max_final_prob) { max_final_prob = T1[s][class="num">49]; last_state = s; } } class=class="str">"cmt">// Initialize the states array ArrayResize(states, class="num">50); states[class="num">49] = last_state; class=class="str">"cmt">// Backtrack for(class="type">int t = class="num">48; t >= class="num">0; t--) { states[t] = T2[states[t+class="num">1]][t+class="num">1]; } class="kw">return class="num">0; class=class="str">"cmt">// Success } class=class="str">"cmt">//+------------------------------------------------------------------+ class=class="str">"cmt">//| Predict Current Hidden State | class=class="str">"cmt">//+------------------------------------------------------------------+ class="type">int PredictCurrentState(class="type">class="kw">double &obs[]) { class=class="str">"cmt">// Define states array class="type">int states[class="num">50]; class=class="str">"cmt">// Apply Viterbi class="type">int ret = Viterbi(obs, states); if(ret != class="num">0) class="kw">return -class="num">1; class=class="str">"cmt">// Error class=class="str">"cmt">// Return the most probable current state class="kw">return states[class="num">49]; }