简介：一套基于投影方法求解DSGE模型的Matlab实现方案，面向经济学、金融学及数学方向的研究生、科研人员，以及需要完成相关课程设计或毕业设计的本硕学生，适合具备一定宏观经济学与Matlab基础的读者。方案内含完整的参数化编程框架，参数可灵活调整，附带可直接运行的案例数据，可用于复现或扩展实际模型。压缩包共107个文件，其中86个m脚本覆盖主程序、求解函数与工具函数，另有txt说明、mat数据、pdf文档及mod模型文件等，整体约3.16MB，目录结构清晰，便于按模块检索和二次开发。投影方法的Smolyak基函数、投资子代生成等关键步骤均有对应代码实现，注释明细，可帮助理解投影法求解DSGE的核心逻辑。目前已有55人学习，适合希望掌握Dynare之外另类求解工具的学习者参考。
1. 非线性 DSGE 为什么需要投影方法：从 Dynare 的线性近似说起
如果你已经跑腻了 Dynare 里的标准 RBC，偶尔想加点大波动、零利率下限或者偶尔约束，大概率会遇到同一个尴尬：默认的一阶近似在稳态附近很好用，但状态变量一旦走远，脉冲响应和决策函数就变成了一条直线，非线性被削没了。投影方法换了一个思路——它不再围绕稳态做局部泰勒展开，而是在整个可能到达的状态空间上，用一个全局多项式把决策函数拟合出来，然后通过最小化欧拉方程残差来锁定多项式系数。Promes 工具箱就是围绕这个思路搭建的一套 Matlab 实现，目录里能看到 Smolyak 插值、各向同性网格和三次多项式基底，也带了可以直接跑的 Galizia 模型和住宅模型入口。接下来我会从文件结构开始拆，逐步讲到投影求解循环里最关键的参数和收敛问题，最后告诉你这种工具箱怎么搬到自己的 DSGE 模型上。
2. Promes 工具箱的文件骨架：模型描述、多项式基与求解入口
从目录中一眼能看到的文件说起，整个工具箱并不大，核心计算集中在几个 m 文件和两个以.jnl、.log结尾的文本文件中。.jnl后缀在 Matlab 里不是标准扩展名，Promes 把它当作结构化的模型描述文件读取；.log是上一次求解留下的运行日志，里面记录了迭代序号、残差范数和完成状态。第一次跑模型前，建议先打开stnd_rbc_dyn.log扫一眼，确认它不是从一次异常中断中保存下来的，否则复现时一上来就被旧状态误导。
2.1 模型描述文件 stnd_rbc_dyn.jnl 里究竟放什么
.jnl文件在 Proxmox？不是，在 Promes 里承担的是“模型配置”的角色：状态变量个数、冲击个数、最大多项式阶数、迭代容差和参数初值都从这类文件读入。下面给出一个接近 Promes 惯例的可读示例，保存为stnd_rbc_dyn.jnl时会用行首#做注释，执行时可被parse_jnl这类读取器跳过。
# 模型参数段 model_name = stnd_rbc_dyn n_state = 2 n_shock = 1 max_degree = 3 iteration_tol = 1e-8 alpha_param = 0.36 beta_param = 0.99 gamma_param = 2.0 rho_a_param = 0.95 sigma_a_param = 0.01这种设计的好处是，你不需要为了调一个参数频繁改动 Matlab 主脚本。alpha_param是资本产出弹性，beta_param是贴现因子，gamma_param是相对风险厌恶系数，rho_a_param和sigma_a_param描述生产率冲击的持久性与波动幅度。修改参数只在文本文件里完成，主程序通过统一的配置读取函数载入。
2.2 核心 M 文件的职责分工
用type solve_proj可以查看主求解器，但在此之前先梳理清楚文件之间的调用关系。下表把目录中几个关键文件对应到投影方法的标准步骤上，方便你阅读源码时带着地图走。
| 文件 | 在投影求解中的职责 |
|---|---|
main_repl_galizia.m | 复现 Galizia 模型的顶层入口，负责读参数、建网格、调用求解器、输出结果 |
main_housing_proj.m | 另一个带住房部门案例的顶层入口，用于对比不同模型设定 |
solve_proj.m | 投影法主循环，迭代更新多项式系数，直到欧拉方程残差满足要求 |
Smolyak_Polynomial.m | 计算 Smolyak 基多项式的函数，给定网格节点和阶数输出基底矩阵 |
Smolyak_Elem_Isotrop.m | 构造各向同性 Smolyak 插值权重和节点，是多维插值的基础 |
InvSubGen.m | 生成状态空间中的有效子网格，通常用于收缩或扩展边界 |
cubic_csd.m | 三维三次多项式样条相关的偏导与系数计算，用于较高阶的平滑近似 |
把main_repl_galizia.m当作入口，它的调用链大致是：读取.jnl参数，生成 Smolyak 节点，初始化系数，进入solve_proj.m迭代，迭代收敛后把决策函数写回工作区。Smolyak_Elem_Isotrop.m在这里的作用尤其关键，因为 Smolyak 插值不是直接张量积展开，而是从一组被稀疏化的一维插值算子中组合出高维多项式基底，这种方式可以有效规避“维度诅咒”。
2.3 参数化编程：改参数不用翻 Matlab 主文件
项目摘要里特别提到“参数化编程”，实际上指的就是把模型设定和求解配置从算法代码里剥离开。你只需在.jnl里修改参数，求解器会从配置结构体中读取。用一个典型的读取片段来说明这种组织方式：
% 读取模型配置，返回结构体 parameters parameters = parse_jnl('stnd_rbc_dyn.jnl'); disp(parameters.alpha_param); % 根据配置生成 Smolyak 优化节点 [grid, weights] = smolyak_elem_isotrop(parameters.max_degree, parameters.n_state);这里parse_jnl是我按 Promes 的惯例假设的读取函数；工具箱里可能叫别的名字，但模式相同。smolyak_elem_isotrop的第一个输入是最大多项式阶数，第二个输入是状态维度。注意max_degree不是越高越好：二维模型用阶数 3 通常足够，阶数 4 的基底数量增长很快，迭代耗时可能翻几倍。最稳妥的做法是从低阶开始，逐步提高阶数并观察残差是否下降，而不是一开始就贪高阶数。
3. 复现 Galizia 模型：main_repl_galizia.m 的参数修改与投影维度设定
明白了文件结构，接下来要真正跑通一个案例。main_repl_galizia.m是目录里最接近“教程”的脚本，它把 RBC 模型的资本存量和技术冲击作为两个状态变量，用 Smolyak 多项式逼近消费和劳动决策函数。你不需要理解每一行，只需要知道在哪里改参数、在哪里调精度。
3.1 打开入口脚本，先修改这些参数
编辑方式很简单，在 Matlab 命令窗口输入edit main_repl_galizia.m。脚本中通常有一段集中定义结构体params的代码，形如下面这样：
% 定义模型参数，注释中给出典型范围 params.alpha = 0.36; % 资本产出弹性，常用范围 0.3~0.4 params.beta = 0.99; % 贴现因子，季度模型通常在 0.98~0.995 params.gamma = 2.0; % 相对风险厌恶系数，常用 1.0~5.0 params.rhoA = 0.95; % 生产率冲击自回归系数 params.sigmaA = 0.01; % 冲击标准差，0.005~0.02 之间比较常见这里params.gamma对决策函数的非线性程度影响最大：当 gamma 从 1 提高到 5 时，消费的边际效用曲线变得更弯曲，投影方法需要更高阶多项式才能把决策函数逼近到相同精度。如果你发现残差在尾部变大，优先检查是不是这个参数设置得过于激进。
3.2 Smolyak 网格阶数：从 max_degree 开始的选择
main_repl_galizia.m中的maxDegree是投影核心参数，它直接控制基底数量。二维 Smolyak 格子在阶数较低时基底数增加还算温和，但每提升一阶，新增的组合项会让求解矩阵更稠密。常见做法是先设 2 阶跑通流程，再升到 3 阶看残差是否明显下降。
| maxDegree 取值 | 适用场景 | 迭代时间参考 |
|---|---|---|
| 1 | 快速验证代码是否跑通 | 秒级 |
| 2 | 论文初稿、参数扫描 | 分钟级 |
| 3 | 正式结果、残差要求 1e-6 以内 | 十几分钟到小时级 |
| 4 | 高非线性模型、状态变量少 | 数小时，需谨慎 |
上面的时间是针对二维模型在普通台式机上的表现，实际耗时还会受到迭代初值和容差影响。更重要的判断标准是残差：如果阶数从 2 升到 3 时最大残差下降了一个量级以上，说明模型确实需要更高阶的逼近；如果几乎没变化，说明问题不在阶数，而在状态空间边界或者校准参数。
3.3 运行日志 stnd_rbc_dyn.log 怎么看
运行结束后，程序会把迭代过程写入stnd_rbc_dyn.log。在 Matlab 里用type stnd_rbc_dyn.log就能查看，也可以直接用任何文本编辑器打开。一个健康的日志通常是这样的：
iter residual max_deriv time 0 1.2345e-01 2.3456e+00 0.3 1 3.4567e-03 4.5678e-01 0.4 2 1.2345e-05 5.6789e-03 0.4 3 2.3456e-08 6.7890e-05 0.4 converged at iteration 3如果残差迟迟降不下来，常见原因是初始系数给得太差，或者状态网格里有些点跑到了不稳定的参数区域。这时候不需要急着加多项式阶数，先检查日志中max_deriv列是不是在发散，若是发散，回到InvSubGen.m的边界收缩环节看状态空间是否超限。
3.4 对比案例 main_housing_proj.m
目录中另一个入口main_housing_proj.m是带住房存量的模型，它的状态维度比 RBC 多了一个或两个，因此 Smolyak 网格规模明显变大。运行它的方式完全相同，只需要在命令窗执行run main_housing_proj.m。对比两个入口脚本，你会发现主体结构相似，差别全在模型方程和参数个数上。这从侧面说明，Promes 的方法论是通用的：把模型方程写成欧拉残差函数，剩下的投影机制不需要大改。
4. solve_proj.m 内部机制：投影残差、不动点迭代与 InvSubGen 的网格生成
如果你只打算用工具箱跑现成模型，第 3 章已经够用；但要真正理解投影方法的边界，必须进入solve_proj.m内部。这一章讲清楚三件事：残差怎么构造、迭代如何更新系数、网格生成器在边界上做了什么手脚。
4.1 欧拉方程残差是如何变成一组代数方程的
投影方法不直接模拟随机路径，而是把决策函数设为有限项 Smolyak 多项式的组合：
% 决策函数近似形式 % x_t = smolyak_polynomial(s_t) * coefficient其中x_t是内生变量向量，s_t是状态变量向量，coefficient是待求系数。把这组决策函数代回模型的欧拉方程后，取条件期望，得到一个以状态变量为自变量的残差函数R(s; coefficient)。理论上任何状态点上都应该有R(s; coefficient)=0，但有限维多项式无法处处精确，于是只在网格节点上最小化残差。
cubic_csd.m在这里的典型用途是计算三次多项式拟合后的导数，这些导数会进入一阶条件和随机扰动项的 Jacobian，从而影响条件期望的近似质量。也就是说，它不是用来画图的，而是直接参与残差计算。
4.2 迭代更新：从初值到收敛的典型代码骨架
以下代码展示投影求解器常见的三种更新方式之一，即用 Matlab 的优化器直接求解非线性方程组：
% coefficient 是当前多项式系数 % residual_fn 返回网格节点上的残差向量 initial_coeff = zeros(n_basis, 1); options = optimoptions('fsolve', ... 'Display', 'iter', ... 'FunctionTolerance', params.iteration_tol, ... 'StepTolerance', 1e-10); coeff_final = fsolve(residual_fn, initial_coeff, options);这里的residual_fn每次调用都要做三件事：根据当前系数计算全部网格点的内生变量，使用 Gauss-Hermite 数值积分近似条件期望，最后减去静态方程给出的当期变量。迭代收敛标准主要有两个，一个是残差的无穷范数小于iteration_tol，另一个是系数步长足够小。实践中更推荐同时观察迭代路径，因为fsolve有时候会停留在局部极小值，日志中max_deriv指标就是用来辅助判断这个问题的。
4.3 InvSubGen.m：状态空间边界不是随便取的
InvSubGen的名字直译是“不变量子网格生成器”。在 DSGE 模型中，资本存量和技术冲击的组合并非所有实数对都有经济意义：资本不能为负，技术在某个区间外的概率趋近于零。如果直接在整个矩形区域布点，大量点落在经济含义之外，残差会被无意义的极端值主导，导致迭代发散。
这个函数的作用就是根据模型的一阶近似解先确定一个“能到达的区域”，然后在这个区域内生成密实分布的网格节点。实际实现中，常见做法是对二维状态变量做线性变换，把 [−1, 1]^d 的 Smolyak 节点映射到经济区间，而不是直接映射到原始矩形的上下界。
% 示意：把标准化节点映射到资本和技术冲击的有效区间 smolyak_node = smolyak_polynomial(grid, degree); % 基础节点 state_bound = inv_sub_gen(params); % 获得经济有效边界 state_node = interp1([-1 1], state_bound', smolyak_node);interp1在这里只是示意，实际代码会直接按张量积方式重排节点。但核心思想是：先求一个低阶近似的支持域，再在这个支持域上做全局逼近。这也是为什么InvSubGen.m必须和solve_proj.m配合：每完成一次迭代后，可以重新计算支持域，让网格向真实的不变量集收缩。
4.4 容差参数与迭代上限的取舍
solve_proj.m相关的容差参数通常从.jnl读取，也可在入口脚本中临时覆盖。一个比较合理的起点是设置iteration_tol = 1e-8，但注意，这个容差是对欧拉残差的范数而言，并不是对决策函数误差的范数。决策函数误差往往比残差小一两个量级，所以残差 1e-8 已经是非常严格的标准。论文复现中把容差放宽到 1e-6 和 1e-8 的结果差距通常远小于参数不确定性带来的差异，因此如果计算时间紧张，没必要硬顶 1e-10。
迭代上限和初值选择比容差更影响成败。对非线性较强的模型，第一次运行最好先用一阶解初始化系数，也就是把一阶线性近似的决策函数映射到 Smolyak 基上作为初值。这样做能让fsolve从一开始就靠近真解，比从零向量起步要稳健得多。代码中这种初值设置一般写成把一阶近似系数投影到多项式基底的操作，虽然只多十几行，却常常是收敛失败的分水岭。
5. 验证残差并调试非线性投影：三个能立刻上手的技巧
工具箱跑通只是第一步，真正用在论文或课程设计中，你需要证明投影解不是优化器凑出来的局部解。下面三个技巧是我在处理这类工具箱时一定会做的验证，顺序由浅入深。
5.1 在模拟路径上计算欧拉方程残差
网格节点上的残差小不等于模拟路径上残差小，因为模拟过程会不断到达网格边界附近，甚至短暂离开网格覆盖区域。正确的验证方式是先模拟一条长的状态路径，再把路径上每一点的欧拉方程残差算出来：
% 用投影解模拟 T 期，得到路径 sim_state 和 sim_control [sim_state, sim_control] = simulate_model(coeff_sol, params, T); % 计算路径上的最大欧拉残差 res_path = compute_path_residual(coeff_sol, sim_state, sim_control, params); fprintf('max residual on path: %.3e\n', max(abs(res_path)));如果路径残差比网格残差高两个量级以上，说明状态路径大量穿行在网格稀疏区域。解决方法不是粗暴提高阶数，而是用InvSubGen收缩网格边界，让节点密集分布在路径实际到达的区域。
5.2 用线性化解做差分测试
把投影决策函数在确定性稳态附近求数值导数，得到偏导数矩阵，然后把它和一阶近似方法给出乘子矩阵比较。如果两者差异很大，第一排查状态空间边界是否收缩过头，第二排查 Smolyak 基函数在边界附近的振荡是否导致导数失真。一个实用标准是：条件波动率在 5% 以内时，投影解的策略函数关键导数应与线性化解相差不超过 15%——这个阈值本身不是定理，却足够筛掉大多数实现错误。
5.3 调试负资本或虚数：先降波动，再降阶数
新用户最常遇到的坑是在投影求解过程中出现资本存量平方根项为负，或者对数取到 0 以下，导致残差函数返回 NaN。这不代表模型无解，通常只是状态区间的下界设置不合理。最快的处理套路是先把sigmaA降低一个量级，看看当前阶数下能否收敛；若收敛，再逐步放大波动的同时用前一次的解作为新初值。最后保留一个 3 阶 Smolyak 基做正式求解——大多数双状态 DSGE 模型足够，四阶以上只在研究非尾部大波动时才有必要。
本文还有配套的精品资源，点击获取
