用Python实战马尔可夫链：从药厂市场预测到股票分析（附完整代码）

最近和几个做数据分析的朋友聊天，发现一个挺有意思的现象：很多朋友对马尔可夫链这个名字不陌生，知道它是种预测方法，但真要动手写代码把它用在实际项目里，又总觉得隔着一层窗户纸。要么被复杂的数学公式吓退，要么不知道如何把理论转化成可运行的脚本。其实，马尔可夫链的核心思想非常直观——未来只取决于现在，与过去无关。这个看似简单的假设，却能帮我们解决从市场份额波动到股价趋势分析等一系列实际问题。

今天，我就抛开那些令人头疼的推导，直接带大家用Python手把手实现两个完整的案例。一个是经典的药厂市场份额预测，另一个是稍微进阶一点的股票价格状态分析。我们的工具箱很简单，主要就是numpy和pandas，目标是让你看完就能把代码复制过去，改改数据就能用在自己的项目里。无论你是想优化产品的用户留存分析，还是试图理解某些序列数据的转移规律，这篇文章提供的思路和代码框架都能给你带来直接的启发。

1. 马尔可夫链：五分钟理解核心概念与数学骨架

在开始敲代码之前，我们得先确保站在同一个理解层面上。马尔可夫链不是什么玄学，你可以把它想象成一个在不同“状态”之间跳来跳去的系统。比如，一个用户今天使用了App（状态A），明天可能继续使用（状态A），也可能流失（状态B）。关键假设在于，他明天是去是留，只取决于今天用没用，而和上周、上个月是否用过无关。这就是所谓的“无记忆性”或马尔可夫性。

这个过程中，最核心的数学对象是状态转移概率矩阵。它就像一个路由表，清晰地规定了从任何一个状态出发，下一步跳到其他各个状态的可能性有多大。

注意：转移概率矩阵的每一行之和必须严格等于1。因为从任何一个状态出发，下一步必然转移到所有可能状态（包括停留在原状态）中的某一个，这是一个完备事件。

我们来形式化地定义一下。假设系统有3个状态：S1, S2, S3。那么转移概率矩阵 P 就是一个 3x3 的矩阵：

P = [[P11, P12, P13],
     [P21, P22, P23],
     [P31, P32, P33]]

其中，P12 表示当前处于状态 S1 时，下一步转移到状态 S2 的概率。计算这个概率最朴实（也最常用）的方法，就是用频率来近似。

如何从历史数据中计算转移矩阵？

假设我们有一串状态序列数据：[A, B, A, C, B, B, A, C, ...]。计算从状态 A 转移到状态 B 的概率 P_AB，我们只需要：

1. 找出所有出现状态 A 的位置（最后一个位置除外，因为它没有“下一个”）。
2. 统计在这些 A 之后，紧接着出现 B 的次数。
3. P_AB = (A->B 的次数) / (A 出现的总次数 - 1)。

用Python实现这个计算逻辑非常直接，这也是我们后续所有案例的基石。理解了这个，复杂的市场预测和股票分析就都只是这个核心思想在不同数据形态上的应用而已。

2. 案例一：药厂市场份额预测实战

让我们从一个经典的、结构清晰的商业案例开始。假设市场上有三家药厂（A, B, C）生产同一种药品，共有1000家规模相当的采购单位（如医院、药店）。当前，采购A、B、C产品的单位分别有400、300、300家。我们通过市场调研，得到了一个季度内采购单位更换供应商的流动情况表。

我们的目标是：预测未来几个季度后，三家药厂的市场份额将如何变化？

2.1 问题定义与数据准备

首先，我们把问题用数据表示出来。初始市场份额，也就是初始状态概率分布，是一个向量：

initial_distribution = [0.4, 0.3, 0.3] # 对应A, B, C

调研得到的客户流动数据，通常以转移数量表的形式呈现。假设我们得到如下数据（单位：家）：

本期采购	下期采购A	下期采购B	下期采购C
A (400家)	160	120	120
B (300家)	90	180	30
C (300家)	60	90	150

这张表怎么读？第一行：当前季度采购A的400家单位中，下个季度有160家继续采购A，120家转去采购B，120家转去采购C。

我们的第一个任务，就是把这张数量表，转换成马尔可夫链所需的状态转移概率矩阵。

2.2 核心代码：计算转移概率矩阵与预测

下面我们用Python来实现这个过程。我会在代码中加入大量注释，确保每一行都清晰可理解。

import numpy as np
import pandas as pd

# 定义转移数量矩阵，行：本期状态（From），列：下期状态（To）
# 顺序为 [A, B, C]
transition_counts = np.array([
    [160, 120, 120],  # 从A转移到 A, B, C的数量
    [90,  180, 30],   # 从B转移到 A, B, C的数量
    [60,  90,  150]   # 从C转移到 A, B, C的数量
])

# 计算转移概率矩阵：每一行的数量除以该行的总和
# 保持精度，使用float类型
row_sums = transition_counts.sum(axis=1, keepdims=True)
# 避免除零错误，但本例中不会出现
transition_matrix = transition_counts / row_sums

print("状态转移概率矩阵 P:")
print(pd.DataFrame(transition_matrix, index=['A', 'B', 'C'], columns=['A', 'B', 'C']))
print("\n验证每行之和是否为1:")
print(transition_matrix.sum(axis=1))

运行这段代码，你会得到一个如下的矩阵：

状态转移概率矩阵 P:
     A    B    C
A  0.4  0.3  0.3
B  0.3  0.6  0.1
C  0.2  0.3  0.5

这个矩阵就是整个模型的发动机。它告诉我们：

• A厂的客户忠诚度是40%，流失到B和C的各30%。
• B厂的客户忠诚度最高，达60%，流失到A占30%，到C仅10%。
• C厂客户忠诚度50%，流失到A占20%，到B占30%。

接下来进行预测。预测未来第k个时期后的市场份额，本质上就是计算初始分布向量与转移矩阵P的k次幂的乘积。

def predict_market_share(initial_vec, trans_matrix, periods):
    """
    预测未来多个时期后的市场分布
    :param initial_vec: 初始状态概率分布，一维数组
    :param trans_matrix: 状态转移概率矩阵，二维数组
    :param periods: 要预测的未来时期数（列表）
    :return: 字典，键为时期，值为对应的预测分布
    """
    predictions = {}
    current_vec = np.array(initial_vec)
    
    for period in range(max(periods) + 1):
        if period in periods:
            predictions[period] = current_vec.copy()
        # 计算下一个时期的状态分布：当前分布 * 转移矩阵
        current_vec = current_vec @ trans_matrix
    
    return predictions

# 初始分布
initial_dist = np.array([0.4, 0.3, 0.3])
# 想预测第1、2、3、4个季度后的情况
periods_to_predict = [1, 2, 3, 4, 8, 12]

results = predict_market_share(initial_dist, transition_matrix, periods_to_predict)

print("未来市场份额预测:")
for period, dist in results.items():
    print(f"第 {period} 季度后: A={dist[0]:.3f}, B={dist[1]:.3f}, C={dist[2]:.3f}")

2.3 结果分析与业务洞察

运行上面的预测代码，我们可能会得到类似下面的结果：

未来市场份额预测:
第 1 季度后: A=0.310, B=0.390, C=0.300
第 2 季度后: A=0.282, B=0.426, C=0.292
第 3 季度后: A=0.267, B=0.445, C=0.288
第 4 季度后: A=0.259, B=0.456, C=0.285
第 8 季度后: A=0.250, B=0.469, C=0.281
第 12 季度后: A=0.250, B=0.469, C=0.281

看到这里，你能得出什么业务结论？

1. 趋势判断：A厂的市场份额从40%开始持续下滑，B厂份额持续上升，C厂份额小幅下降后趋于稳定。
2. 稳态预测：在大约第8个季度后，市场份额分布基本稳定在 [0.250, 0.469, 0.281] 附近。这个稳定的分布被称为马尔可夫链的平稳分布。这意味着在当前的客户流动规律不变的前提下，长期来看，B厂将占据接近47%的市场，成为主导者。
3. 行动启示：对于A厂和C厂的管理者来说，这个预测是一个强烈的预警。他们需要分析客户流失的原因（是价格、服务还是产品力？），并重点研究B厂高客户留存率（60%）和从A、C厂吸引客户的能力，从而调整自身策略。

这个案例完美展示了如何将一个具体的业务问题（市场份额预测）转化为马尔可夫链问题，并通过不到50行的Python代码获得有指导意义的量化洞察。

3. 案例二：股票价格变动分析

把马尔可夫链用在金融时间序列上，想法很直接：我们不去预测明天的具体股价（那太难了），而是预测明天股价相对于今天，是上涨、下跌还是持平的概率。我们把“涨”、“跌”、“平”定义为三种状态，分析它们之间相互转换的规律。

3.1 数据预处理与状态定义

我们以一支股票的历史日度收盘价数据为例。首先，需要将连续的价格数据，离散化成我们定义的状态序列。

import yfinance as yf # 需要安装：pip install yfinance
import pandas as pd
import numpy as np

# 获取股票数据（这里以苹果为例）
ticker = 'AAPL'
start_date = '2023-01-01'
end_date = '2024-01-01'

stock_data = yf.download(ticker, start=start_date, end=end_date)
prices = stock_data['Close']

# 计算每日收益率
returns = prices.pct_change().dropna()

# 定义状态：根据收益率阈值划分
# 这里是一个示例，阈值可以根据实际情况调整
threshold = 0.005  # 0.5%
states = []
for r in returns:
    if r > threshold:
        states.append('Up')   # 上涨
    elif r < -threshold:
        states.append('Down') # 下跌
    else:
        states.append('Flat') # 持平

# 创建状态序列的DataFrame
state_series = pd.DataFrame({
    'Date': returns.index,
    'Return': returns.values,
    'State': states
})
print(state_series.head(10))

3.2 构建股价变动的马尔可夫链

有了状态序列 ['Up', 'Down', 'Flat', 'Up', 'Flat', ...]，我们就可以像药厂案例一样，计算状态转移概率矩阵。

def build_transition_matrix(state_sequence):
    """
    根据状态序列构建转移概率矩阵
    :param state_sequence: 状态列表，如 ['Up', 'Down', 'Flat', ...]
    :return: 转移概率矩阵（DataFrame）和状态列表
    """
    # 获取所有唯一状态并排序，确保矩阵顺序一致
    unique_states = sorted(list(set(state_sequence)))
    state_index = {s: i for i, s in enumerate(unique_states)}
    n_states = len(unique_states)
    
    # 初始化计数矩阵
    count_matrix = np.zeros((n_states, n_states), dtype=int)
    
    # 遍历序列，统计转移次数
    for i in range(len(state_sequence) - 1):
        current_state = state_sequence[i]
        next_state = state_sequence[i + 1]
        row = state_index[current_state]
        col = state_index[next_state]
        count_matrix[row, col] += 1
    
    # 将计数转换为概率
    prob_matrix = count_matrix.astype(float)
    row_sums = prob_matrix.sum(axis=1, keepdims=True)
    # 处理除零情况（某状态未出现）
    np.divide(prob_matrix, row_sums, out=prob_matrix, where=row_sums!=0)
    
    # 转换为易于阅读的DataFrame
    prob_df = pd.DataFrame(prob_matrix,
                           index=unique_states,
                           columns=unique_states)
    return prob_df, unique_states

# 应用函数
transition_df, state_list = build_transition_matrix(states)
print("股票价格变动状态转移概率矩阵:")
print(transition_df.round(3))

运行后，你可能会得到一个这样的矩阵：

股票价格变动状态转移概率矩阵:
      Down  Flat    Up
Down  0.450  0.300 0.250
Flat  0.250  0.500 0.250
Up    0.200  0.350 0.450

这个矩阵揭示了有趣的规律：

• 持续性：对角线上的值（Down->Down, Flat->Flat, Up->Up）相对较高，说明市场存在一定的“动量”或“反转”惰性，当前状态倾向于持续。
• 对称性？ 观察 Down->Up (0.25) 和 Up->Down (0.20) 的概率，它们不一定相等，这可以用于分析上涨后下跌与下跌后上涨的难易程度。
• 平稳状态：Flat状态转移到Up和Down的概率都是0.25，自身保持的概率是0.5，说明横盘整理后，选择方向的概率大致相等。

3.3 基于当前状态的短期预测

假设今天是交易日结束，我们计算出今日状态为 Up。我们可以直接利用转移矩阵的第一行（对应Up状态）来预测明天的状态概率。

# 假设今日状态为 'Up'
current_state = 'Up'
print(f"\n基于今日状态为 '{current_state}' 的明日状态概率预测:")
tomorrow_probs = transition_df.loc[current_state]
print(tomorrow_probs)

# 找出最可能的下一个状态
most_likely_next = tomorrow_probs.idxmax()
print(f"最可能出现的明日状态是: '{most_likely_next}' (概率：{tomorrow_probs.max():.2%})")

更进一步，我们可以计算未来多日的状态概率分布，就像药厂案例那样，通过连续乘以转移矩阵来实现。

def predict_future_states(start_state, trans_matrix_df, steps):
    """
    预测从某个初始状态开始，未来多步的状态概率分布
    :param start_state: 初始状态（字符串）
    :param trans_matrix_df: 转移概率矩阵DataFrame
    :param steps: 预测步数列表
    :return: 包含各步预测分布的DataFrame
    """
    states = trans_matrix_df.index.tolist()
    # 创建初始分布向量（one-hot编码）
    current_vec = np.zeros(len(states))
    current_vec[states.index(start_state)] = 1.0
    
    trans_matrix = trans_matrix_df.values
    results = []
    
    for step in range(max(steps) + 1):
        if step in steps:
            results.append(current_vec.copy())
        current_vec = current_vec @ trans_matrix
    
    result_df = pd.DataFrame(results, index=[f'Step_{s}' for s in steps], columns=states)
    return result_df

# 预测如果今天大涨('Up')，未来1、2、3天后的状态概率
future_pred = predict_future_states('Up', transition_df, [1, 2, 3, 5])
print("\n未来多步状态概率预测（初始状态为'Up'）:")
print(future_pred.round(3))

这个模型能做什么？

1. 量化市场情绪转换：提供一种量化“上涨动能能否持续”或“下跌后反弹概率”的视角。
2. 策略回测辅助：可以基于“当出现X状态序列后，下个状态为Y的概率显著提升”来设计简单的交易信号。
3. 风险管理：评估资产价格陷入连续下跌（Down->Down）状态的概率。

提示：金融数据噪声极大，直接用历史数据算出的转移矩阵预测未来，效果往往有限。更稳健的做法是使用滚动窗口计算动态的转移矩阵，或者结合其他因子（如成交量、市场情绪指数）来构建更复杂的状态（例如“放量上涨”、“缩量下跌”）。

4. 模型评估、局限性与高级应用思路

通过两个案例，我们已经掌握了马尔可夫链建模的基本流程：定义状态 -> 从数据计算转移矩阵 -> 进行预测。但在实际应用中，直接套用可能会掉进坑里。这一节我们聊聊怎么评估模型好坏，它有哪些天生的局限，以及如何把它用得更高阶。

4.1 如何验证你的马尔可夫模型？

模型建好了，但它靠谱吗？这里有几个实用的检查方法：

• 历史数据拟合度：将模型预测的历史状态分布与实际历史分布进行比较。可以使用卡方检验等统计方法评估差异是否显著。
• 样本外预测：将数据分为训练集和测试集。用训练集计算转移矩阵，去预测测试集的状态序列，计算预测准确率。
• 稳定性检验：计算转移矩阵的幂次 P^n。当n增大时，如果矩阵的行向量趋于一致，说明存在平稳分布，模型在长期预测上可能有一定意义。反之，如果剧烈波动，则长期预测不可信。

# 简单示例：检查多步转移后矩阵是否稳定（收敛）
def check_matrix_convergence(P, max_power=20, tolerance=1e-6):
    """
    检查转移矩阵P的幂次是否收敛
    """
    prev_Pn = P
    for n in range(2, max_power+1):
        Pn = np.linalg.matrix_power(P, n)
        # 检查连续两次幂的矩阵是否几乎相同
        if np.allclose(Pn, prev_Pn, atol=tolerance):
            print(f"转移矩阵在 n={n} 次幂后趋于稳定。")
            print(f"稳定后的矩阵（近似）:")
            print(pd.DataFrame(Pn, index=['A','B','C'], columns=['A','B','C']).round(4))
            return Pn
        prev_Pn = Pn
    print(f"在 {max_power} 次幂内未观察到明显收敛。")
    return None

# 使用药厂案例的矩阵
P_pharma = np.array([[0.4, 0.3, 0.3],
                      [0.3, 0.6, 0.1],
                      [0.2, 0.3, 0.5]])
check_matrix_convergence(P_pharma)

4.2 马尔可夫链的局限性

没有完美的模型，马尔可夫链的局限性也很明显：

1. 马尔可夫性假设过强：现实世界中，很多过程都有“记忆”。用户明天的购买决策，很可能受他过去一周的体验影响，而不仅仅是今天。股票价格更是受长期趋势、基本面等多种历史因素驱动。
2. 状态定义的主观性：如何划分“上涨”、“下跌”、“持平”？阈值设为0.5%还是1%？不同的定义会得到完全不同的转移矩阵和结论。在药厂案例中，如果客户流动数据是按月统计而非按季度，结果也会不同。
3. 静态转移矩阵：我们假设转移概率 P 是固定不变的。但在现实中，市场环境、公司策略、宏观经济都在变化，转移概率本身应该是时变的。
4. 不适用于长期预测：正如原始资料中指出的，马尔可夫链在短期预测中可能有效，但对于长期预测，由于误差累积和模型假设的偏离，其准确性会迅速下降。它更适合分析趋势和稳态，而非精确预测遥远未来的具体点位。

4.3 进阶应用方向与技巧

知其局限，方能突破。下面是一些让马尔可夫模型更强大的思路：

• 引入高阶马尔可夫链：如果觉得一阶（只依赖前一个状态）不够，可以尝试二阶或更高阶。例如，定义状态为“连续两天上涨后”，预测第三天的状态。这相当于扩大了状态空间（从Up/Down变为(Up,Up), (Up,Down)等），能捕捉更长的历史依赖。
• 结合隐马尔可夫模型（HMM）：这是马尔可夫链的超级升级版。它假设我们观察到的数据（如股价）是由一个我们看不见的、遵循马尔可夫链的“隐藏状态”（如“牛市”、“熊市”、“震荡市”）生成的。HMM可以用来推断这些隐藏的市场状态，是非常强大的时间序列分析工具。Python的hmmlearn库可以方便地实现。
• 动态转移矩阵：使用滚动时间窗口重新计算转移矩阵，让模型能够适应变化。例如，用过去60天的数据计算当前的转移矩阵，用于预测明天。
• 与机器学习模型融合：将马尔可夫链的状态或转移概率作为特征，输入到梯度提升树（如XGBoost）或神经网络中，与其他技术指标一起进行预测。

我在一个用户留存分析的项目中就用过动态转移矩阵的思路。我们不再简单定义“活跃”和“流失”两个状态，而是细分为“新用户”、“低频活跃”、“高频活跃”、“风险用户”、“流失用户”五个状态，并每周更新转移矩阵。这样，我们不仅能预测整体的流失率，还能清晰地看到用户在不同生命周期阶段间的流转路径，精准定位流失漏斗的瓶颈环节。这套方法比单纯看日活、月活数字要深刻得多。