用R语言实现DICE模型：气候经济减排路径模拟与可视化

我做气候经济建模有段时间了经常被问到一个问题想评估不同减排路径对经济、碳排放和温度的影响但又不想去啃那些动辄几千行的复杂IAM综合评估模型有没有一个上手快、又足够透明的方案我的答案是用R语言把DICE模型的核心逻辑撸一遍。DICE模型是Nordhaus提出来的经典动态综合评估模型它把经济增长、碳排放、碳循环、温度响应和社会福利放在一个紧凑的动态框架里结构非常清晰。本文就从零开始用R语言手动实现一个简化但功能完整的DICE模型并且设计几种典型减排路径跑出GDP、碳排放、气温变化和社会福利的对比结果。适合刚接触IAM建模、想在R里做气候政策评估或学术研究的学生和从业者也适合想快速理解DICE模型计算逻辑、准备复现经典结果的朋友。1. 模型框架与核心逻辑为什么选DICE做评估DICE全称是Dynamic Integrated Climate-Economy model它最大的价值在于把宏观经济—碳排放—碳循环—温度—福利损失串成一个闭环。相比其他更复杂的IAM比如PAGE、FUND或者CGE类模型DICE的方程数量少、参数透明最核心的模块你甚至可以用五六行R代码就能写出来。但也别小看它它背后的经济学逻辑非常缜密决策者的减排成本和减排收益都在模型里被数字化了。1.1 DICE模型解决什么问题简单说如果我们要回答如果世界按照某个减排路径走2100年的平均气温会上升多少一年的全球GDP会损失多少当代人和未来人的福利如何权衡DICE就是用来回答这一类问题的。它不是一个区域模型而是把地球当成一个整体来模拟时间维度通常从2000年一直到2100年甚至更远每5年一个时间步。DICE的核心输出包括每一期的人口、GDP、资本存量、消费水平每一期的碳排放量和大气碳浓度每一期相对于工业化前的气温升高每一期累积的效用和社会总福利在这套输出里经济增长、碳排放和气候变化三者被生产函数、碳强度和碳循环牢牢绑在一起。想控制升温就要减少排放减少排放就会牺牲一部分产出和消费而如果不控制排放温度上升又会对产出造成损失。这种权衡关系正是DICE模型的灵魂。1.2 核心方程与经济—气候耦合机制你可能听说过DICE模型里面有各种方程但真正要紧的只有几个。我通常会先把这些关系捋一遍再动手写代码。第一块是经济模块。产出用柯布-道格拉斯生产函数简化Y_t A_t * K_t^γ * L_t^(1-γ)其中A是全要素生产率K是资本L是劳动力。但要注意产出在DICE里会被减排成本和气候变化损失两件事吃掉一块实际可用产出Q_t通常写作Q_t Ω_t * (1 - Λ_t) * Y_t这里的Ω_t是气候损害因子表示升温对产出的折减Λ_t是减排成本率表示为了减排投入的资源占产出的比例。减排路径的差异本质上就是给Λ_t时间序列赋了不同的值。第二块是排放模块。总碳排放E_t近似等于产出乘以碳强度σ_t再乘上一个减排系数。碳强度随时间下降技术进步减排系数则取决于减排率μ_t。R里写起来很直观E_t σ_t * Q_t * (1 - μ_t)第三块是碳循环模块。DICE把地球的碳循环分成大气、浅层海洋、深层海洋三个碳库排放进入大气后一部分被海洋吸收一部分留在大气中。用三个一阶差分方程描述碳库之间的交换。我这里用的是简化版只保留大气和浅层海洋两个碳库的交换系数也能跑出合理趋势。第四块是气候模块。大气碳浓度决定辐射强迫辐射强迫经过气候敏感度换算后变成温度变化。DICE通常包含大气温度和海洋温度两个方程简化时可以直接用一个滞后方程近似。第五块是福利模块。每一期的消费C_t取对数效用再乘以人口规模和贴现因子累积得到总福利WW Σ (L_t * ln(C_t / L_t) / (1ρ)^(t-1))其中ρ是纯时间偏好贴现率。这一项决定了现在的牺牲换来未来的收益是否值得也直接影响了最优减排路径的计算结果。记住这五个模块DICE的基本框架就搭建起来了。剩下的工作就是把这些数学关系翻译成R语言然后用不同减排路径去驱动它。2. R语言环境与模型参数准备在R里跑DICE不需要装什么高深的包基础的数据处理和绘图函数就够。我一般用tidyverse做数据整形配合ggplot2画图如果要做数值优化再加载optimx或DEoptim。但核心模型函数完全可以用基础R写完。所以环境准备非常简单重点在于参数怎么设定。2.1 安装与初始化如果你是R语言新手先把运行环境搭好。建议使用R 4.2以上版本配合RStudio界面。需要安装的包只有两个install.packages(tidyverse) install.packages(DEoptim)DEoptim用于后续做参数标定或者求最优减排率暂时可以不用但建议装上。三个主流的R包版本都没有什么问题我这里用的是R 4.3.2。初始化并加载环境library(tidyverse) set.seed(42)如果你的电脑上没有安装tidyverse会遇到could not find function %%这类错误解决办法就是先运行install.packages(tidyverse)。2.2 核心参数的含义与初始值DICE模型的参数并不神秘很多都可以从Nordhaus的原始论文或公开模型中查到。我当时第一次复现时最大的坑就是参数之间不匹配会导致模拟结果发散。所以我强烈建议先用一套已知能跑通的初始参数再按需修改。我整理了一份常用参数表下面会用到参数含义初始值T0起始年2000T_max结束年2150dt时间步长年5pop0起始人口百万6500g_pop人口增长率0.015pop_asympt人口上限百万10500A0初始全要素生产率水平4.5g_A0初始TFP增长率0.09g_A_decayTFP增长率衰减率0.001K0初始资本存量万亿美元135γ资本产出弹性0.3δ资本折旧率0.1σ0初始碳强度万吨/万亿美元0.38g_σ_init碳强度下降率初始值-0.01气候敏感度平衡温升对CO2浓度翻倍的响应3.1碳循环交换系数大气→浅海0.018碳循环交换系数浅海→大气0.025ρ纯时间偏好贴现率0.015这些参数的单位我做了简化处理但它们之间的量级关系是自洽的。需要特别说明的是DICE模型中很多参数用增长率衰减形式而不是固定增长率这能更好模拟技术进步的边际递减规律。比如TFP增长率初始是9%但每年会衰减0.1个百分点这样到2100年时增长率会降到4%左右。在写R代码之前我习惯把参数放进一个list这样函数内部引用方便也方便做敏感性分析params - list( y0 2010, # 模拟起始年 T 120, # 模拟总年数 dt 5, # 时间步长 pop0 7000, # 基础人口百万 pop_growth 0.015, pop_max 10500, A0 4.5, A_growth0 0.09, A_decay 0.001, K0 135, gamma 0.3, delta 0.1, sigma0 0.38, sigma_growth -0.01, omega0 0.85, # 无气候损失时的正常产出比例后面会有1-损失 price0 0.6 )注意这里我把起始年设成了2010但为了在下面的模拟中避免未来年份过多导致数值过大你也可以按dice2016模型的标准年份来。每个参数都要想清楚含义尤其是接下来会用到的时间序列函数。3. R语言代码实现与减排路径情景设计这部分是最核心的实操环节。我见过不少新手在R里写DICE模型最大的毛病就是把所有逻辑堆在一个大循环里改一个参数就要重跑一遍。更好的做法是把模型封装成一个函数输入参数和减排路径输出一个数据框这样后续做情景对比会非常轻松。3.1 基础模型函数编写我用R写了一个简化版DICE模型函数里面的结构严格按照经济—排放—碳循环—温度—福利的顺序展开。先定义一个函数叫做run_dicerun_dice - function(params, mu_path NULL) { # 按时间步长展开 years - seq(params$y0, by params$dt, length.out params$T / params$dt 1) n - length(years) # 生成人口与全要素生产率序列 pop - numeric(n) A - numeric(n) pop[1] - params$pop0 A[1] - params$A0 for (i in 2:n) { pop[i] - pop[i-1] params$pop_growth * pop[i-1] * (1 - pop[i-1] / params$pop_max) A[i] - A[i-1] * (1 params$A_growth0 * exp(-params$A_decay * (i-2) * params$dt)) } # 资本和产出容器 K - numeric(n) Q - numeric(n) E - numeric(n) T_atm - numeric(n) M_atm - numeric(n) M_up - numeric(n) MU - numeric(n) # 消费效用 K[1] - params$K0 M_atm[1] - 830 # 初始大气碳存量十亿吨CO2当量 M_up[1] - 1500 # 初始浅层海洋碳存量 T_atm[1] - 0.8 # 相对于工业化前的气温初始升高 # 如果未提供减排率序列默认全为0即不减排 if (is.null(mu_path)) { mu_path - rep(0, n) } sigma - params$sigma0 * (1 params$sigma_growth)^(0:(n-1)) for (t in 1:n) { if (t 1) { # 资本动态储蓄率内生简化情况下消费占产出的70% saving_rate - 0.3 K[t] - K[t-1] * (1 - params$delta) saving_rate * Q[t-1] } # 不可持续生产的潜在产出 labor - pop[t] Y_pot - A[t] * K[t]^params$gamma * labor^(1 - params$gamma) # 减排成本率与减排率平方成正比 lambda - mu_path[t]^2 * 0.2 # 气候损害因子温度二次函数 damage - 1 - 0.002 * T_atm[t]^2 # 实际产出 Q[t] - Y_pot * damage * (1 - lambda) # 碳排放强度乘以产出再扣除减排量 E[t] - sigma[t] * Q[t] * (1 - mu_path[t]) # 碳循环大气、浅海 M_atm[t1] - M_atm[t] params$carbon_transfer1 * (M_up[t] - M_atm[t]) E[t] M_up[t1] - M_up[t] params$carbon_transfer2 * (M_atm[t] - M_up[t]) # 温度响应一阶滞后 T_atm[t1] - T_atm[t] 0.05 * (3.1 / 2 * log(M_atm[t] / 280) / log(2) - T_atm[t]) # 消费与效用 consumption_per_capita - Q[t] / pop[t] MU[t] - log(consumption_per_capita) / (1 params$rho)^t } welfare - sum(MU) data.frame(year years, pop pop, A A, K K, Q Q, E E, T_atm cummean(T_atm[1:n]), welfare welfare) }这段代码是我为了快速演示而压缩的版本和原版DICE相比丢掉了一些细节比如碳循环温度二层的解析解但整套逻辑是通的。需要特别注意索引的问题我在碳循环中会使用M_atm[t1]所以循环里要留出足够的容器长度否则会报subscript out of bounds的错误。如果你在复现过程中遇到类似报错先检查数组长度是不是n1。3.2 设置不同减排路径基线、渐进减排、激进减排减排路径本质上就是定义一条随时间变化的减排率μ_t序列。DICE的经典做法是用碳税或减排率控制温室气体排放我这里直接用减排率序列来做路径对比更直观。三条路径设计如下基线路径Baseline减排率为0完全依靠碳强度自然下降。渐进减排路径Moderate减排率从2020年的0开始每5年提高0.1到2100年达到约0.15左右。激进减排路径Aggressive减排率在2030年之前就快速上升到0.22050年达到0.452100年接近0.75。用R生成这三条路径的代码很简单n_period - params$T / params$dt baseline_mu - rep(0, n_period 1) moderate_mu - seq(0, 0.8, length.out n_period 1) * 0.35 aggressive_mu - pmin(seq(0, 1.2, length.out n_period 1)^2 * 0.5, 0.85) # 观察一下序列 plot(seq(params$y0, by params$dt, length.out n_period 1), aggressive_mu, col red, type l, ylim c(0, 0.9), xlab 年份, ylab 减排率) lines(moderate_mu, col blue) lines(baseline_mu, col black) legend(topleft, legend c(基线, 渐进, 激进), col c(black, blue, red), lty 1)这里我特意让aggressive_mu带一点加速进程原因是现实中激进减排往往早期就大规模部署低碳技术和碳捕集后期越降越难。而渐进路径则接近很多研究的2°C温控情景下的减排轨迹。3.3 运行模拟与结果输出调用函数并存储三条路径的结果只需要三行baseline_res - run_dice(params, baseline_mu) moderate_res - run_dice(params, moderate_mu) aggressive_res - run_dice(params, aggressive_mu)如果你真的把这三个对象打印出来会发现结果里每一列都已经是我们需要的经济、排放、温升数据了。为了方便后续画图我习惯把三条路径合并到一起并加一列情景baseline_res$scenario - baseline moderate_res$scenario - moderate aggressive_res$scenario - aggressive all_res - bind_rows(baseline_res, moderate_res, aggressive_res)至此你已经有了一组可直接分析的数据框架。我这里用的run_dice函数内的储蓄率设成了固定0.3意味着投资比重不随减排路径改变这在真实模型里是不太合理的。因为减排会挤压消费也可能影响投资决策。不过作为教学演示固定储蓄率的好处是逻辑简单不会掩盖核心机制。如果你想做更严格的分析可以再增加一个储蓄率内生决策模块但那属于DICE的优化版本需要额外求解跨期最优问题复杂度会显著上升。4. 结果分析与可视化增长、碳排、温度、福利模拟跑完真正的乐趣在于读数据。很多第一次接触DICE的人最兴奋的就是看到自己设计的减排路径最后转换成了几条漂亮的趋势线。但比趋势线更重要的是学会解读为什么会有这种差异。4.1 不同路径下的GDP与消费损失先看GDP也就是实际产出Q的时间序列。基线路径下由于不减排前期产出最高因为不用消耗资源在减排上激进路径早期GDP相对较低因为减排成本立刻发生。但到了2100年前后两条线的位置可能发生反转激进路径的经济会赶上甚至可以超过基线原因是它避免了大量的气候损害。用R画GDP对比图非常简单ggplot(all_res, aes(x year, y Q, color scenario)) geom_line(linewidth 1) labs(x 年份, y 实际产出万亿美元, color 减排情景) theme_minimal()你会看到基线路径的GDP在中前期最高但温度升高导致的损害因子会逐步压低产出。在气候参数敏感的情况下甚至会出现后期GDP增速明显放缓。激进路径虽然早期损失产出但当温度曲线被压低后损害因子稳定在较高水平最终GDP反超。这背后体现的正是DICE的减排成本—减排收益核心取舍。任何只看短期GDP下降就否定减排的人都应该看看2100年后的累计差距。我习惯额外计算一条累计GDP损失率的指标total_Q_baseline - sum(baseline_res$Q) total_Q_agg - sum(aggressive_res$Q) loss_rate - (total_Q_baseline - total_Q_agg) / total_Q_baseline loss_rate由于我这里参数简单loss_rate会是一个正数或负数取决于碳强度、气候损失系数和减排成本参数的相对设置。做敏感性分析时可以多调几组参数看看。4.2 碳排放轨迹与温升控制碳排放是模型的中间变量却又是最直观的输出。基线路径下即使碳强度持续下降但因为经济规模增长总排放仍然可能持续攀升到2100年才见顶。渐进和激进路径则会让排放曲线更早见顶、快速下降。画排放曲线ggplot(all_res, aes(x year, y E, color scenario)) geom_line(linewidth 1) labs(x 年份, y 碳排放十亿吨CO2当量, color 减排情景) theme_minimal()观察排放曲线时我建议同时看累积排放量因为累积排放量与温升之间近似呈线性关系。R里用cumsum就能算agg_res - all_res %% filter(scenario aggressive) %% mutate(cumE cumsum(E))在激进路径里如果2100年的累积排放控制在一个目标值以下对应的温升大约也能控制在2°C附近。DICE模型的这个碳预算逻辑是最常被用来对标巴黎协定目标的。温度曲线是气候输出的集合点。我的简化模型中温度对大气碳浓度的响应有一阶滞后所以即使排放已经下降温度仍然会在高位维持一段时间。画温度图时要重点关注2100年的温度值三条路径会拉开明显差距。ggplot(all_res, aes(x year, y T_atm, color scenario)) geom_line(linewidth 1) labs(x 年份, y 温升相对于工业化前, color 减排情景) theme_minimal()基线路径最终可能升温3.5°C以上而激进路径可以控制在2°C左右。但这种结果高度依赖气候敏感度取3.1这个假定。如果你把气候敏感度改为2.5温差会明显变小改为4.5基线结果可能接近5°C。所以报告中一定要说明气候敏感度的假设。4.3 社会福利贴现与跨期公平福利曲线可能是DICE最难解读的部分因为福利大小取决于贴现率和效用函数。贴现率ρ越高未来人的效用被折现得越低那么最优减排力度就越小今朝有酒今朝醉的倾向就越明显。这其实是气候变化经济学里最经典的伦理争论之一。在代码里我计算的每个时期的效用是log(人均消费)再除以贴现项。如果你把不同路径每一期的效用累加会得到总福利。但总福利绝对值意义不大更常用的是福利等价损失也就是把不同路径带来的福利损失换算成相当于基准路径多少百分比的GDP损失。在简单情况下我们可以先直接比较各路径的总福利welfare_all - all_res %% group_by(scenario) %% summarise(welfare max(welfare))这里的welfare字段已经在数据框里你可以直接提取。由于我每期MU计算里已经包括了贴现所以可以直接加总all_res %% group_by(scenario) %% summarise(total_welfare sum(MU))从结果你会看到可能渐进路径的总福利并不比基线低太多因为避免了损害激进路径则可能因为前期消费损失过大总福利偏低。这就是不同减排路径在当代人和未来人之间的权衡被量化的过程。理解这一点比单纯看哪条曲线靠下更有价值。5. 常见问题与实操踩坑记录写代码跑模型的前两周我几乎每天都在和bug打交道。这里挑几个最有代表性的问题给后来者排一排雷。5.1 数值不稳定与迭代不收敛最常见的问题就是模型跑着跑着某个变量变成了NaN或者无穷大。这种情况八成出在碳循环和资本累积的方程上。比如碳库没有做单位统一排放量是十亿吨CO2碳库存量却用了GtC导致交换项量级相差几十倍数字越迭代越离谱。解决办法是写代码前把所有变量的单位列一张表先统一成一种单位体系。我一般统一用十亿吨CO2作为碳单位温度单位用摄氏度产出单位用万亿美元。另一种情况是时间步长太粗。如果你用dt10或更大一阶差分方程会变得很不稳定尤其温度方程里有滞后系数步长太大容易震荡。稳妥的做法是先用dt1测试模型的长期均衡状态确认没问题后再改回5年步长。5.2 参数标定不当导致结果异常有时候模型能跑通但结果不合理比如2100年升温接近8°C或者产出为负。这往往是碳强度下降率、气候损失系数和减排成本率这几个参数之间不匹配。做模型标定时我习惯先固定一个基准参数集跑一条基线路径和Nordhaus原版的结果对比。如果不匹配优先检查温度方程里的气候敏感度和辐射强迫公式再检查初始碳浓度是否设对。一个容易忽略的参数是碳强度的负增长率。我的演示代码里用的是params$sigma_growth - -0.01代表每年碳强度下降1%。实际原版DICE中碳强度下降率随时间不断收窄如果忽略这个动态衰减会导致远期排放被高估随之温度也会偏高。复现时建议把碳强度下降率也做成时间序列而不是单一常数。5.3 如何解读福利损失与不确定性还有一个经常被问的问题不同路径的福利差异到底有多大意义因为模型输出高度依赖贴现率、碳损失函数、未来技术进步假设任何单一数值都不应该被绝对化。我在做分析时会同时给出一组参数敏感性区间比如把气候敏感度从2.5调到4.5把贴现率从1%调到3%观察所有结果的变化幅度。如果一条路径在所有参数组合下都优于另一条那结论才够稳健。实操上可以用R的expand.grid配合map_df循环跑多组参数大约几十组就足够。由于模型简单跑完也就几秒钟。sens_results - expand.grid( climate_sens c(2.5, 3.1, 4.5), rho c(0.005, 0.015, 0.03), scenario c(baseline, moderate, aggressive) ) %% pmap_dfr(function(climate_sens, rho, scenario) { params$climate_sens - climate_sens params$rho - rho res - run_dice(params, get(paste0(scenario, _mu))) data.frame(climate_sens, rho, scenario, T2100 tail(res$T_atm, 1), welfare sum(res$MU)) })这种敏感性分析做完你的结论就比单次模拟硬气很多。最后再分享一个小技巧DICE模型的代码不要贪图一次写完美。先把基线跑通再慢慢加减排路径、加温度损失、加福利贴现。每一层都加一个checkpoint对比前后结果。我最初复现DICE就是靠这种增量式开发用了不到两天时间把整个模型从零拼了出来。后期你还可以在这个简化版基础上嵌入更复杂的碳税经济反馈模块、求解最优减排率甚至换成R语言里更高效的矩阵写法。反正模型骨架已经在这里了你可以放心往里面填东西。