4. 基于马尔可夫状态切换的价差策略

任务目标

本案例实现一个基于马尔可夫状态切换的价差交易策略：将市场划分为低波动（状态0）和高波动（状态1）两种隐状态，利用价差观测序列估计状态转移概率与观测参数，并通过维特比算法在线预测当前最可能状态，据此动态调整价差交易的入场阈值与仓位。

环境依赖

需要 numpy（矩阵运算）与 scipy（正态分布计算）。安装命令：

pip install numpy scipy

核心实现

以下代码实现完整的马尔可夫状态切换模型：generate_data 生成可复现的模拟价差序列（含真实状态标签用于验证），HMMSwitching 类通过EM算法（Baum-Welch）估计参数，viterbi 方法输出每个时刻的最可能状态序列。策略逻辑根据预测状态调整交易阈值：低波动时使用较紧的入场阈值（0.5倍标准差），高波动时放宽至1.5倍标准差。

import numpy as np
from scipy.stats import norm

class HMMSwitching:
    def __init__(self, n_states=2):
        self.n = n_states
        self.A = np.array([[0.95, 0.05], [0.10, 0.90]])  # 状态转移矩阵
        self.mu = np.array([0.0, 0.0])                    # 各状态均值
        self.sigma = np.array([0.5, 2.0])                 # 各状态标准差
        self.pi = np.array([0.8, 0.2])                    # 初始状态分布

    def _e_step(self, obs):
        T = len(obs)
        # 前向算法
        alpha = np.zeros((T, self.n))
        alpha[0] = self.pi * norm.pdf(obs[0], self.mu, self.sigma)
        alpha[0] /= alpha[0].sum()
        for t in range(1, T):
            for j in range(self.n):
                alpha[t, j] = norm.pdf(obs[t], self.mu[j], self.sigma[j]) * \
                              np.sum(alpha[t-1] * self.A[:, j])
            alpha[t] /= alpha[t].sum()
        # 后向算法
        beta = np.ones((T, self.n))
        for t in range(T-2, -1, -1):
            for i in range(self.n):
                beta[t, i] = np.sum(self.A[i, :] * norm.pdf(obs[t+1], self.mu, self.sigma) * beta[t+1])
            beta[t] /= beta[t].sum()
        # 后验概率 gamma 与 xi
        gamma = alpha * beta
        gamma /= gamma.sum(axis=1, keepdims=True)
        xi = np.zeros((T-1, self.n, self.n))
        for t in range(T-1):
            for i in range(self.n):
                for j in range(self.n):
                    xi[t, i, j] = alpha[t, i] * self.A[i, j] * \
                                  norm.pdf(obs[t+1], self.mu[j], self.sigma[j]) * beta[t+1, j]
            xi[t] /= xi[t].sum()
        return gamma, xi

    def _m_step(self, obs, gamma, xi):
        self.pi = gamma[0]
        self.A = xi.sum(axis=0) / xi.sum(axis=0).sum(axis=1, keepdims=True)
        for k in range(self.n):
            w = gamma[:, k]
            self.mu[k] = np.sum(w * obs) / np.sum(w)
            self.sigma[k] = np.sqrt(np.sum(w * (obs - self.mu[k])**2) / np.sum(w))
            self.sigma[k] = max(self.sigma[k], 1e-4)

    def fit(self, obs, max_iter=50, tol=1e-4):
        prev_loglik = -np.inf
        for _ in range(max_iter):
            gamma, xi = self._e_step(obs)
            self._m_step(obs, gamma, xi)
            # 计算对数似然（用前向概率近似）
            alpha = np.zeros((len(obs), self.n))
            alpha[0] = self.pi * norm.pdf(obs[0], self.mu, self.sigma)
            alpha[0] /= alpha[0].sum()
            for t in range(1, len(obs)):
                for j in range(self.n):
                    alpha[t, j] = norm.pdf(obs[t], self.mu[j], self.sigma[j]) * \
                                  np.sum(alpha[t-1] * self.A[:, j])
                alpha[t] /= alpha[t].sum()
            loglik = np.sum(np.log(alpha[-1].sum() + 1e-12))
            if abs(loglik - prev_loglik) < tol:
                break
            prev_loglik = loglik

    def viterbi(self, obs):
        T = len(obs)
        delta = np.zeros((T, self.n))
        psi = np.zeros((T, self.n), dtype=int)
        delta[0] = self.pi * norm.pdf(obs[0], self.mu, self.sigma)
        for t in range(1, T):
            for j in range(self.n):
                trans = delta[t-1] * self.A[:, j]
                psi[t, j] = np.argmax(trans)
                delta[t, j] = trans[psi[t, j]] * norm.pdf(obs[t], self.mu[j], self.sigma[j])
        states = np.zeros(T, dtype=int)
        states[-1] = np.argmax(delta[-1])
        for t in range(T-2, -1, -1):
            states[t] = psi[t+1, states[t+1]]
        return states

def generate_data(n=500, seed=42):
    """生成模拟价差数据：状态0为低波动，状态1为高波动"""
    rng = np.random.default_rng(seed)
    true_states = np.zeros(n, dtype=int)
    obs = np.zeros(n)
    state = 0
    for t in range(n):
        true_states[t] = state
        obs[t] = rng.normal(0.0, 0.5 if state == 0 else 2.0)
        # 状态切换：低波动→高波动概率0.03，高波动→低波动概率0.10
        if state == 0 and rng.random() < 0.03:
            state = 1
        elif state == 1 and rng.random() < 0.10:
            state = 0
    return obs, true_states

# ========== 策略执行与验证 ==========
obs, true_states = generate_data(500)
model = HMMSwitching()
model.fit(obs, max_iter=30)
pred_states = model.viterbi(obs)

# 策略：状态0阈值0.5σ，状态1阈值1.5σ（σ取全序列标准差）
sigma_global = np.std(obs)
thresholds = np.where(pred_states == 0, 0.5 * sigma_global, 1.5 * sigma_global)
spread = obs  # 价差序列即obs

# 生成交易信号：价差超过阈值做空价差（卖价差），低于负阈值做多价差
signals = np.zeros(len(spread))
signals[spread > thresholds] = -1   # 价差过高，做空
signals[spread < -thresholds] = 1   # 价差过低，做多

# 状态预测准确率（用于验证模型有效性）
accuracy = np.mean(pred_states == true_states)
print(f"状态预测准确率: {accuracy:.2%}")
print(f"交易信号数量: {int(np.sum(signals != 0))}")

运行与验证

直接运行上述脚本即可完成训练、预测与信号生成。验证路径如下：

模型收敛性：fit 方法在30次迭代内应收敛（对数似然变化小于 1e-4），若未收敛可增大 max_iter。
状态区分度：训练后打印的 sigma 应明显区分（低波动约0.5，高波动约2.0），且 mu 接近0。
预测准确率：由于模拟数据状态切换概率较低（平均每个状态持续约20~30个时间步），维特比解码的准确率通常高于85%。若准确率低于80%，可检查 generate_data 中的切换概率是否被EM算法正确捕获。
策略信号合理性：输出的信号数量应在合理范围（通常为总样本的5%~20%），且信号多出现在高波动状态（状态1）期间——因为高波动时阈值更宽，但价差波动也更大，实际触发频率需结合具体序列观察。

判定标准：状态预测准确率 > 85% 且模型参数（sigma）能清晰区分两种波动水平，即视为策略状态切换机制有效。

常见问题

sigma 估计值趋同：若两个状态的 sigma 在训练后几乎相等，说明EM陷入局部最优。可重置 model.sigma = np.array([0.5, 2.0]) 后重新 fit，或增加 max_iter。
状态翻转过于频繁：维特比路径中状态每1~2步就切换，说明转移矩阵对角元素被低估。可对 A 添加先验：训练后执行 model.A = 0.9 * model.A + 0.1 * np.eye(2) 平滑。
实盘应用差异：本案例使用模拟数据验证机制，实盘价差序列需先做平稳性检验（如ADF测试），且参数需滚动重估以适应市场变化。

📚 课程目录