简介一份面向经济学研究人员、政策制定者及高校师生的论文复现资料聚焦1988—2016年中国省级经济周期的一致波动、区域协同与异质分化运用层级动态因子模型分析全国、区域与省级三层传导机制并据此提出区域联动协调的政策建议。共1个PDF文件大小921KB除完整论文外还配有Python代码及逐步解释覆盖数据预处理、共同因子与区域因子提取、格兰杰因果检验、K-means聚类和结果可视化等环节方便读者对照复现。论文揭示全国共同因子主导多数省份的中高速、低波动新增长模式同时呈现东部引领—中西部跟随—东北异化的三元格局东部、中部、西部存在“俱乐部收敛”而东北三省走势分化作者还提出“三力均衡指数”以量化区域协同状态。目前已有68人学习适合需要掌握省级经济周期量化分析方法、开展区域经济实证研究或政策评估的读者直接参考代码逻辑与结果解释。1. 层级动态因子模型论文复现从数据构造到政策结论的完整链路这篇文章要拆的资源是一份可运行的论文复现代码主题是《中国省级经济周期的一致波动、区域协同与异质分化》。你拿到手的是完整的 Python 实现包括模拟数据生成、层级动态因子模型Hierarchical Dynamic Factor Model构建、因子贡献度分解、格兰杰因果检验、K-means 俱乐部收敛分析以及区域同步性的滚动相关性计算。这套代码覆盖了从 1988 到 2016 年中国 31 个省级行政区经济周期波动研究的主要实证链路核心结论指向一个判断全国共同因子主导多数省份的经济周期波动区域层面呈显「东部引领、中西跟随、东北独立」的三元结构而省级层面存在「俱乐部收敛」但东北三省异化。对于做区域经济学量化研究、需要复现权威论文实证流程、或者正在写政策建议报告的人来说这份代码直接省掉了从论文到可执行代码的翻译成本。我拆完一遍后确认它不只是把 statsmodels 的 DynamicFactor 跑通那么简单里面数据生成逻辑、因子解释度算法、格兰杰检验的滞后设定每一步都有值得细抠的细节。2. 模拟数据生成与预处理为什么用合成数据权重参数怎么设2.1 数据生成的分层逻辑论文原文没有提供现成的省级季度 GDP 数据文件所以代码的第一步是用 numpy 构造一个符合层级因子结构的合成数据集。这个设计不是偷懒而是为了在拿到真实数据前先把模型链路跑通——如果连已知因子结构的数据都拟合不出来那换到真实数据上只会更难排查问题。np.random.seed(123) years range(1988, 2017) provinces [北京, 天津, 河北, 山西, 内蒙古, 辽宁, 吉林, 黑龙江, 上海, 江苏, 浙江, 安徽, 福建, 江西, 山东, 河南, 湖北, 湖南, 广东, 广西, 海南, 重庆, 四川, 贵州, 云南, 西藏, 陕西, 甘肃, 青海, 宁夏, 新疆] data pd.DataFrame(indexpd.date_range(start1988-01-01, end2016-12-31, freqQ)) national_factor np.sin(np.linspace(0, 10*np.pi, len(data))) np.random.normal(0, 0.2, len(data)) east_factor 0.8*national_factor np.random.normal(0, 0.3, len(data)) central_factor 0.6*national_factor np.random.normal(0, 0.4, len(data)) west_factor 0.5*national_factor np.random.normal(0, 0.5, len(data)) northeast_factor 0.3*national_factor np.random.normal(0, 0.6, len(data))这段代码的关键在于np.sin(np.linspace(0, 10*np.pi, len(data)))这个表达式。linspace(0, 10*pi, 116)生成 116 个等间距的点1988Q1 到 2016Q4 共 116 个季度正弦函数保证序列是周期波动的5 个完整周期对应约 5.6 年一个经济周期这和现实中的朱格拉周期尺度接近。加上np.random.normal(0, 0.2, len(data))的噪声让因子不会太完美符合真实宏观数据的粗糙感。各区域因子的构造使用「全国因子乘以系数加噪声」的形式这个系数代表区域因子对全国因子的敏感度——东部 0.8 最高、东北 0.3 最低暗示东北与全国周期联动弱。噪声标准差从东部的 0.3 递增到东北的 0.6进一步拉开区域差异。这样生成的数据在后续模型拟合时如果能够识别出全国因子主导、东北独立的结构就说明建模链路没有系统性偏差。2.2 省级数据合成与标准化陷阱for i, prov in enumerate(provinces): if i in [0,1,8,9,10,13,14,18,20]: region_factor east_factor elif i in [3,4,11,15,16,17]: region_factor central_factor elif i in [21,22,23,24,25,26,27,28,29,30]: region_factor west_factor else: region_factor northeast_factor prov_factor np.random.normal(0, 0.2, len(data)) data[prov] 0.6*national_factor 0.3*region_factor 0.1*prov_factor np.random.normal(0, 0.1, len(data)) scaler StandardScaler() data_scaled pd.DataFrame(scaler.fit_transform(data), columnsdata.columns, indexdata.index)省级 GDP 的合成公式是0.6*全国因子 0.3*区域因子 0.1*省份特有因子 噪声。这里 60%、30%、10% 的权重分配是理解全文的主线——它预设了「全国共同因子是主要方差来源」的结论后续的方差分解本质上是在验证这个数据生成过程。你换成真实数据时这个权重关系会变成模型估计的结果但代码逻辑完全复用不需要改结构。注意代码里区域划分有个容易看走眼的地方[0,1,8,9,10,13,14,18,20]是东部省份的索引位置但海南索引 20在论文中的区域归属其实有时会被归到西部[3,4,11,15,16,17]是中部六省[21,22,23,24,25,26,27,28,29,30]是西部十二省含内蒙古。如果你直接跑论文附带的后续分析代码南、海南这些省份的分组会决定因子提取结果建议对照自己研究的行政区划口径做调整。标准化用StandardScaler().fit_transform()对整个 DataFrame 一次完成每一列每个省份被减去均值、除以标准差。这一步对后续DynamicFactor模型的收敛很重要——不平滑的量纲差异会导致状态空间模型在数值优化时陷入局部最优。2.3 数据准备的三个硬性要求用这份代码处理真实数据时请务必满足以下条件数据频率必须一致代码里freqQ生成季度索引真实数据如果是年度数据要么重采样为季度但年度数据插值到季度会引入伪信息要么把date_range参数改成freqA并相应调整周期数。面板必须平衡31 个省份在 116 个季度里都不能有缺失值。真实数据中重庆 1997 年才直辖、西藏部分年份数据缺失这些都会直接导致DynamicFactor报错或因子估计失真。时间跨度要有足够周期数论文覆盖 29 年的季度数据约 5 个完整经济周期。如果只有 10 年数据周期数不足动态因子模型很难把周期成分从趋势和噪声中分离出来。3. 层级动态因子模型实现DynamicFactor 参数、因子提取与贡献度计算3.1 模型构建与参数选择from statsmodels.tsa.statespace.dynamic_factor import DynamicFactor east_provinces [北京, 天津, 河北, 上海, 江苏, 浙江, 福建, 山东, 广东, 海南] central_provinces [山西, 安徽, 江西, 河南, 湖北, 湖南] west_provinces [内蒙古, 广西, 重庆, 四川, 贵州, 云南, 西藏, 陕西, 甘肃, 青海, 宁夏, 新疆] northeast_provinces [辽宁, 吉林, 黑龙江] national_model DynamicFactor(data_scaled, k_factors1, factor_order1) national_results national_model.fit(maxiter1000) east_model DynamicFactor(east_data, k_factors1, factor_order1) east_results east_model.fit(maxiter1000)DynamicFactor是 statsmodels 中专门处理动态因子模型的状态空间实现。k_factors1表示提取一个共同因子factor_order1表示这个因子服从 AR(1) 过程。这两个参数不是随便选的——论文假设存在一个全国共同因子主导所有省份所以k_factors1对应这个理论前提而经济周期本身具有惯性AR(1) 的设定允许因子在相邻季度之间存在自相关比静态因子模型更贴合宏观数据的时间依赖特征。maxiter1000是最大迭代次数。状态空间模型的参数估计用最大似然法通过数值优化算法迭代求解。真实数据中如果迭代到 1000 次还没收敛你会看到Warning: Maximum iterations reached这时需要把maxiter提到 2000 或 5000但也要警惕是否模型设定有问题比如该用k_factors2。模型拟合后因子提取用的是卡尔曼滤波的结果national_factor national_results.factors.filtered[0] east_factor east_results.factors.filtered[0]factors.filtered返回的是滤波后的因子估计值[0]取第一个因子因为k_factors1。filtered表示只用了截至当前时点的信息做估计对应的时间点是每季度末如果要用的信息可以用factors.smoothed。论文图表中展示的因子序列通常是smoothed版本因为平滑估计利用了全部样本信息真实波动模式还原度更高。3.2 因子贡献度分解的核心算法def calculate_variance_decomposition(results, province_data): loadings results.params[results.param_names[results.param_names.str.startswith(loading)]] total_variance np.var(province_data, axis0) factor_variance (loadings * np.var(results.factors.filtered[0]))**2 explained_ratio factor_variance / total_variance return explained_ratio这个函数是整个复现过程中最容易理解错的部分。results.params包含所有估计参数param_names里以loading开头的就是因子载荷——也就是每个省份对共同因子的敏感度系数。loadings * np.var(results.factors.filtered[0])这一步计算的是因子载荷乘以因子方差然后整体平方得因子贡献的方差。这里其实有个统计学上值得商榷的地方严格来说因子解释的方差应该是loadings² * var(factor)而不是(loadings * var(factor))²。代码里的写法是(loadings * np.var(...))**2数学上相当于把因子方差的平方又乘了一次数值会被放大。实际跑出来的解释度比例不一定会超过 1因为后面还要除以总方差但这个计算方式确实和标准公式有偏差。我建议改成factor_variance (loadings**2) * np.var(results.factors.filtered[0])这样才是「因子载荷平方乘以因子方差」的正确表达。如果按原代码跑全国因子解释度普遍偏高会让你误以为共同因子主导一切。论文结论本身是这个方向但拿错误公式去验证结论就失去复现的意义了。3.3 因子贡献度结果的解读方式results_df[全国因子解释度] national_explained results_df[区域因子解释度] east_explained # 按区域填充 results_df[剩余波动] 1 - results_df[全国因子解释度] - results_df[区域因子解释度]剩余波动的计算逻辑是1 - 全国因子解释度 - 区域因子解释度它代表省份特有因子加上随机扰动解释的部分。这个值越高说明该省份越「与众不同」。在论文的叙事框架下东北三省的剩余波动显著高于其他区域——辽宁、吉林、黑龙江的省级特有成分占主导意味着这些省份的经济周期和全国及区域层面都不同步。你在复现时可重点关注这个指标的区域分布它是判断「异质分化」是否成立的直接依据。输出结果建议用堆叠柱状图展示代码里results_df[[全国因子解释度, 区域因子解释度, 剩余波动]].plot(kindbar, stackedTrue)已经写好了。但注意原代码在绘图前没有剔除标准化导致的负值——StandardScaler会让某些省份的标准化序列有负值方差分解结果理论上应该是非负的但如果你改用了修正后的公式而数据本身有异常柱状图可能出现负高度视觉上会很难看。遇到这种情况直接把负值截断为 0并在论文里注明是调整后的结果。4. 区域协同与异质分化检验格兰杰因果、滚动相关与 K-means 聚类4.1 格兰杰因果检验的滞后阶数选择from statsmodels.tsa.stattools import grangercausalitytests def analyze_granger_causality(factors, names, maxlag4): for i in range(len(factors)): for j in range(len(factors)): if i ! j: test_data pd.DataFrame({x: factors[i], y: factors[j]}) result grangercausalitytests(test_data, maxlagmaxlag, verboseFalse) p_values [round(result[lag1][0][ssr_ftest][1], 4) for lag in range(maxlag)] if any(p 0.05 for p in p_values): print(f{names[i]}是{names[j]}的格兰杰原因)格兰杰因果检验在这里的作用是验证「东部引领、中西跟随」的领先-滞后关系。maxlag4意味着最多检验 4 个季度的滞后效应——恰好对应一年。经济周期传导通常需要 2 到 8 个季度4 作为默认值是个折中方案。注意代码里对 p 值的判断用了any(p 0.05)也就是只要 4 个滞后阶数里任意一个显著就判定存在格兰杰因果关系。这有个多重比较的问题——4 次检验里出现一次假阳性的概率约为 18.5%1 - 0.95⁴。严谨做法是对 p 值做 Bonferroni 校正p 0.05/4 0.0125或者只看最小 p 值对应的滞后阶是否显著。原代码的宽松标准在复现论文结论时问题不大但如果你要拿这个结果发表建议用校正后的阈值再跑一遍。另一点容易被忽略grangercausalitytests要求输入的序列是平稳的。模型提取的因子是 AR(1) 过程只要自回归系数小于 1 就是平稳的一般没问题。但如果你换用自己的因子提取方法比如用 PCA 提取静态因子得到的因子可能有单位根跑格兰杰检验前必须做 ADF 检验。4.2 滚动相关性观察同步性的时间演变window_size 20 # 5年窗口季度数据 rolling_corr region_factors_df[东部].rolling(window_size).corr(region_factors_df[中部])滚动相关性的窗口设为 20 个季度正好 5 年。这个窗口长度兼顾了平滑性和时效性——太短如 8 个季度相关性波动太剧烈太长如 40 个季度又会把 2008 年金融危机前后的结构变化抹平。代码注释里的「5年窗口」对应的就是这个参数。画出来的图是一条随时间变化的曲线你需要重点关注两个时间节点2008 年金融危机后东部-中部的滚动相关性是否显著上穿 0.6以及 2012 年后是否出现下滑。论文中的「东部-中部同步性在 2008 年后显著增强」这个结论就是靠这条曲线撑起来的。如果你用真实数据复现发现滚动相关性没有在某个事件节点出现明显跳变那要么是数据口径不同要么是区域分组的假设需要调整。4.3 K-means 聚类与俱乐部收敛检验from sklearn.cluster import KMeans idio_components results_df[剩余波动].values.reshape(-1, 1) kmeans KMeans(n_clusters3, random_state42).fit(idio_components) results_df[俱乐部] kmeans.labels_n_clusters3来源于论文「俱乐部收敛」的理论预设——预期东部、中部、西部各自形成一个收敛俱乐部东北是异类。random_state42固定了随机种子保证每次运行聚类结果一致这是可复现研究的基本要求。这里有个方法论层面的问题值得你注意K-means 聚类用「剩余波动」作为唯一特征只把每个省份浓缩成一个数字平均剩余波动比例丢掉了时序信息。严格来说这不是标准的俱乐部收敛检验——Phillips-Sul 收敛检验需要基于横截面方差的时序演变来判断而这里只是「对离散度聚类再比较组内方差」。在论文复现的语境下这样做够用能直观展示「哪些省份的波动模式更接近」但如果你的目标是发表实证论文建议用 Phillips-Sul 方法做补充验证。输出结果会有一个交叉表results_df.groupby([区域, 俱乐部]).size()。预期的理想结果是东部省份集中在俱乐部 1、中部在俱乐部 2、西部在俱乐部 3而东北三省的标签比较分散。如果聚类结果出现「东部省份被拆到两个俱乐部」的情况说明东部内部也有不小分化这本身也是一个可报告的发现。5. 动态因子模型避坑五个高频翻车点及排查方案5.1 模型报错Cholesky 分解收敛失败现象运行national_model.fit(maxiter1000)时出现ValueError: The state covariance matrix is not positive definite或ConvergenceWarning模型无法正常收敛。原因这通常是两类问题导致的。第一输入数据存在线性相关比如两个省份的 GDP 增长率完全一致导致方差协方差矩阵奇异第二factor_order1设定下因子方差的初始值不合适数值优化时状态协方差矩阵在迭代过程中失去了正定性。解决先检查数据相关性用data_scaled.corr()找出相关系数超过 0.99 的省份对考虑删掉一个或做去均值处理。如果数据没问题改参数初始值——DynamicFactor接受initialization参数可以尝试initializationstationary强制使用平稳初始化或者手动设置更保守的初始方差。5.2 因子解释度超过 100%现象计算出来的全国因子解释度 区域因子解释度总和大于 1剩余波动出现负值。原因如我在 3.2 节提到的原代码的方差分解公式(loadings * np.var(factor))**2在数学上不正确它会高估因子贡献。另外如果不同层级模型全国与区域的数据标准化方式不一致也会引入额外差异。解决改用(loadings**2) * np.var(factor)的标准公式。如果修正后仍有超过 100% 的省份检查该省份是否在区域分组中被重复纳入比如海南同时在东部和西部列表里。5.3 格兰杰检验结果不稳定现象同一组因子数据跑两次格兰杰因果检验结论有时显著、有时不显著。原因grangercausalitytests的 F 检验依赖残差平方和的精确计算如果你的因子序列长度较短少于 30 个观测检验统计量的分布不稳定p 值在小样本下会出现较大波动。解决先确认样本量116 个季度数据在maxlag4时每个回归的有效样本约 112问题不大。但如果换成了年度数据29 个观测maxlag4就完全不合适了——这时改maxlag2或 1同时用 AIC/BIC 选择滞后阶数而不是固定一个值。5.4 K-means 聚类结果每次跑都不一样现象不设置random_state直接跑 K-means聚类结果在不同运行中会变化有时东三省被分到不同俱乐部。原因K-means 的初始质心是随机选择的不同的初始点可能收敛到不同局部最优。如果数据分布不够分离比如剩余波动在各省份间差距不大聚类结果就很不稳定。解决设置random_state42或其他固定值保证可复现。更稳妥的做法是跑多次聚类比如n_init100sklearn 默认是 10选惯性最小的结果。另外如果用 3 个俱乐部聚类效果不好轮廓系数小于 0.3考虑改用层次聚类看树状图判断俱乐部数量是否为 3。5.5 滚动相关性曲线锯齿状剧烈跳动现象画出来的rolling_corr曲线噪声很大无法分辨结构变化2008 年前后看不出明显差异。原因rolling(window20).corr()对每个窗口内的相关性做等权平均如果窗口内出现极端值某个季度某省数据异常波动相关性会被拉偏。另外如果两个区域的因子同时包含趋势成分滚动窗口内的子样本相关性会受趋势影响出现「伪相关」波动。解决先对因子做 HP 滤波或差分处理提取周期成分再进行滚动相关计算。窗口大小也可以调整——按论文口径用 20 个季度5 年合理但如果曲线仍然太噪可以试window287 年做平滑。6. 结果验证与拓展三力均衡指数和对照实验的设计技巧论文摘要中提到了「三力均衡指数」的概念原代码里没有给出具体实现。我在复现时自己补了一个版本这可以作为拓展实验的方向定义均衡指数为全国因子解释度、区域因子解释度、剩余波动三者之间的均衡程度用三维向量的归一化熵来度量。具体做法是把每个省份的三个解释度归一化为概率分布计算信息熵熵越高说明三方解释力越均衡熵越低说明一个层级主导。这个指数可以直接嵌入现有代码from scipy.stats import entropy def calculate_balance_index(row): components np.array([row[全国因子解释度], row[区域因子解释度], row[剩余波动]]) components np.maximum(components, 0) if components.sum() 0: return 0 probs components / components.sum() return entropy(probs) results_df[均衡指数] results_df.apply(calculate_balance_index, axis1)这个指数的优势在于把三个维度的解释力压缩到一个标量可以直接做省份间排序和区域间比较。期望结果是东北三省的均衡指数偏低因为剩余波动主导而东部和中部省份的均衡指数偏高说明全国因子和区域因子的解释力相对平衡。这和你前面做方差分解的结果互为印证比单纯堆叠柱状图更有政策解释力。在验证模型有效性时我建议你做两组对照实验。第一组是把k_factors从 1 改成 2重新拟合全国模型然后对比两个因子分别解释的省份集合。如果第二个因子主要解释东北省份那说明论文「东北独立」的结论其实可以通过双因子结构来捕捉单因子模型未必是最佳设定。第二组是随机打乱省份的区域分组后重跑模型如果随机分组下的区域因子解释度显著下降说明区域层级的划分是有实际信息量的不是数据噪声驱动的结果。验证的最后一步是看残差的自相关图。DynamicFactor拟合后用national_results.resid做 Ljung-Box 检验p 值应大于 0.05说明残差没有显著自相关。如果检验失败说明模型遗漏了某些动态结构考虑增加factor_order或引入观测方程的自回归项。这套代码我用完整流程跑过一遍从数据生成到结论输出耗时约三分钟主要时间花在DynamicFactor的最大似然迭代上。从那以后我每次做区域经济周期分析都强制走一遍「方差分解 → 格兰杰检验 → 滚动相关 → 聚类验证」的完整链路先把合成数据的已知结构拟合出来再换真实数据找差异。这套思路帮我省了至少两周的排错时间希望帮到你。本文还有配套的精品资源点击获取
