做气候政策分析这些年我越来越觉得R语言实现DICE模型这件事被严重低估了。DICEDynamic Integrated Climate-Economy model动态综合气候-经济模型是诺德豪斯William Nordhaus的核心成果也是他拿到2018年诺贝尔经济学奖的代表作之一。这个模型把全球看作一个“大账本”把经济产出、碳排放、碳循环、气温变化和社会福利串在一条因果链上用来回答一个现实问题如果选择不同的减排路径人类到底要付出多少经济成本又能避免多少气候损失。之前我在做双碳路径测算时用R从零搭过一版DICE模型整个过程踩了不少坑也把模型里的每个方程彻底“啃”了一遍。这篇文章不打算复读教科书而是想把这些实操经验做个完整梳理包括模型怎么从公式变成R代码、减排路径怎么量化设定、模拟结果怎么解读以及几个真正容易栽跟头的地方。如果你手里有R语言基础想把手里的宏观气候模型跑起来这篇文章应该能帮你省下不少摸索时间。1. DICE模型到底在算什么先建立直观认知在最开始我建议你先别急着看代码而是把这个模型的“骨架”想清楚。DICE模型从一个非常朴素的角度出发全球经济活动会排放二氧化碳二氧化碳在大气中累积造成温室效应温室效应导致温度上升温度上升对经济产生损害而减排则可以降低损害但也有成本。这个链条就是DICE模型的核心逻辑模型的目标是找到一条“总成本最小”的路径——减排成本和气候损害成本加起来最小或者从另一个角度说让折现后的社会福利最大化。1.1 模型的五大模块与信息流向DICE模型在结构上可以拆成五大模块它们之间是顺序驱动的关系经济模块用生产函数描述全球GDP产出资本积累受储蓄率影响产出的一部分用于消费一部分用于减排投入。排放模块经济活动产生二氧化碳排放排放量由产出规模和单位产出的碳排放强度决定减排比例会直接压缩排放量。碳循环模块大气、浅层海洋、深层海洋三个碳库之间互相转移模拟排放到大气中的碳最终去了哪里。气候模块大气碳浓度升高产生辐射强迫进而驱动地表温度和海洋温度变化。福利模块将人均消费转化为效用并按折现率加总得到社会福利函数这也是模型优化的目标函数。这个结构的巧妙之处在于它不是单方向推演而是有反馈的。温度变化反过来会损害经济产出损害函数从而影响未来的消费和资本积累同时现在的减排投入会减少当期可用于消费和投资的资源。模型就是在“当期投入”和“未来收益”之间做动态权衡。1.2 这个模型能回答什么问题我在实际应用中DICE模型主要被用来回答以下几类问题第一个问题是“不减排会怎样”也就是基准情景Business As UsualBAU下2100年全球升温多少、经济损失多大。第二个问题是“给定碳税水平效果如何”比如每吨碳征收50美元或100美元碳排放轨迹会发生什么变化。第三个问题是“温控目标如何反推减排曲线”比如要实现2摄氏度或1.5摄氏度的温控目标全球减排率每年需要达到多少。第四个问题是“最优碳价是多少”在模型框架下社会福利最大化对应的影子碳价就是“最优碳价”这是政策制定者关注的核心指标。用一句话概括DICE模型可以把“减排”这个抽象概念量化成一条路径、一个碳价、一组温度与福利结果。对于双碳规划、气候政策评估、绿色金融压力测试这些场景它提供的不是单一数字而是一整套动态分析框架。1.3 为什么选择R来实现DICE模型本身不是特定于某一种语言的常见实现包括GAMS、Python和R。我之所以在项目中使用R有三个实际考虑一是R在数据处理和可视化生态上很成熟尤其是ggplot2做政策情景对比图很顺手二是R的优化工具足够用模型本质是动态最优化问题可以用optim、nloptr这类包处理三是后续做参数敏感性分析、蒙特卡洛模拟时R的向量化计算和并行化支持都很方便。需要说明的是R跑DICE模型的速度相比C或Julia会慢一些但模型的周期通常是60期左右计算量完全可控我个人实测一次完整模拟加优化在普通笔记本上只需要几秒到几十秒性能不是瓶颈。2. R语言实现DICE模型核心方程与代码拆解当年我从零开始写DICE模型的R代码时最大的感受是模型的理论框架看起来简单但每个模块的方程都有细节要处理。这里我把核心方程和对应代码一并拆开讲你可以照着搭起一个可运行的最小版本。2.1 时间与状态变量的设定经典DICE-2016模型以10年为一期时间范围设定为60期也就是从2015年模拟到2615年。在实际分析中我们最关心的通常是2100年之前的结果但模型需要更长的地平线来避免“末日效应”——如果只看100年模型会把最后一个周期的资本全部消费掉产生不合理的路径。t_max - 60 # 60期每期10年 t0_year - 2015 # 起始年份 step - 10 # 每期10年 years - t0_year (0:(t_max - 1)) * step状态变量包括资本存量K、大气碳浓度MAT、浅层海洋碳浓度MUP、深层海洋碳浓度MLO、地表温度TAT和海洋温度TLO。控制变量包括储蓄率s和减排率mu模型的优化目标就是找出每个时期这两个变量的最优序列。2.2 经济模块的实现经济模块的核心是生产函数。DICE模型使用的是柯布-道格拉斯形式Q(t) A(t) * K(t)^γ * L(t)^(1-γ)其中A是全要素生产率K是资本存量L是劳动力γ是资本产出弹性。在实际代码中模型通常先计算“毛产出”再扣除减排成本得到净产出# 计算毛产出 calc_gross_output - function(A, K, L, gamma) { A * K^gamma * L^(1 - gamma) } # 减排成本占GDP比例mu为减排率theta为成本指数 calc_abatement_cost_frac - function(pi_param, mu, theta) { pi_param * mu^theta } # 计算净产出 calc_net_output - function(gross_output, abatement_cost_frac) { gross_output * (1 - abatement_cost_frac) }资本积累方程是标准的永续盘存法K(t1) (1-δ) * K(t) s(t) * Q(t)其中δ是资本折旧率s是储蓄率。代码上就是简单的递推for (t in 1:(t_max - 1)) { K[t 1] - (1 - delta) * K[t] s[t] * net_output[t] }在这一步我踩过的一个坑是储蓄率s被设定为控制变量后优化过程容易让储蓄率冲击上下边界导致路径震荡。后来我参考诺德豪斯原文的做法把储蓄率处理成常数或仅在有限范围内浮动模型的稳定性显著提升。2.3 排放模块的实现排放模块的逻辑并不复杂总排放等于工业排放加上土地排放其中工业排放由经济活动驱动E_Ind(t) σ(t) * Q(t) * (1 - μ(t))σ是碳排放强度表示单位产出对应的碳排放量它会随时间下降技术进步μ是减排率代表从基准路径中削减的排放比例。calc_industrial_emissions - function(sigma, net_output, mu) { sigma * net_output * (1 - mu) } calc_total_emissions - function(e_industrial, e_land) { e_industrial e_land }土地排放是一个外生序列早期DICE版本把这个量设为常数或缓慢下降的路径可以直接用一条预设的序列填充。这部分有一个容易忽略的细节公式中的使用产出到底是毛产出还是净产出。诺德豪斯的原始模型中排放是基于毛产出来算的也就是说减排成本本身也是“有碳排放”的活动忽略这一点会低估总排放。我自己实际对比过用净产出算排放和用毛产出算排放2100年温度结果会差零点几度这个误差在政策讨论中并不小。2.4 碳循环与气候模块的实现碳循环模型采用三层碳库结构大气碳、浅层海洋碳和深层海洋碳之间以线性转移系数相互交换。我写成代码后最直观的感受是这就是一个三阶线性系统calc_carbon_cycle - function(MAT, MUP, MLO, E_total, phi) { dt - 1 new_MAT - phi[1,1] * MAT phi[1,2] * MUP E_total * dt new_MUP - phi[2,1] * MAT phi[2,2] * MUP phi[2,3] * MLO new_MLO - phi[3,2] * MUP phi[3,3] * MLO c(new_MAT, new_MUP, new_MLO) }这里phi是3乘3的转移矩阵各元素表示碳库之间的年转移比例。这些系数取值在诺德豪斯的开源代码和论文附录里都有我建议直接照抄原始校准结果不要自己去拟合否则容易得到不稳定的系统。气候模块先将大气碳浓度转化为辐射强迫F(t) η * log2(MAT(t) / MAT_1750) F_EX其中η是辐射强迫系数MAT_1750是工业革命前的大气碳浓度F_EX是其他温室气体的外生辐射强迫。然后温度方程是这个样子calc_temperature - function(F, TAT, TLO, lambda, c1, c2, c3) { new_TAT - TAT c1 * (F - lambda * TAT - c2 * (TAT - TLO)) new_TLO - TLO c3 * (TAT - TLO) c(new_TAT, new_TLO) }温度方程是整个模型中最容易出现数值震荡的地方。我调试时有一次把时间步长从10年改成1年重新模拟发现温度路径出现了周期摆动。排查后发现是因为参数是从以10年为单位校准的直接缩小时间步长而不重新校准参数就不稳定。这也是一个典型“不要随便改模型尺度”的教训。2.5 福利模块与目标函数福利模块把所有时期的人均消费效用按折现率加总。DICE模型使用对数效用函数W Σ β^t * L(t) * ln(C(t) / L(t))其中β 1 / (1 ρ)ρ是纯时间偏好率。在R中目标函数可以写成一个数组乘以权重的加权和calc_welfare - function(consumption, population, rho, t_max) { beta - 1 / (1 rho) years_idx - 0:(t_max - 1) utility - population * log(consumption / population) welfare - sum(beta^years_idx * utility) * step welfare }这里我对每一个核心模块在优化循环内的写法做了整理一个完整的DICE模拟函数可以封装成传入减排率和储蓄率序列、返回福利值和各状态变量的形式。有了这个封装后面的优化和情景分析就方便很多。3. 减排路径的量化设计从政策蓝图到模型参数DICE模型的真正价值不在于跑一个基准情景而在于回答“如果这样减排会怎样”。但“减排路径”这个说法什么都解释不了必须转成模型可识别的量化参数。这一节我重点讲怎么把不同的减排路径落实到模型代码里。3.1 四种典型减排路径的模型表达在DICE框架下减排路径主要通过减排率μ的路径和碳税来体现。我把实操中常用的四类路径整理成了一个表方便对照情景名称核心设定模型表达方式典型用途基准情景BAU不采取额外减排措施μ恒定为0或极小值作为对照基线碳税情景设定碳税水平由碳税反推μ或用碳税作为影子减排约束评估碳税政策效果碳强度约束设定排放强度下降速度每期给定最大μ使排放强度沿预设斜率下降模拟碳中和目标温控目标反推温升不超过给定上限在优化中加入温度约束反解最优μ路径评估2度/1.5度目标可行性BAU情景最省事把μ序列设成全0就可以。但这里有个细节即使μ为0排放强度σ仍然会因技术进步而下降所以BAU情景下碳排放并不是一条直线上升而会有一定程度的自动减速。碳税情景的实现要复杂一些。在DICE框架内可以直接把碳税当作减排的影子价格通过模型的边际减排成本函数反推减排率碳税等于边际减排成本时就是最优减排率。我的实现方式是先计算边际减排成本曲线再根据给定碳税水平插值得到对应减排率。碳强度约束的实现我以“2030年碳达峰、2060年碳中和”这个目标做了个示例。把这个目标分解成模型语言2030年之前排放强度下降速度保持历史趋势2030到2060年排放总量从峰值逐年线性下降至接近零对应地每年需要达到的减排率可以反推出来# 反推满足碳中和目标的减排率 calculate_mu_carbon_neutral - function(e_baU, target_emissions) { mu_needed - 1 - target_emissions / e_baU mu_needed[mu_needed 0] - 0 mu_needed }3.2 关键参数的校准与敏感性DICE模型的参数很多但真正对结论影响最大的是这么几个气候敏感度平衡升温对二氧化碳加倍的响应、减排成本函数的系数π和θ、纯时间偏好率ρ、以及资本产出弹性γ。诺德豪斯原文给出的参数值可以作为基准参数符号典型取值影响方向资本产出弹性γ0.3越高则减排对产出影响越大资本折旧率δ0.1/10年影响资本积累速度纯时间偏好率ρ0.015/年越高则越不重视未来损失气候敏感度S3.1摄氏度越高则损害越严重减排成本指数θ2.6-3.0越高则减排成本越大我强烈建议你拿到模型后先做一轮参数敏感性分析只变动一个参数其他保持不变然后观察温度、GDP、福利这三个关键结果的变化幅度。我做了这个分析后发现ρ和气候敏感度对最优碳价的影响非常大这也是为什么不同机构用同一套DICE框架会得出不同政策建议的根本原因。3.3 用R代码实现情景切换为了让情景切换灵活我在代码里设计了一个参数列表结构每次跑模拟前只需要更新这个列表就行scenario_params - list( name carbon_tax_100, carbon_tax 100, # 美元/吨CO2 mu_path NULL, # 若为NULL则用碳税反推 max_temp_cap NULL, # 温控目标如2.0 start_year 2025 )然后写一个统一的run_dice函数根据情景参数生成减排路径后调用模拟。这个设计的好处是后面对比十几种情景时每个情景只需要一行参数配置不需要改动模型主体代码。这也是一个典型的“把变化的部分和不变的部分分开”的工程化思路。4. 模拟结果解读经济增长、碳排放、温度与福利的权衡模型跑完之后最重要的环节是结果解读。这里我用一组教学演示参数跑了四种情景下面要讨论的数值请作为方法论参考而不是权威预测。但解读数据的方式方法是完全通用的。4.1 碳排放轨迹峰值的“尖”与“钝”我分别跑了无减排BAU、碳税100美元、碳中和目标和2度温控目标四个情景把它们在2100年之前的碳排放路径放在一起对比。BAU情景的全球碳排放一路上升2070年左右才因技术自然进步而触顶2100年排放水平大约是工业排放的2倍左右。碳税100美元的情景下排放曲线明显下压2080年以后开始持续下降但峰值出现得并不早——碳税虽是逐年递增的一开始税低时效果有限。碳中和目标这条路径很有意思排放轨迹有一个非常尖锐的顶点因为“2030达峰、2060中和”相当于人为设定了一个陡峭的反转结构。2度温控目标反推出来的路径则更平滑因为加了跨期优化后模型会把减排压力部分前置避免后期“急刹车”。我有一个判断路径形状比某个单点峰值更重要。尖峰意味着转型冲击集中在短时间内资本折旧、就业转移、技术替代的压力都会放大。如果分析政策可行性不能只盯着“能不能达峰”还要看峰值前后十年的变化斜率。4.2 温度变化2100年升温的“终局”温度是气候模型最直观的输出。BAU情景跑出来的2100年升温大约是4度以上碳税100美元约3.2度碳中和目标约2.4度2度温控目标则真的可以压在1.9度左右。这个顺序符合直觉减排力度越大终局温度越低。但需要注意一点DICE模型里的温度方程是包含海洋热惯性的即使碳排放已经下降温度仍然会继续爬升一段时间直到辐射强迫和热吸收达到平衡。看温度曲线时不要只看“什么时候达峰”还要看“峰后回落的速度”。碳中和情景下温度大约在2080-2090年之间触顶之后缓慢回落这个“滞后顶”是气候系统的惯性决定的。4.3 GDP代价与福利变化钱花得值不值接下来就是政策讨论中最核心的模块了。每个减排情景都要付出GDP代价——减排成本直接挤占当期消费和投资。我把四种情景相对BAU的累计GDP损失折现到现值算出来看碳税100美元情景的累计GDP损失大约在1.5%左右碳中和目标大约2.8%2度温控目标大约3.5%。这些数字单独看很抽象但换算成绝对量就是几万亿到几十万亿美元的规模。此时福利函数的作用就出来了。DICE模型的社会福利是人均消费效用的折现加总它会同时惩罚“消费下降”和“未来气候损失”。某些情景下GDP损失高但避免了更大的气候损害福利反而更高。这正是DICE模型最有教育意义的输出它让“减排成本”和“气候损害”在同一个度量下比较。我跑下来的结果显示碳税100美元情景相对BAU的福利是提升的因为减排带来的气候损失减少超过了减排成本本身。这不是说碳税越高越好——如果直接跳到500美元每吨过高的成本会让福利重新恶化。福利最大的区域在中间某处这就是模型给出的“最优碳价区间”。4.4 用ggplot2做出能直接进报告的结果图结果可视化我建议直接用ggplot2完成。核心就是整理出每个情景的温度、排放、GDP的逐年数据框然后用分面图把四张图拼在一起library(ggplot2) library(dplyr) results_df %% ggplot(aes(x year, y temp_at, color scenario)) geom_line(size 1.1) theme_minimal() labs( title 不同减排路径下的全球升温模拟, x 年份, y 升温摄氏度, color 情景 )我实际经验中最难处理的不是绘图本身而是把模型输出的长格式数据准备好。模拟函数每跑一次会生成几十列数据我建议在封装模拟函数时就顺便把结果整理成“长表”每行是“情景-年份-指标-数值”这样绘图和后续统计分析会省掉大量数据处理的麻烦。5. 实操中踩过的坑与改进思路每次做R语言DICE模型实操多多少少都会在几个地方卡住。这里我把自己踩过的坑和一个改进思路整理出来算是给准备入手的你一点预警。5.1 优化算法不收敛从初始值到罚函数第一版代码用的是optim()函数优化减排率和储蓄率序列几十个变量同时优化经常遇到不收敛或收敛到局部最优的问题。我一开始以为是算法设置不对后来才发现是初始值距离最优路径太远迭代过程中目标函数值出现NaN。解决办法是做一个“两阶段优化”第一阶段用遗传算法在全局范围内粗搜索找到近似最优的减排率路径第二阶段再用梯度类算法从该初始值出发做精细优化。另一个有效技巧是给控制变量加边界约束用L-BFGS-B方法并显式指定lower和upper。对减排率μ设定0到1的区间对储蓄率s设定0.05到0.4的区间可以避掉一大半数值问题。5.2 折现率的选取模型里最“脆弱”的旋钮折现率这个话题我在实操中觉得它已经不止是技术问题而是一个哲学问题。把纯时间偏好率从每年0.015降到0.001模拟出的最优碳价会翻倍还不止。低折现率意味着更“在乎”子孙后代的福利自然会得出更激进的减排路径。我的处理方式是不要纠结于“哪个折现率是对的”而是做区间分析分别用0.5%、1.5%和3%跑一组情景把结果作为一个区间呈现给决策者。模型的价值在于揭示敏感性和权衡关系而不是给出唯一“真理”。5.3 参数不确定性从点估计到分布模拟DICE模型输出的所有结果都是建立在参数点估计之上的。但参数本身就有不确定性比如气候敏感度的置信区间很宽减排成本函数的形状参数也是估计出来的。更稳健的做法是给关键参数设定先验分布然后用蒙特卡洛模拟跑几百组参数组合得到结果的概率分布。这一步在R里实现成本很低用mvtnorm包生成参数向量然后并行跑模拟就行。我跑过500组参数组合普通笔记本大概需要十几分钟。输出不再是单条温度曲线而是一个扇形图——中位数曲线加上置信区间这样的图放到政策报告里会专业得多也会回答很多来自决策者的质疑。5.4 代码架构从单人研究工具到多人协作项目最后一条是我自己的经验教训。最早的模型代码是一个300行的脚本文件函数和全局变量混在一起改一个参数要搜索三四个地方。后来因为要加情景对比、敏感性分析、蒙特卡洛模拟我重构了代码结构按功能拆成了几个模块dice_params.R参数定义与默认取值dice_model.R模型主体函数包括经济、碳循环、气候、福利模块dice_scenarios.R情景设定与路径生成dice_optim.R优化部分包括目标函数和求解设定dice_plot.R结果可视化这种结构未必是标准答案但把“参数-模型-情景-优化-可视化”五层分清楚之后代码复用率高了很多新接手的同事理解起来也快。如果只是自己跑一两次实验脚本式写法没问题但只要想长期用这个模型做跟踪分析代码结构越早理清越好。写在最后一点关于模型边界的提醒跑通了DICE模型之后你会越来越清晰地感受到这个模型的边界它的简洁既是优点也是局限。全球被压缩成一个“代表性地区”没有区域差异消费者是完全理性的跨期最优没有行为摩擦损害函数是总产出的一个比例但真实世界的气候损害是非线性的、不均匀的、可能带突变风险的。这些局限并不等于模型没有用关键在使用者是否清醒地知道模型在哪里简化了现实。我的建议是把DICE模型当作“战略级”的思考工具而不是“战术级”的精准预测工具。它最适合回答的是全局性的比较问题——哪个政策方向更优、碳价区间大致在什么范围、时间路径应该前置还是后置。如果要做区域性分析或行业影响分析DICE这套框架需要拓展成区域化模型或IAM更细致版本那已经是另外一个层次的工作了。最后分享一个实操小技巧跑出结果后先不说结论把各个情景的碳价、温度、GDP损失三个指标列成一张表然后请团队里最毒舌的同事来“审问”这张表。如果他能问出“为什么这个情景的温度下降比排放下降慢那么多”“为什么碳税情景的GDP损失比碳中和情景低”这类问题说明你的结果已经足够引发深度讨论而这才是DICE模型真正应该发挥的作用。