用马尔可夫链预测城市人均GDP：状态转移建模实战解析

“把人均GDP数据变成离散档位，然后算一个城市在一个档位待多久、跳到另一个档位的概率，这听起来像是做天气预测，但我拿它跑了一轮城市经济数据，结果比我预想的要有意思得多。”这是我这几个月做区域经济分析时的一个真实感受。如果你也烦恼于传统的线性回归和ARIMA模型解释不了“城市经济水平突然跳档”之类的现象，那这篇文章大概能给你一个新的工具视角。

我最近完成了一个模拟项目：用马尔可夫链预测模型，对某省各城市的人均GDP做状态转移分析。为了不涉及具体地域的敏感性，我统一把目标省份称为“某省”，下辖城市的数据也做了脱敏处理，但整体的建模流程和思路是完整且可复现的。这篇博文会从模型原理、数据预处理、Python实现，到状态数选择、平滑技巧、结果解读，一步一步讲清楚，并把我踩过的坑也一并交代。

1. 为什么我会选马尔可夫链来预测人均GDP：反直觉背后的合理性

1.1 传统回归模型在城市经济预测上的尴尬

很久以前我在做城市经济预测时，第一反应就是上ARIMA或者线性回归。毕竟人均GDP是一个典型的连续时间序列，直接套回归模型似乎顺理成章。但实际跑起来，问题一堆：数据不平稳，需要差分；经济结构突变，导致拟合参数不稳定；城市之间的差异大，单个城市样本量又少，拟合出某个城市的模型后，换一个城市又完全不能套用。

更麻烦的是，传统回归给出的是一个“期望值”，比如“明年某市人均GDP约为6.2万元”，但这个期望值没有给出“增长落在哪个区间”的概率有多大。而领导或业务方更想知道的是“某市一年从中等收入区间跳入高收入区间的可能性有多大”，这种问题回归模型很难直接回答。这时候，马尔可夫链的“状态转移”思路就变得很有吸引力。

1.2 马尔可夫链不要求平稳，也不要求固定趋势

马尔可夫链的核心假设是“下一时刻的状态只取决于当前状态”，也就是所谓的无后效性。这听起来像是一个很强的简化，但在经济领域，很多短期预测恰恰可以用这种简化来抓住主要矛盾——比如某个城市当前处于低收入水平，两年后它升到中等偏上水平的概率，靠历史同类的转移频次就能估算出来。

与传统时间序列模型不同，马尔可夫链不要求数据平稳，也不需要对趋势做显式建模。你只需要把人均GDP映射成有限个离散状态，然后统计状态转移的频数，就能得到一张转移概率矩阵。这个矩阵本身就是一个解释性极强的产出：它告诉你“整体区域的经济发展通道”长什么样。

我在这个模拟项目中，把某省的城市人均GDP按全省年份数据的四分位数切分成四个状态：低、中低、中高、高。这样一来，每个城市在某一年都有一个状态标签，相邻两年的状态变化就构成一条转移记录。整个省份几十个城市、近二十年数据，累计拿到了数百条转移记录，足以支撑一个稳定的转移概率矩阵估计。

1.3 状态跃迁比数值本身更有决策价值

从业务角度看，城市经济部门关心的问题往往是“能不能从当前档位升到更高档位”。例如一个经济中等的城市，想要预测未来五年内进入高收入行列的可能性，或者一个经济较高的城市，评估它“倒退到中等”的风险。这些问题的本质，就是对“状态空间”的转移概率进行计算，而马尔可夫链恰好是干这个的。

在项目里，我进一步把预测分为两个层次：第一步是用整体转移矩阵预测单个城市在下一期的状态分布；第二步是计算多步转移后，看每个城市进入高状态的概率。这一步在决策上特别有用，因为它给出了一个概率区间，而不是一个看似精确的点值。比如某市当前处于中高状态，通过矩阵连乘，我可以告诉团队：它五年后仍然保持在中高及以上的概率约是78%，而回落到中低及以下的概率只有6%左右——这比一个回归点预测要直观得多。

2. 马尔可夫链预测模型的数学原理：从状态划分到稳态分布，我尽量说人话

2.1 状态空间的定义与档位划分

用马尔可夫链做经济预测，第一步就是定义“状态”。在本项目里，状态就是人均GDP的档位。我没有直接使用原始数值，而是把全省历年人均GDP取值从小到大排列，按四分位数切分。注意，这里的“全省历年”是混合城市、混合年份的数据，这样做的好处是保证不同年份、不同城市的档位具有可比性，坏处是如果某一年全省经济水平大幅增长，同一个状态对应的原始数值可能偏高，但这不影响状态转移的统计逻辑。

具体到模拟数据，我生成了某省20个城市过去20年的人均GDP序列，范围大致从2万元到15万元。取四分位数后，四个状态的门槛大约是4.2万、6.5万和9.4万。低于4.2万的记为状态1（低），4.2万到6.5万记为状态2（中低），6.5万到9.4万记为状态3（中高），高于9.4万记为状态4（高）。这种离散化处理方式损失了一部分数值精度，但换来了概率建模的便利，对经济预测场景来说，档位比精确数值更容易形成判断。

2.2 转移概率矩阵的估计：数数就行，但要解决零频问题

转移概率矩阵P是一个4×4的矩阵，P[i, j]表示从状态i转移到状态j的概率，估计方式就是统计历史转移记录中的频次比例。假设我把所有城市所有相邻年份的状态变化都统计出来，得到频数矩阵N，其中N[i, j]是从状态i转移到状态j的次数，那么P[i, j] = N[i, j] / sum(N[i, :])。

在项目实际运行中，我发现一个常见坑：如果某个状态在城市间极少出现，比如几乎没有城市处于最低状态，那么该行对应的样本量会很小，估计出的概率会很不稳定。这时候需要做平滑处理，比较常用的是拉普拉斯平滑，也就是在计数上加一个小的常数α（比如0.5），避免零概率的出现。这点我在后面专门讲。

2.3 多步预测与稳态分布：矩阵连乘和特征向量

有了单步转移概率矩阵之后，做未来多期预测并不复杂。如果当前状态分布是一个行向量π0（在某个城市，π0是一个one-hot向量，代表它当前处于某个状态），那么t期之后的状态分布就是πt = π0 × P^t。这里P^t就是转移矩阵的t次幂，可以用矩阵乘法连续计算。

稳态分布更有意思。当t足够大时，πt会收敛到一个与初始状态无关的分布，这个分布就是马尔可夫链的稳态分布。在极端情况下，它可以反映区域经济的长期平衡状态——比如在没有任何外力干预时，某省各城市人均GDP的档位最终会以何种概率分布。我算了一下模拟数据，稳态分布大致是状态1占8%，状态2占22%，状态3占35%，状态4占35%，说明长期来看，大多数城市会停留在中高或高状态。

不过要提醒一句，稳态分布是基于现有转移规律外推的结果，如果未来出现政策激励、产业转移等重大变化，转移概率矩阵本身就可能改变，稳态分布自然也就失真了。所以多步预测我通常只做3到5年，不主张推演太远。

3. 数据准备与状态划分：一份干干净净的转移序列是怎么来的

3.1 数据获取和清洗：没有完整的历史数据，后面全是白做

在模拟项目里，我构造了一套“某省城市人均GDP面板数据”，包含城市名、年份、人均GDP三个字段。虽然数据是模拟生成的，但在处理真实统计数据时，数据清洗的逻辑是一样的。第一步是补齐缺失值。统计年鉴里常有某几年数据空缺的情况，我的处理方法是如果相邻两年都有值，就用线性插值补上；如果首尾缺失，就沿用最近一年的值。缺失值处理一定要先于状态划分，否则会导致状态标签错乱。

第二步是统一计量单位。不同来源的数据可能有的用元、有的用万元，一定要换算成同一单位。我在项目中统一为“万元”，方便做四分位数比较。

第三步是剔除异常年份。比如某市某年GDP数据因为行政区划调整而出现突增或突降，这类异常点会扭曲转移概率，需要在划分状态前单独处理。我的做法是计算每个城市人均GDP的年度环比变化，如果变化超过100%的极值，就认为该点是数据质量问题，直接用前后两年均值替代。

3.2 人均GDP的离散化：不要拍脑袋，先看分布形态

离散化方式直接影响预测效果。我尝试过两种方案：第一种是固定阈值，比如直接设4万元、6万元、9万元为门槛；第二种是全省分位数，比如每25%分位切一个状态。对比下来，分位数的好处是不受个别极值影响，能动态适应整体经济发展水平。

具体操作时，我把全省所有年份、所有城市的人均GDP合并成一个一维数组，然后计算其25%、50%、75%分位数。比如在模拟数据中，这三个分位数分别是4.0万、6.2万、9.0万。于是状态划分规则如下：

• 状态1：< 4.0万元
• 状态2：4.0万 ~ 6.2万元
• 状态3：6.2万 ~ 9.0万元
• 状态4：> 9.0万元

我会检查每个状态对应的样本量，确保每个状态至少占全部样本的10%以上。如果有状态样本量过小，就适当调整分位数比例，比如改成20%、40%、60%。这种动态调整在后续估计转移概率时能显著减小方差。

3.3 从面板数据到转移序列：生成训练样本的关键代码逻辑

把离散化规则应用到每一年的数据后，每个城市每年都会得到一个状态标签。随后我按城市分组，对每个城市按年份排序，然后相邻年份的状态变化就是一条转移记录。以城市A为例，如果它2005年状态是1，2006年状态是2，那么记录就是“1 -> 2”。把所有城市的所有相邻年份转移记录汇总，就形成了频数矩阵的输入。

这里有一个容易被忽略的细节：必须是“同一城市”相邻两年的状态变化，不能跨城市或跨年配对。否则相当于打破了马尔可夫链对对象连续性的前提。我最初在生成转移序列时，曾不小心把不同城市同一年的观测当成前后关系，结果转移矩阵明显失真，后来加了分组逻辑才修正。

为了便于后续建模，我把转移记录整理成一个DataFrame，字段包括：from_state、to_state、转移前城市名、转移前年份。这样我既可以汇总全局转移矩阵，也可以按城市类型或时间段做子集分析，比如只使用最近10年的数据进行预测，和全量20年数据做对比。

4. 建模实操：用Python完整实现马尔可夫链预测的每一步

4.1 环境准备与工具选择

整个项目我用Python实现，核心库只有pandas和numpy。没有必要引入重量级的机器学习框架，因为转移概率矩阵本质上就是计数加归一化。为了可复现，我先构建了一个函数，输入是包含“城市、年份、状态”三列的DataFrame，输出是转移概率矩阵和频数矩阵。

项目文件夹结构很简单：

• data/ 存放原始和中间数据
• src/ 存放转换脚本
• models/ 存放输出结果

如果你直接在Jupyter Notebook里分步跑，也完全够用。

4.2 核心实现：状态划分、转移矩阵、多步预测一次性搞定

状态划分我已经在上一节提过，这里直接给出核心代码。第一步是根据全省列人均GDP计算分位数，然后把每个年份的GDP映射到对应状态：

import pandas as pd import numpy as np def discretize_gdp(gdp_series, quantiles=[0.25, 0.5, 0.75]): """ 根据分位数把GDP划分为4个状态。 返回状态标签（1-4） """ thresholds = gdp_series.quantile(quantiles).values def assign_state(value): if value < thresholds[0]: return 1 elif value < thresholds[1]: return 2 elif value < thresholds[2]: return 3 else: return 4 return gdp_series.apply(assign_state)

调用时，我先合并所有城市所有年份的GDP值来计算全局分位数，再逐城市逐年份打标签：

# 假设df是原始数据：city, year, gdp all_gdp = df["gdp"] df["state"] = discretize_gdp(all_gdp) df.sort_values(["city", "year"], inplace=True)

接下来生成转移序列，并统计频数：

transfer_rows = [] for city, group in df.groupby("city"): group = group.sort_values("year") states = group["state"].values for i in range(len(states) - 1): transfer_rows.append([states[i], states[i+1]]) transfer_df = pd.DataFrame(transfer_rows, columns=["from_state", "to_state"]) # 统计频数矩阵 freq_matrix = transfer_df.groupby(["from_state", "to_state"]).size().unstack(fill_value=0) # 转为概率矩阵（每行归一化） prob_matrix = freq_matrix.div(freq_matrix.sum(axis=1), axis=0)

这里要注意，freq_matrix的列可能不完整。比如某个状态从未出现过，导致某一行没有对应的列，这时需要预先定义好完整的4×4索引，用reindex补齐。我在第一次跑代码时就是因为索引不一致，导致矩阵运算报错，所以建议一开始就设定好state_index = [1, 2, 3, 4]，然后对freq_matrix做reindex。

多步预测部分也很直接，我写了一个函数，输入初始状态分布向量（可以用one-hot），以及步数t，输出每条路径的概率分布：

def predict_distribution(init_state, steps, prob_matrix): """ init_state: 初始状态向量，比如 [0, 1, 0, 0] steps: 预测期数 返回: 第1到第steps期的状态分布 """ prob_matrix = prob_matrix.T # 注意行列对齐 res = [] current = np.array(init_state, dtype=float) for _ in range(steps): current = current @ prob_matrix res.append(current) return res

这里我把prob_matrix转置了一下，因为pandas的groupby得到的是“from行，to列”，而矩阵乘法需要“当前状态向量（行向量）乘以转移矩阵”，所以直接用原矩阵的话，行表示当前状态，列表示转移目标，current @ prob_matrix就是对的，不需要转置。我建议在代码里做好注释说明，避免自己下次看的时候绕晕。

4.3 案例输出：某城市未来五年的状态分布到底怎么读

我用某省的中心城市C做了个演示。该城市当前处于状态3（中高），即初始向量为[0, 0, 1, 0]。代入模型，预测未来1到5年的状态分布如下：

从这个结果能看出来，城市C大概率会继续保持在中高或高状态，且进入高状态的概率逐年上升，但上升幅度在放缓。这个“从跃升迅速到逐渐饱和”的过程，就是转移矩阵连乘的结果。你可以把这些概率直接拿去做一个堆积条形图，管理层一看就懂，根本不用解释数学概念。

我在模型里还顺手加了一个“跃迁升级概率”指标，即从状态3升到状态4的累计概率：第1年是0.25，到第5年已经累积到0.41左右，再往下会逼近稳态分布中的0.35附近。这个指标在区域经济评估里很常用。

5. 踩过的坑和优化思路：预测结果不一定稳，但可以想办法更稳

5.1 状态数量的选择：四个状态是我试出来的平衡点

一开始我尝试过6个状态，等于把GDP切成六档，这样信息保留得更多。但问题是，某省城市样本量和时间跨度有限，6个状态下许多转移组合根本没有出现，频数矩阵有一大半是零，估计出的概率方差极大。相反，我后来又试了3个状态，转移矩阵是稳定了，但预测结果过于粗糙，连“中低与中高的区别”都看不出来。

最终我选定4个状态。做这个决定前，我对比了不同状态数量下的BIC统计量（贝叶斯信息准则），4个状态的BIC最低，也就是说在拟合优度和复杂度之间它最平衡。如果你想自行调整，建议先看频数矩阵的稀疏程度：如果非零元素占比小于60%，就考虑减少状态数；如果每个状态样本量都超过50条转移记录，再考虑增加。

5.2 拉普拉斯平滑：样本量不足时的救命稻草

在实际统计中，即使4个状态，也可能会出现“从状态1转移到状态4”的零频现象。这不等于是不可能，只是样本量没覆盖到。如果不做平滑，这个概率就是0，直接导致预测模型中该转移项的引擎完全失效。我做了一个平滑处理：在频数矩阵每个单元格上加一个伪计数α，然后重新归一化。

alpha = 0.5 # 平滑系数 freq_smooth = freq_matrix + alpha prob_smooth = freq_smooth.div(freq_smooth.sum(axis=1), axis=0)

α取多少合适？我试验了0.1、0.5、1.0、3.0。越小越接近原始数据，但零频问题越明显；越大越平滑，但会扭曲真实转移规律。综合来看，0.5在模拟数据上表现最好，既能把零频概率拉到0.05以下，又不至于让低频状态赔率猛涨。

5.3 和回归模型对比：马尔可夫链的优势是“给概率、给通道”

在项目最后，我把马尔可夫链预测结果和一份传统的多元线性回归预测做了对比。回归模型给出某市未来一年人均GDP预测值是6.8万元，而状态划分下，6.5万到9.4万属于状态3。马尔可夫链则告诉我：该市明年位于状态3的概率是63%，位于状态2的概率是10%，还有25%的概率进入状态4。

你可以看到，两者并不矛盾，但信息的维度不同。回归模型适合回答“具体数值大概是多少”，马尔可夫链适合回答“处于某个发展档位的可能性有多大”。在做区域经济比较或者政策决策时，后者其实更有用。因为它能直观反映经济“跨越式发展”的概率，比如从状态2一步跳到状态4，虽然概率只有2%，但这种极端情形势头对决策者极其重要。

还有一个额外收益：转移矩阵本身可以用来做城市间“发展动能”的比较。比如我分别计算了某省北部和南部的城市群转移矩阵，发现南部城市从状态3升到状态4的概率是0.3，而北部只有0.15。这说明南部区域经济活力和联动性更强，可以把更多资源配置给南部潜在跃迁城市。这个视角在传统回归模型里完全看不到。

最后说一个我自己的体会：做这类马尔可夫链预测，最忌讳的就是把模型结果当成必然结论。它的本质是用历史频次推演未来概率，如果未来政策环境出现重大变化，比如大规模产业转移、行政区划调整，转移矩阵必须重新训练。我在项目里保留了按时间段分窗口建模的方法，比如用第一个十年训练、用第二个十年验证，结果显示第二个十年预测准确率明显低于第一个十年，这也提醒我：经济规律本身也在不断演进，定期用滚动窗口更新模型，比一次建完用到底要可靠得多。