ARTICLE DETAIL

建站实战干货

来自一线的建站与推广经验沉淀,每一条都经过真实交付验证。

用Python实现动态CGE模型:从静态均衡到递归动态模拟

2026/10/3 11:16:03 拓冰建站 浏览量
用Python实现动态CGE模型:从静态均衡到递归动态模拟 简介面向经济学研究者与政策分析师的Python动态CGE模型完整实现方案覆盖数据清洗、柯布-道格拉斯生产函数设定、市场均衡求解、结果可视化及tkinter图形界面设计可用于宏观经济政策、贸易政策与环境经济分析等场景。资源包为1个docx文档压缩后仅31KB正文包含分模块代码讲解、求解思路、可能遇到的问题及未来改善路径便于按段落逐步复现。目前已有238人学习浏览。文档从数据导入到GUI运行给出体系化代码并提醒数据质量与计算效率的权衡适合希望快速上手动态一般均衡建模的程序员、研究生及研究人员。整体结构清晰兼具理论说明与工程实现细节可帮助读者理解政策冲击对经济系统的动态影响并在此基础上扩展模型、支持更多数据集或集成高级计算技术从而提升模型准确性与效率。1. 动态CGE模型为什么用Python重写一遍是值得的动态CGE模型可计算一般均衡模型在贸易、投资、税收政策评估中几乎是标准工具传统上被GAMS和GTAP生态垄断。但如果你不想被黑匣子绑架或者只是做课程设计、政策模拟的原型Python反而是最合适的选择Numpy做数值运算、Scipy的求解器做均衡搜索、pandas做结果汇总一套几十行的代码就能把静态CGE的引擎跑起来再套一个时间循环就是动态CGE。这篇文章用一个小型两部门递归动态CGE做例子从方程写到可运行的代码让读者看完能动手复现也能按自己的数据换掉参数。适合已经有Python基础和微观经济学概念、但被GAMS语法劝退的人。我会把每一段代码拆开讲清楚参数含义和坑在哪里你照着敲一遍就能理解动态CGE的整个落点。2. 从静态CGE到动态CGE模型设计与核心方程2.1 为什么选择“递归动态”而不是“跨期优化”动态CGE模型分成两类一类是无限期跨期优化家庭在完美预期下选择消费和储蓄路径需要求解欧拉方程和终值条件另一类是递归动态把每一期当作独立静态均衡来解期与期之间只通过外生更新的劳动供给和资本积累方程连接。递归动态虽然不像跨期优化那样有坚实的微观福利基础但它可解释性好、收敛容易、对初学者友好而且绝大多数政策模拟项目用的是它。本文选择递归动态因为Python写起来最直接静态均衡是一个方程组动态只是把这个方程组放在循环里反复解而不是去解一个大规模的动态规划或最优控制问题。递归动态里储蓄率通常是外生常数这是一个刻意的简化。如果你需要把储蓄率内生化可以把模型升级成拉姆齐式但那需要引入Bellman方程或打靶法代码复杂度会上升一个量级。对第一次用Python实现动态CGE的人先跑通递归动态再往复杂方向扩是最稳妥的技术路径。2.2 两部门CGE的方程结构与均衡条件这里的模型包含两个生产部门消费品部门C和投资品部门I。每个部门使用劳动和资本两种要素技术用Cobb-Douglas生产函数描述Y_s A_s * K_s^α_s * L_s^(1-α_s)其中s∈{C,I}α_s是部门s的资本产出弹性。假设规模报酬不变因此企业零利润条件成立商品价格等于单位成本。单位成本函数的推导结果是c_s(w,r) (1/A_s) * (w/(1-α_s))^(1-α_s) * (r/α_s)^α_s这里w是工资率r是资本租金率。零利润条件写成P_C c_C(w,r) P_I c_I(w,r)家庭拥有全部劳动和资本存量其收入是Y_H w * L_t r * K_t家庭把固定比例s储蓄剩余用于消费。于是消费品的需求量和投资品的需求量分别是C_d ((1-s) * Y_H) / P_C I_d (s * Y_H) / P_I商品市场出清条件是Y_C C_d Y_I I_d要素市场出清条件是l_C(w,r)*Y_C l_I(w,r)*Y_I L_t k_C(w,r)*Y_C k_I(w,r)*Y_I K_t其中l_s和k_s是单位产出的劳动和资本需求由成本最小化推得l_s (1/A_s) * ((1-α_s)r / (α_sw))^α_s k_s (1/A_s) * (α_s*w / ((1-α_s)*r))^(1-α_s)这里有一个瓦尔拉斯红利五个方程中有一个是冗余的所以固定P_C1作为价格基准numeraire剩下的未知数是w、r、P_I、Y_C、Y_I。我一般保留劳动市场出清方程资本出清方程留作事后校验这样求解器不会因为方程数超过变量数而抱怨。实际求解时资本市场出清的误差通常在1e-8以下如果误差偏大就说明参数校准有毛病。2.3 参数校准让基准年经济被模型“复制”出来动态CGE的起点是基年社会核算矩阵SAM。如果拿不到完整SAM也可以用简化校准用基年的生产数据和要素收入份额反推α_s。常见做法是假设基年利润和工资构成增值那么α_s的估计值就是资本报酬占部门产出的比例。用最小二乘或直接代入都可以。比如基年部门s的产出Y_s、劳动投入L_s、资本投入K_s已知工资w和租金r从SAM看那么α_s的校准公式是α_s (r * K_s) / (P_s * Y_s)在代码里我会先用一套假定的基年数据做示范把上面公式代入保证基准年模型的产出、就业、要素收入完全等于输入数据。这一步不通过后面任何动态情景模拟都是地基不牢。3. 用Python实现静态均衡求解器最小可运行代码3.1 数据结构与参数定义先把参数集中在字典里方便后续校准和情景修改。这样做的好处是动态循环里换参数非常容易不需要改动求解函数本体。import numpy as np from scipy.optimize import root, fsolve # 基础参数 params { alpha: {C: 0.35, I: 0.25}, # 各部门资本产出弹性 A: {C: 1.0, I: 1.0}, # 全要素生产率系数 s_rate: 0.25, # 家庭储蓄率 delta: 0.06, # 资本折旧率 g_L: 0.02, # 劳动供给年增长率 g_A: 0.015, # 全要素生产率年增长率 L0: 10.0, # 基年劳动供给 K0: 25.0 # 基年资本存量 }上面的α值不是拍脑袋它们应该来自基年SAM校准。这里直接拿去运行会得到一个“假想经济”但代码结构是真能跑的。如果需要自己的数据就把α替换成校准值。3.2 单位成本函数与要素需求函数我把单位成本和单位要素需求封装成独立函数让均衡方程清晰可读也方便后面做成本分解或替代弹性扩展。def unit_cost(w, r, sector, params): a params[alpha][sector] A_s params[A][sector] return (1.0 / A_s) * (w / (1.0 - a))**(1.0 - a) * (r / a)**a def unit_factor_demand(w, r, sector, params): a params[alpha][sector] A_s params[A][sector] l (1.0 / A_s) * ((1.0 - a) * r / (a * w))**a k (1.0 / A_s) * (a * w / ((1.0 - a) * r))**(1.0 - a) return l, k注意这里的指数不能写错。Cobb-Douglas成本函数对CD生产函数而言必然满足这一形式。要是α或A修改后求解失败第一步先检查这个函数的计算结果是否为正w和r必须是严格正数否则指数会生成nan。3.3 均衡方程组与求解器的构造静态均衡的核心是一个五个未知数、五个方程的系统。我使用scipy.optimize.rootmethod选择lmLevenberg-Marquardt因为它对初始猜测的敏感度比默认的hybr低一些尤其适合这种非线性价格方程。P_C 1.0 # 价格基准消费品价格定为1 def static_equilibrium(vars, L_supply, K_supply, params): w, r, P_I, Y_C, Y_I vars # 零利润条件 eq1 P_C - unit_cost(w, r, C, params) eq2 P_I - unit_cost(w, r, I, params) # 家庭收入与储蓄 Y_H w * L_supply r * K_supply C_d (1.0 - params[s_rate]) * Y_H / P_C I_d params[s_rate] * Y_H / P_I # 商品市场出清 eq3 Y_C - C_d eq4 Y_I - I_d # 劳动市场出清资本出清留作校验 l_C, k_C unit_factor_demand(w, r, C, params) l_I, k_I unit_factor_demand(w, r, I, params) eq5 l_C * Y_C l_I * Y_I - L_supply return [eq1, eq2, eq3, eq4, eq5] def solve_equilibrium(L_supply, K_supply, params): # 初始猜测w1, r0.1, P_I1, 产出按劳动供给占大头估计 x0 np.array([1.0, 0.1, 1.0, L_supply*0.8, L_supply*0.2]) sol root(static_equilibrium, x0, args(L_supply, K_supply, params), methodlm) if not sol.success: raise RuntimeError(均衡求解失败: sol.message) w, r, P_I, Y_C, Y_I sol.x # 事后校验资本出清 l_C, k_C unit_factor_demand(w, r, C, params) l_I, k_I unit_factor_demand(w, r, I, params) K_demand k_C * Y_C k_I * Y_I capital_check K_demand - K_supply return { w: w, r: r, P_I: P_I, Y_C: Y_C, Y_I: Y_I, K_demand: K_demand, capital_check: capital_check }初始猜测是一个容易翻车的地方。我给的x0里r0.1是凭经验资本租金率远高于折旧率但又不是激进到让单位成本变成负数。如果你改了α或A最好先跑一次基准年求解把输出的下一期均衡解作为新的初始猜测。更稳妥的做法是动态循环里把上一期的解作为当前期的x0后面可以看到。3.4 基准年校准与验证脚本在动态模拟前先确认静态求解器能复现基准年。假如基年L10、K25那么求解出的要素收入和产出应该落在合理范围。def calibrate_baseline(params): L_base params[L0] K_base params[K0] eq solve_equilibrium(L_base, K_base, params) print(基年静态均衡结果) print(工资率 w , eq[w]) print(资本租金率 r , eq[r]) print(投资品价格 P_I , eq[P_I]) print(消费品产量 Y_C , eq[Y_C]) print(投资品产量 Y_I , eq[Y_I]) print(资本出清偏差 , eq[capital_check]) calibrate_baseline(params)这段脚本的意义是如果capital_check超过1e-6说明零利润方程、要素需求函数或家庭收入公式里有bug。我见过有人在单位成本函数里把α和1-α写反结果资本出清偏差巨大而求解器仍然能返回一组数因为‘lm’求的是最小二乘解。所以校准校验必须独立保留。4. 动态递推与情景模拟把静态求解器装进时间循环4.1 资本积累与要素更新方程递归动态的“动态”体现在期与期之间资本存量的更新上。每期静态均衡解出投资品产量Y_I就是当期的实际投资I_t。期末的资本存量按下式更新K_{t1} (1 - δ) * K_t I_t劳动供给按外生人口增长率增长L_{t1} (1 g_L) * L_t全要素生产率A也可以逐年增长A_{t1,s} A_{t,s} * (1 g_A)这三条规则决定了整个动态路径。投资品价格P_I会影响名义投资额但实物资本积累用的是数量Y_I所以更新方程里的I_t是投资品数量而不是投资额。这里很容易混淆特别是你从GAMS代码转过来时GAMS里常用价格乘数量而这里我把价格和数量拆开。4.2 动态模拟主循环逐年递推与结果保存把静态求解器放进一个循环每期都用最新的L_t和K_t求解然后更新。为了稳定把上一期的解作为下一期的初始猜测。def simulate_dynamic(years, params): L params[L0] K params[K0] A_base {s: params[A][s] for s in params[A]} results [] # 上一期均衡解用于初始猜测 prev_x None for t in range(years): # 当期生产率基年A乘以增长率 for s in params[A]: params[A][s] A_base[s] * (1.0 params[g_A])**t eq solve_equilibrium(L, K, params) I_t eq[Y_I] # 收集结果 results.append({ year: t, L: L, K: K, Y_C: eq[Y_C], Y_I: eq[Y_I], w: eq[w], r: eq[r], P_I: eq[P_I], K_demand: eq[K_demand] }) # 更新资本和劳动 K (1.0 - params[delta]) * K I_t L (1.0 params[g_L]) * L # 准备下一个周期的初始猜测 prev_x [eq[w], eq[r], eq[P_I], eq[Y_C], eq[Y_I]] return results results simulate_dynamic(20, params) for rec in results[:5]: print(rec)这个循环有几点需要注意。第一params[A]在循环内被直接修改所以A_base必须提前拷贝否则第二次循环时生产率会重复累加最终A变成(1g_A)^(t*(t1)/2)而不是(1g_A)^t。第二prev_x这里只是摆在那里真正要用它当初始猜测需要改solve_equilibrium让它接受x0参数。下面给出升级版def solve_equilibrium_with_guess(L_supply, K_supply, params, x0None): if x0 is None: x0 np.array([1.0, 0.1, 1.0, L_supply*0.8, L_supply*0.2]) sol root(static_equilibrium, x0, args(L_supply, K_supply, params), methodlm) # ... 后半段同上一版然后在simulate_dynamic里每期求解前把prev_x传给x0。实际经验是用上一期解当初始猜测绝大多数年份一次收敛不用的话遇到较大的技术进步冲击就可能出现“求解成功但数值异常”的假收敛。4.3 情景模拟储蓄率冲击与技术冲击动态CGE最常见的用法是做政策或环境参数冲击分析。比如想模拟储蓄率从0.25永久提高到0.35对经济的影响def simulate_scenario(years, params, scenario): p params.copy() p[s_rate] scenario.get(s_rate, params[s_rate]) p[g_A] scenario.get(g_A, params[g_A]) return simulate_dynamic(years, p) # 基准情景 base simulate_dynamic(30, params) # 高储蓄情景 params_high_saving params.copy() params_high_saving[s_rate] 0.35 high_saving simulate_dynamic(30, params_high_saving) # 对比期末资本存量 print(基准期末K , base[-1][K]) print(高储蓄期末K , high_saving[-1][K])这里需要留意params.copy()是浅拷贝嵌套字典沿用旧引用。如果你在scenario里改s_rate不会影响params[s_rate]但你若在循环里修改params[A]由于浅拷贝共享那个子字典会串数据。安全做法是用copy.deepcopy。这是一个容易踩的接口设计坑后面会专门说。4.4 结果输出与增长率计算动态模拟得到的是逐年结果需要计算增长率来判断模型是否走在平衡增长路径上。我通常用numpy的diff求对数增长率ys np.array([r[Y_C] r[Y_I] for r in results]) gdp_growth np.diff(np.log(ys)) print(GDP对数增长率前10期, gdp_growth[:10])这里把消费品产出和投资品产出直接相加。严格的实际GDP应该用基准年价格加权即固定P_C1和P_I,base计算链式加权数量。但小模型里如果价格P_I变化不大简单加总用作路径增长率是够的。正式报告时建议用Laspeyres数量指数。5. 动态CGE建模的避坑与常见问题排查5.1 求解器不收敛初始猜测与参数范围现象scipy.optimize.root返回“The iteration is not making good progress”或者抛出RuntimeError。原因最常见的是初始猜测离真解太远尤其是r和P_I的量纲不对。第二个常见原因是参数α或A更新后单位成本函数计算出现负数或nan比如w或r为负导致指数运算错误。解决我在前面已经给了两个措施——用上一期解作为当前期初始猜测给r和w加正数约束。还可以在static_equilibrium函数开头做一次变量检查如果vars里有非正值直接返回一个很大的残差列表让求解器离开非法区域def static_equilibrium(vars, L_supply, K_supply, params): w, r, P_I, Y_C, Y_I vars if min(vars) 0: return [1e6, 1e6, 1e6, 1e6, 1e6] ...这招在Levenberg-Marquardt下特别有效因为它本质是数值梯度搜索用大残差可以挡住负值方向。注意这不能保证最优解正性但能避免搜索过程崩掉。5.2 价格归一化与名义变量漂移现象动态模拟里P_I逐年上升或下降但实物变量增长正常。原因价格基准只固定了P_C1没有固定货币总量。在递归动态里如果劳动生产率和资本存量都在增长名义收入也会增长但P_C固定意味着总价格水平会变。这不一定是错误但如果你把名义工资w当作实际工资来看就会误判。解决明确区分实际变量和名义变量。如果想看实际工资就用w/P_CP_C1所以这里w就是实际工资如果想看实际资本租金用r/P_C。投资品价格P_I的变化会影响投资名义量但不影响实物投资数量Y_I。我一般会额外存储w/P_C和r/P_C作为输出变量避免分析时写错。另外不要在代码里同时固定P_I和P_C那样方程系统会被瓦尔拉斯定律打乱求解器反而更容易失败。5.3 负投资或负资本存量出现现象某年的Y_I变成负值或者K在后期变成负数。原因递归动态模型里投资I_t s*Y_H / P_I。如果参数设置导致储蓄率过高或折旧率过高而投资品价格又太低那么名义储蓄不足以维持正的净投资。更隐蔽的原因是要素市场出清条件被删掉后资本需求远大于供给导致资本租金r飙高家庭收入变大但投资品生产消耗太多劳动挤占了消费品生产使得Y_I为负。解决先检查参数量纲。储蓄率一般不超过0.4折旧率在0.04到0.15之间是合理区间。如果Y_I为负把s_rate调低或delta调高再观察资本出清偏差。如果想约束投资数量非负可以在均衡方程中加不等式但那需要换成scipy.optimize.minimize和约束优化本文不展开。5.4 基年校准与SAM数据不一致现象基准年求解出的要素收入份额和SAM里对不上比如劳动收入占总产出比重不是1-α。原因α是校准参数不是外生给定的。很多初学者直接把文献里的α抄进来却忘了α必须和基年SAM的资本-劳动比以及要素价格一致。解决用基年数据反推α。假设基年工资总额W_total、资本总额R_total部门s的资本产出弹性α_s R_total_s / (P_s * Y_s)在代码里校准自动完成def calibrate_alpha(Y_base, K_base, L_base, w_base, r_base, sector): return r_base * K_base / (w_base * L_base r_base * K_base)如果手头有完整SAM这一步更简单α_s 部门s资本报酬 / 部门s总产出。校准完必须用基准年求解器验证资本出清偏差小于1e-6才继续。5.5 技术进步参数g_A被重复计入现象模拟第10年A值比理论值大得多增长路径明显加速异常。原因动态循环内直接修改了params[A]而params[A]又被下一次循环继续叠加形成了连乘中的连乘。这是所有递归动态模拟最容易犯的错。解决在动态模拟函数开头保存A_base循环内临时计算当期A不要直接写入paramsA_curr {s: A_base[s] * (1.0 params[g_A])**t for s in A_base}然后在调用静态求解器前把params[A]替换成A_curr但循环结束后要还原。如果用copy.deepcopy生成params副本也能避免污染外部状态。我个人的习惯是凡是要在模拟中修改的参数一律先deepcopy绝不原地改动外部传入的字典。6. 验证模型动态特性稳态检验与收敛性诊断写完动态循环不要急着用来出政策结论。第一件事是验证模型能不能走到平衡增长路径。递归动态CGE的理论稳态是资本存量、产出、消费都以同一个增长率g*增长这个增长率由劳动增长率和外生技术增长率决定。对于CD生产函数人均产出增长率等于g_A/(1-α)总量增长率再加g_L。你可以用这个公式检验模拟结果。我在模拟结束后会这样验证import numpy as np def check_balanced_growth(results, params): # 取后十年增长率 ys np.array([r[Y_C] r[Y_I] for r in results]) growth_end np.mean(np.diff(np.log(ys[-10:]))) # 理论增长率a是针对总经济的聚合资本份额取两个部门加权 a_agg (params[alpha][C] params[alpha][I]) / 2 g_theory params[g_A] / (1.0 - a_agg) params[g_L] return growth_end, g_theory g_sim, g_theory check_balanced_growth(results, params) print(f模拟末期增长率: {g_sim:.4f}) print(f理论平衡增长率: {g_theory:.4f})如果两者偏差超过千分之一说明模型还没有收敛到稳态或者部门加权α的估计太粗糙。此时不要怀疑理论而是检查动态循环是否真的让要素市场出清以及投资品价格P_I是否出现了异常趋势。除了增长率我还会检查资本-产出比K/Y是否最终水平稳定。在递归动态模型里K/Y会逐步收敛到一个常数k_y np.array([r[K] / (r[Y_C] r[Y_I]) for r in results]) print(K/Y 序列后5期:, k_y[-5:])如果K/Y还在明显上升或下降通常意味着折旧率和储蓄率的组合与稳态不匹配。我自己的经验是先跑100期看最后20期的K/Y波动如果波动幅度小于0.01%才认为模型合格。这也是为什么我动态模拟总喜欢用deepcopy——跑情景对比时不会污染基准参数。最后再提醒一个小习惯永远保留基年校准脚本。每次改参数后直接跑一遍基年求解如果capital_check明显变大就不要继续做情景分析回过头查单位成本函数和要素份额公式。这比在100期模拟结果里找问题省时间得多。希望这些代码和排查思路能帮你在Python里亲手把动态CGE跑通并少走一点弯路。本文还有配套的精品资源点击获取