使用Matlab和Simulink模拟Ramsey-Cass-Koopmans模型万博1manbetx

许多经济和金融模型，如资源分配或最优增长模型，都涉及没有显式解析解的微分方程组。对经济学家和其他金融分析师来说，用数字和可视化的方法解决这些系统是很重要的任务。

使用基本Ramsey-Cass-Koopmans（RCK）模型作为示例，本文介绍了用于创建，模拟和可视化常微分方程（ODES）系统的两个工作流程。一种方法是基于matlab^®，另一个在simulink上万博1manbetx^®．MATLAB方法使用在技术计算环境中工作的金融专业人员熟悉的编程技术。Simu万博1manbetxlink方法提供了一个可视化的建模环境和系统的图形表示。

万博1manbetxSimulink是一个用于建模带有反馈的时变系统的框图环境。这些系统在控制工程应用中是典型的，多年来影响了经济建模[1]。经济模型可以涉及具有许多方程和依赖关系的大规模ode系统。万博1manbetxSimulink提供了方便的特性，如子系统和模型引用，以支持大型系统的建模。万博1manbetx

本文中使用的代码和模型可供选择下载．

Ramsey-Cass-Koopmans模型

RCK模型旨在解释资本积累和消费增长的长期经济增长[2-4]。核心RCK模型是二维的，包括每个人均财富的两个耦合杂志（k)及人均消费量(c）。图1显示了系统的相位肖像。

人均财富的核心RCK模型方程(k)及人均消费量(c):

\[\压裂{dk} {dt} = f (k) - c (\ \ xi + \δφ+)k \四\压裂{直流}{dt} = c左(\ \压裂{f ' (k) - \θ-习\ \φδ}{\ρ}- \ \)\]

自从k和c出现在两个方程中，两个ode是耦合的。其中各项为:

• (f(k) = k^{\alpha}\)是一个生产函数，用\(k\)和一个资本弹性参数\(\alpha\)(产出生产对投入资本变化的响应)来衡量相对经济产出。
• (phi\)是劳动生产率的增长率(例如，由于技术创新或效率提高)。
• \（\ xi \）是劳动力供应的增长率（例如，由于迁移或人口增加）。
• (delta\)是资本的折旧率(例如，由于通货膨胀)。
• \（f'（k）= \ alpha k ^ {\ alpha - 1} \）为生产函数(f(k))的变化率(导数)。
• 是一个弹性参数，表示消费者在一段时间内倾向于平稳消费。
• 是消费者对其未来消费进行折现的比率(例如，表示对当前消费的偏好或试图保持其长期平均消费)。

使用MATLAB创建RCK模型

我们可以使用MATLAB功能直接解决许多杂散系统ODE45.，只要在时间区间\([t_0, t_f]\)上用标准形式表示系统\(frac{dY}{dt} = G(t,Y)\)，且满足初始条件\(Y(0) = Y_0\)。注意，当方程组中存在多个未知时间函数时，\(Y\)和\(G(t,Y)\)是向量值。

我们首先定义结构变量中的模型参数参数个数，然后写一个向量值函数RCK_Equations表示标准微分方程的右边\(\frac{dY}{dt} = G(t,Y)\)。这个函数返回一个包含两个元素的向量，每一步都包含\(\frac{dk}{dt}\)和\(\frac{dc}{dt}\)。

我们还写了两个辅助函数RCK_f和rck_df.分别返回生产函数\(f(k)\)及其导数\(f'(k)\)的值。我们将每个函数保存在单独的文件中，这样就可以很容易地研究不同的生产函数对数值结果的影响。

定义定义Ramsey-Cass- %Koopmans模型的两个%耦合常微分方程右边的函数。%提取k和c k = Y(1);c = Y (2);写出dk/dt和dc/dt的方程。RCK_f(k, params) - c -…(参数。φ+参数。* k;% dk/dt dY_dt(2,1) = ((RCK_df(k, params) - params. log . log . log . log . log . log . log . log . log . log。θ——参数。xi -…) / params.delta。ρ-参数。* c; % dc/dt end % RCK_Equations

下一步是创建一个函数句柄(＠的输入函数ODE45.通过参数化RCK_Equations使用预定义的参数结构参数个数．这个函数必须是时间的函数(t)及州(Y只)。

RCK_Fun = @(t, Y) RCK_Equations(t, Y, params);

我们通过修改求解器选项来确保人均财富和消费随时间的推移保持非负odeset．

选择= odeset（'nonnegative'，[1,2]）;

取初始条件为\(k_0 = 25\)和\(c_0 = 2\)并创建一个从0到1.5单位的时间向量，我们现在可以用数值方法求解方程组ODE45.．输出ODE45.是时间和状态。注意，因为我们正在创建一个时间向量t为ODE45.直接，只需要返回第二个输出参数Y从ODE45.．

y0 = [25;2];t = linspace（0,1.5,5000）;[~， Y] = ode45(RCK_Fun, t, Y0, opts);k_out = y（：，1）;％产出人均财富C_OUT = Y（:,2）;人均消费％产量

我们使用MATLAB的可视化函数彗星创建解决方案路径的动画轨迹，显示资本和消费的时间演变。这个动画的最后一帧叠加在相平面上，如图2所示。红线是图1中曲线\(\frac{dk}{dt} = 0\)的一小部分。

用时间消元法求解系统的稳态

通过求解方程组\(\frac{dk}{dt} = \frac{dc}{dt} = 0\)，可以得到系统的稳态(s)。解\(frac{dk}{dt} = 0\)将得到曲线\(c^{*}(k)=f(k)-(\phi + \xi + \delta)k\(图1中的红色曲线)。，即。

\ [\ hat {k} = \ left（\ frac {\ alpha} {\ rho \ phi + \ theta + \ xi + \ delta} \ revaly）^ {\ frac {1} {1 - \ alpha}} \]

这个值定义了图1中的垂直线\(k = \hat{k}\)。我们可以使用Meshgrid.函数创建晶格K和C资本/消费网格点。计算差异后dK和直流如RCK方程所定义，我们可以使用streamslice可视化函数以渲染\（（k，c）\）平面中的Streamlines。

流切片(K, C, dK, dC)

叠加曲线\(c^{*}\)和\(k = \hat{k}\)创建如图1所示的相位画像。

如Carrol所示[3]，模型过渡到稳态没有解析解。但是，我们可以使用时间消去技术得到:

\ [\ frac {dc} {dk} = \ frac {dc / dt} {dk / dt} = \ frac {c（f'（k） - \ theta - \ xi - \ delta - \ phi \ rho）}{rho（f（k） - c - c - （\ phi + \ xi + \ delta）k）}}} \]

对k给出了\(c)作为\(k)的函数的解轨迹。为了避免在\(frac{dc}{dk})的分子或分母为零时计算\(frac{dc}{dk}\)的问题，我们将\(k\)域分成两部分，一部分在\(hat{k}\)的左边，另一部分在[2]的右边。应用ODE45.给出了解决方案轨迹(图3)。

请注意，上层解决方案路径平滑，而较低的解决方案路径在\（\ FRAC {DK} {dt}×0 \）附近的数值不稳定性，这是某些刚性系统的特征。在我们的示例中，而不是使用ODE45.，我们将使用一个求解器来处理刚性系统，例如ode15s．为了提高求解轨迹的可靠性，我们计算了系统的雅可比矩阵，并将其传递给求解器odeset．得到的平滑路径如图4所示。对于更大或更复杂的系统，我们可以使用Symbolic Math Toolbox™来计算解析雅可比矩阵，而无需手工计算。

使用Simulink创建RCK模型万博1manbetx

我们使用财务分析师熟悉的技术在MATLAB中解决了RCK系统。现在我们将采用一种不太熟悉但更直观的方法:我们将在Simulink中使用预定义块的库动态地实现系统。万博1manbetx

在Si万博1manbetxmulink中，数据表示为信号和参数。信号是连接块的线条，并且表示时变数据，例如衍生物\（\ frac {dk} {dt} \）。参数是存储在块内的系统值 - 例如，初始条件\（k_0 \）。

在Simulink中对ode建模时，我们从C万博1manbetxontinuous库中的Integrator块开始。这个块对它的输入信号(导数)进行数值积分。因为系统有一阶导数k和c，我们使用两个集成商块。初始条件\（k_0 = 5 \）和\（c_0 = 3 \）被分配为块内的参数（图5）。

我们使用来自Simulink库的块来实现RCK方程的右侧(表1)。万博1manbetx

表1。万博1manbetxSimulink块用来实现右边的RCK方程。

随着我们的模型的大小和复杂性增加，我们可以通过使用子系统块将块分组为子系统来简化它。我们在子系统中封装了生产函数\（f（k）\）（图6）。

我们使用导数块来近似导数\(f'(k)\)。注意，Derivative块没有初始条件参数，因此不应该用作建模ode的起点。

图7显示了完整的RCK模型。

在建立模型后，我们指定仿真停止时间为500个时间单位，并开始仿真。图8显示了从点\((k_0,c_0) =(5,3)\)开始的最终轨迹。

并行运行模拟

作为模型分析的一部分，我们可能希望通过使用不同的参数值集运行模拟来研究模型对其参数的依赖关系。每个仿真可以独立运行，并使用parfor从并行计算工具箱™中构造。

从基于MATLAB的模型实现开始，我们创建网格点的格K0和C0表示我们要研究的不同初始条件。的每个迭代parfor循环，我们选择不同的初始条件Y0并存储输出k_out.和c_out使用细胞阵列。

RCK_Fun = @(t, Y) RCK_Equations(t, Y, params);opts = odeset('非负'，[1,2]);T = linspace(0, 1.5, 1000);parfor k = 1:numel(K0) %人均财富和消费的初值。Y0 = [K0 (k);C0 (k)];解决耦合系统。[~， Y] = ode45(RCK_Fun, t, Y0, opts);k_out{k} = Y(:， 1);%人均财富产出c_out{k} = Y(:， 2); % Output per-capita consumption end % parfor

使用100乘100的初始条件格意味着我们要对模型进行10,000次并行模拟。这些模拟生成了如图9所示的解决轨迹。

为了并行地模拟Simulink 万博1manbetxRCK模型，我们做了以下工作:

1. 在每个worker上加载模型load_system.和spmd构造。
2. 定义数组K0和C0对于我们想要模拟模型的初始条件。
3. 方法编写函数以编程方式模拟模型sim卡命令。(注意，要在模型中使用set_param，我们需要将数值转换为文本。这try / catch为孤立的初始条件集构造防止意外收敛问题的防护措施。)
4. 的每个迭代parfor循环时，使用不同的初始条件对调用函数。

每个worker加载一次模型。spmd load_system(“RCK_Model”);end % spmd %%并行执行模拟。parfor k = 1:numel(K0) simout(k) = runSim(K0(k)， C0(k));函数模拟不同初始值的模型，人均财富和消费使用刚性系统解算器% ode15s和45个时间单位的停止时间。%将人均财富和消费的初始值格式化为文本。k0 = num2str (k0);c0 = num2str (c0);set_param('RCK_Model/capital'， 'InitialCondition'， k0) %运行模拟。try simout = sim('RCK_Model'， 'Solver'， 'ode15s'， 'StopTime'， '45'); catch % If a simulation run fails to converge, assign an empty output. simout = Simulink.SimulationOutput; end % try end % runSim

图10显示了10,000个模型并行仿真的解决轨迹。

概括

在本文中，我们演示了一个耦合ode系统可以使用MATLAB或Simulink进行模拟，并且每种方法都有好处。万博1manbetx

利用MATLAB，将方程转化为所需的形式ODE45.并使用功能处理。在Simulink中使用Integrator块，我们以其本地万博1manbetx形式实现了差分方程。

MATLAB和Simulink使万博1manbetx我们能够自动计算导数。在MATLAB中，我们使用梯度在符号数学工具箱中的功能。在Si万博1manbetxmulink中，我们使用衍生块和求解器Jacobian配置设置ode15s．

万博1manbetxSimulink简化了大型或复杂的ode系统的建模。子系统通过对功能相关的块进行分组来帮助我们组织模型。万博1manbetx

这两种方法都支持并行化。在MATLAB和Simulink中，我万博1manbetx们使用parfor构造。在Si万博1manbetxmulink中，我们使用sim卡命令以编程方式模拟模型。

MATLAB和Simulink都万博1manbetx为ode的求解和可视化提供了一个集成的建模环境。使用的方法取决于ode系统的大小和复杂性。使用Simulink布局模型是快速、可视化和直观的。万博1manbetx如果您是Simulink的新手，可以万博1manbetx在mathworks.com上找到大量的资源来帮助您入门。