使用Matlab和Simulink模拟Ramsey-Cass-Koopmans模型金宝app

由索尼娅桥和肯德埃利,Mathworks


许多经济和金融模式,例如资源分配或最佳增长,涉及差分方程式,没有明确的分析解决方案。在数字上解决这些系统和可视化结果是经济学家和其他金融分析师的重要任务。

使用基本Ramsey-Cass-Koopmans(RCK)模型作为示例,本文介绍了用于创建,模拟和可视化常微分方程(ODES)系统的两个工作流程。一种方法是基于matlab®,另一个在simulink上金宝app®。Matlab方法使用在技术计算环境中工作的金融专业人士熟悉的编程技术。Simu金宝applink方法提供了一种可视化建模环境和系统的图形表示。

金宝appSimulink是一个框图环境,用于建模具有反馈的时变系统。这种系统在控制工程应用中是典型的,这多年来影响了经济建模[1]。经济模型可以涉及具有许多方程和依赖性的大规模杂散系统。金宝appSimulink提供方便的功能,如子系统和模型引用,以支持大型系统的建模。金宝app

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

Ramsey-Cass-Koopmans模型

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

图1:常微分方程Ramsey-Cass-Koopmans系统的相位画像。

人均财富的核心RCK模型方程(K.)和人均消费(C) 是:

\ [\ frac {dk} {dt} = f(k) - c - c - (\ phi + \ xi + \ delta)k,\ quad \ frac {dc} {dt} = c \ left(\ frac {f'(k) - \ theta - \ xi - \ delta} {\ rho} - \ phi \右)\]

自从K.C在两个方程中,两个ode是耦合的。这些方程中的项如下:

  • \(f(k) = k^{\alpha}\)是一个以\(k)衡量相对经济产出的生产函数和一个资本弹性参数\(\alpha\)(产出生产对投入资本变化的反应)。
  • (\phi\)是劳动生产率的增长率(例如,由于技术创新或效率的提高)。
  • \(\ xi \)是劳动力供应的增长率(例如,由于迁移或人口增加)。
  • (\delta)是资本的折旧率(例如,由于通货膨胀)。
  • \(f'(k)= \ alpha k ^ {\ alpha - 1} \)是生产函数\(f(k)\)的变化率(导数)。
  • \(\ theta \)是一种弹性参数,表明消费者随着时间的推移平滑消费量的趋势。
  • \(\ rho \)是消费者折扣其未来消费的速度(例如,通过表示偏好立即消费或试图维持其长期平均消费)。

使用MATLAB创建RCK模型

我们可以直接使用MATLAB功能来解决许多ODES系统ODE45.如果我们在时间间隔\([t_0,t_f] \)上以标准表单\(\ frac {dy} {dt} {dt} = g(t,y)\)表达系统,并且受初始条件\(y(0)= y_0 \)。请注意,如果在方程式系统中存在多个未知的时间功能,则向\(y \)和\(g(t,y)\)是矢量值。

我们首先定义结构变量中的模型参数参数个数,然后编写矢量值函数RCK_Equations表示标准差分方程的右侧\(\ frac {dy} {dt} = g(t,y)\)。此函数返回一个两个元素向量,其中每个时间步在每个时间步骤都包含\(\ frac {dk} {dt} {dt}×{dt} {dt})的值。

我们还编写了两个辅助功能RCK_fRCK_DF.返回生产函数\(f(k)\)的值及其导数\(f'(k)\)。我们将每个函数保存在单独的文件中,以便很容易调查不同生产功能对数字结果的影响。

函数dy_dt = rck_equations(t,y,params)%rck_equations函数定义两个%耦合的常微分方程的右手侧,定义Ramsey-Cass%Koopmans模型。%提取物K和C.k = y(1);c = y(2);%写下DK / DT和DC / DT的方程。dy_dt(1,1)= rck_f(k,params) -  c  -  ...(params.phi + params.xi + params.delta)* k;%dk / dt dy_dt(2,1)=((rck_df(k,params) -  params.theta  -  params.xi  -  ... params.delta)/ params.rho  -  params.phi)* c;%DC / DT末端%RCK_EQUATION

下一步是创建功能句柄(@)包含输入功能ODE45.通过参数化RCK_Equations使用预定义的参数结构参数个数。这个函数必须是时间(T.)和国家(y) 只要。

RCK_Fun = @(t, Y)

我们通过修改求解选项,确保每人的财富和消费都随着时间的推移而仍然是非线性的odeset.

选择= odeset('nonnegative',[1,2]);

将初始条件置于\(k_0 = 25 \)和\(c_0 = 2 \)并创建从0到1.5单位的时间向量,我们现在可以用数字地解决系统ODE45.。输出ODE45.是时间和国家。注意,因为我们正在创建时间向量T.为了ODE45.直接,只需要返回第二个输出参数yODE45.

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)。

图2.从点\((k_0,c_0)=(25,2)=(25,2)\)开始的解决方案轨迹通过直接使用的等式的耦合系统获得ODE45.

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

通过求解方程可以找到系统的稳态\(\ FRAC {DK} {dt} = \ FRAC {DC} {dt} = 0 \)。解决\(\ frac {dk} {dt} = 0 \)产生曲线\(c ^ {*}(k)= f(k) - (\ phi + \ xi + \ delta)k \)(红色图1中的曲线)。解决\(\ frac {dc} {dt} = 0 \)产生单个值为\(k \)即。

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

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

Streamslice(K,C,DK,DC)

覆盖曲线\(c ^ {*} \)和\(k = \ hat {k})创建图1所示的相位画象。

正如carol所示[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)}} \]

integrK.为\(c \)提供一个解决方案轨迹作为\(k)的函数。避免问题评估\(\ frac {dc} {dk})当其分子或分母为零时,我们将\(k \) - 域分为两个部分,左侧为\(\ hat {k}\)和一个到右[2]。申请ODE45.给出解决方案轨迹(图3)。

图3.使用时间消除方法获得的消耗策略的上下解决方案路径。

请注意,上解决方案路径平滑,而较低的解决方案路径在\(\ FRAC {dk} {dt} {dt} {dt}×0 \)附近的数值不稳定性,而某些刚性系统的特征则受到数值不稳定性。在我们的示例中,而不是使用ODE45.,我们将使用一个设计来处理刚性系统的解算器,例如ode15s.。为了提高解决方案轨迹的可靠性,我们计算系统的雅蟒并通过求助者odeset.。所得到的平滑路径如图4所示。对于更大或更复杂的系统,我们可以使用Symbolic Math Toolbox™来计算解析雅可比矩阵,而无需手动计算。

图4.使用僵硬求解器获得的较低解决方案路径ode15s.

使用Simulink创建RCK模型金宝app

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

在Si金宝appmulink中,数据表示为信号和参数。信号是连接块的线条,并且表示时变数据,例如衍生物\(\ frac {dk} {dt} \)。参数是存储在块内的系统值 - 例如,初始条件\(k_0 \)。

在Simulink中建模ode时,我们从Co金宝appntinuous库中的Integrator模块开始。该模块对其输入信号(导数)进行数值积分。因为系统有一阶导数K.C,我们使用两个集成商块。初始条件\(k_0 = 5 \)和\(c_0 = 3 \)被分配为块内的参数(图5)。

图5. K和C的集成器块。红线表示尚未连接到其他块的信号。

我们使用Simulink库中的块来实现RCK方程的右侧(表1)。金宝app

堵塞 象征 目的
不变 参考模型参数(例如,params.phi
产品 乘法信号
添加或减去信号
获得 将一个信号乘以或除以一个常数
数学函数 数学运算(例如,权力和对数)
外港 将结果传递到MATLAB工作空间(例如,K.C
衍生物 数值近似衍生物
子系统 对功能相关的块进行分组

表1。金宝appSimulink块用于实现RCK方程的右侧。

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

图6.作为子系统实现的生产函数\(f(k)= k ^ {\ alpha})。

我们使用衍生块近似衍生\(f'(k)\)。请注意,衍生块没有初始条件参数,因此不应用作建模ODES的起点。

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

图7.完整的Simulink RCK模型。金宝app

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

图8.从点开始的解决方案轨迹((k_0,c_0)=(5,3)\)。

并行运行模拟

作为模型分析的一部分,我们可能希望使用不同的参数值运行模拟来研究模型对其参数的依赖性。每个模拟都可以独立运行并使用该模拟并行实现议案构造从并行计算工具箱™。

从基于MATLAB的模型实现开始,我们创建网格点的格子K0C0.代表我们想要调查的不同初始条件。在每次迭代中议案循环,我们选择不同的初始条件Y0并存储输出k_out.c_out使用细胞阵列。

RCK_Fun = @(t, Y)opts = odeset('NonNegative', [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所示的解决方案轨迹。

图9。不同初始条件下RCK模型的求解路径。

要在并行模拟Simulink R金宝appCK模型,我们执行以下操作:

  1. 使用每个工人加载模型load_system.SPMD.构造。
  2. 定义阵列K0C0.对于我们想要模拟模型的初始条件。
  3. 控件编写函数以编程方式模拟模型sim卡命令。(注意,以编程方式在模型中设置参数使用set_param,我们需要将数值转换为文本。这试着抓构造针对孤立初始条件集的意外收敛问题的防护措施。)
  4. 在每次迭代中议案循环,使用不同的初始条件调用函数。
%%每个工作人员加载一次模型。spmd load_system('rck_model');终端%SPMD %%并行执行模拟。k = 1:numel(k0)simout(k)= runsim(k0(k),c0(k));终端%Parcon FILE SIMOUT = runsim(k0,c0)%运行函数模拟模型,用于使用刚性系统求解器%ode15的每个人均财富和消费的不同初始值%和45个时间单位的停止时间。%格式化人均财富的初始值和作为文本的消费。k0 = num2str(k0);c0 = num2str(c0);set_param('rck_model / capital','initialcondition',k0)set_param('rck_model / cosption','initialcondition',c0)%运行模拟。尝试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个平行模拟的解决方案轨迹。

图10.从不同初始条件开始的RCK模型的解决方案路径。

概括

在本文中,我们演示了可以使用MATLAB或Simulink对耦合ode系统进行模拟,并且每种方法都有好处。金宝app

使用MATLAB,我们将方程转换为所需的形式ODE45.并使用功能处理。在Simulink中使用Integrator块,我们以其本地金宝app形式实现了差分方程。

Matlab和Simulink都金宝app让我们自动计算衍生产品。在Matlab中,我们使用坡度在符号数学工具箱中的功能。在Si金宝appmulink中,我们使用衍生块和求解器Jacobian配置设置ode15s.

金宝appSimulink简化了大型或复杂的ODES系统的建模。子系统帮助我们通过分组功能相关的块来组织模型,Simulink模型窗口提供了模型的直观视觉布局。金宝app

两种方法都能平行化。在Matlab和Simulink中,我金宝app们使用议案构造。在Si金宝appmulink中,我们使用sim卡命令以编程方式模拟模型。

MATLAB和SIMULINK都金宝app提供了一个集成的建模环境,用于解决和可视化ODES系统。您使用的方法取决于您的ODES系统的大小和复杂性。使用Simulink铺设模型是快速,视觉和直观的。金宝app如果您是Simulink的新功能,M金宝appathWorks.com上有充足的资源可以帮助您开始。

发布2016年 - 93062V00

查看相关功能的文章

查看相关行业的文章