马尔可夫调制的随机时滞IS-LM模型的随机稳定性仿真

Simulation of Stochastic Stability of Delayed IS-LM Model with Markovian Switching

摘要: 本文研究了标准随机IS-LM模型处于税收过程中存在时滞及二元离散状态值的马尔可夫链随机性波动环境中，在给定参数条件下通过利用随机稳定性方法探讨了马尔可夫调制的标准随机时滞IS-LM模型的指数稳定性行为，得到了解的仿真结果。

Abstract: This paper studies that the economic system is in the random environment which has delayed time in taxation and Markov chain with two discrete values. By using the method of stochastic stability and the given parameters, this paper analyzes the stability and simulation of the solution of the standard random IS-LM model with Markovian Switching and delayed time.

文章引用：周路军. 马尔可夫调制的随机时滞IS-LM模型的随机稳定性仿真[J]. 统计学与应用, 2018, 7(1): 1-6. https://doi.org/10.12677/SA.2018.71001

1. 引言

随机微分方程一般常用来描述经济系统的不确定性因素的影响，而如自然灾害的突然发生，或者经济政策的改变，或者企业生产技术的革新等会带来经济系统结构性的变化，因此具有离散状态特征的马氏链也会引入到随机微分方程中。

所以，具有马尔可夫链调制的随机微分方程与随机时滞微分方程的研究大量应用于经济、金融、控制等领域，其相应的稳定性理论、控制理论与数值仿真得到了极大的发展。

在标准的IS-LM经济周期模型应用方面，马氏链调制的随机IS-LM经济周期模型的研究尚少，尤其是模型中考虑到税收过程中存在时滞的情形，研究更少。本文考虑两状态的马尔可夫链到随机时滞IS-LM模型中，研究它的随机稳定性并给出数值仿真例子。

2. 预备知识及模型构建

Mao (2001)研究了如下的马尔可夫调制的随机微分方程 [1] ：

$\text{d}x\left(t\right)=f\left(x\left(t\right),x\left(t-\tau \right),t,r\left(t\right)\right)\text{d}t+g\left(x\left(t\right),x(t-\tau ),t,r\left(t\right)\right)\text{d}w(\; t\; )$

这里 $t\ge 0$ ，具有初始值 $x_0\in C_{F_0}^b\left(\Omega ;R^n\right),f:R^n\times R_{+}\times S\to R^n,g:R^n\times R_{+}\times S\to R^{n\times m}$ 。

假如A是一个向量，则用 $|\cdot |$ 表示Rⁿ中欧氏范数。假如A是一个矩阵(用A^T表示A的转置)，则定义 $|A|=\sqrt{trace\left(A^{\text{T}}A\right)}$ 为迹范数，定义 $\Vert A\Vert =\sup \left\{|Ax|:|x|=1\right\}$ 为算子范数。设 $r\left(t\right),t\ge 0,$ 是定义在完备概率空间 $\left(\Omega ,F,{\left\{F_t\right\}}_{t\ge 0},P\right)$ 上的右连续马尔可夫链，取两个值的状态空间 $S=\left\{1,2\right\}$ ，且具有生成元 $\Gamma ={\left({\gamma }_{ij}\right)}_{2\times 2}$ ，其中

$P\left\{r\left(t+\Delta \right)=j|r\left(t\right)=i\right\}=\{\begin{array}{l}{\gamma }_{ij}\Delta +o\left(\Delta \right),\\ 1+{\gamma }_{ij}\Delta +o\left(\Delta \right),\end{array}$ $\begin{array}{l}i\ne j\\ i=j\end{array}$

这里 $\Delta >0$ ， ${\gamma }_{ij}\ge 0$ 是由i状态转移到j状态的转移速率，且 ${\gamma }_{ii}=-{\displaystyle \underset{j}{\sum }{\gamma }_{ij}},i\ne j$ 。

马尔可夫链 $r\left(t\right)$ 与布朗运动 $w(•)$ 假设是相互独立的，则 $r\left(t\right)$ 的每一个样本轨迹都是一个右连续的阶梯函数，且在每一个有限子区间上只有有限个跳跃点。

V. Torre (1977)用分岔理论分析了非均衡市场的完全凯恩斯系统 [2]

$\{\begin{array}{l}\stackrel{\dot{}}{Y}\left(t\right)=\alpha \left(I\left(Y\left(t\right),R\left(t\right)\right)-S\left(Y\left(t\right),R\left(t\right)\right)\right),\\ \stackrel{\dot{}}{R}\left(t\right)={\beta }_0\left(L\left(Y\left(t\right),R\left(t\right)\right)-L_s\right),\end{array}$

其中 $\frac{\partial I}{\partial Y}>0,\frac{\partial I}{\partial R}<0,\frac{\partial S}{\partial Y}>0,\frac{\partial S}{\partial R}<0,\frac{\partial L}{\partial Y}>0,\frac{\partial L}{\partial R}<0$ ；Y，R，I(Y,R)，S(Y,R)，L(Y,R)，Ls， $\alpha >0$ 和 $\beta >0$ 分别表示产品供给函数、利率、投资(产出和利率的函数)、储蓄(产出和利率的函数)、货币流动需求量、常数货币流动供给量、产品市场和货币市场调整参数。第一个方程表示产品市场的非均衡调整，第二个表示货币市场的非均衡调整。

本文中的第一个方程利用萨谬尔森的动态乘数-加速数模型 [3] 的第一个方程即 $\stackrel{\dot{}}{Y}\left(t\right)=\alpha \left(D\left(t\right)-Y\left(t\right)\right)$ ，其中需求函数 $D\left(t\right)=C\left(t\right)+I\left(t\right)+G$ 表示成消费函数、投资函数及政府支出的和。而消费函数 $C\left(t\right)=\sigma Y^d\left(t\right)=\sigma \left(Y\left(t\right)-T\left(t\right)\right)$ 是关于可支配收入的函数。税收函数是产品供给函数的正常数倍，且税收过程中存在时间延迟现象，因此引入时滞 $\tau$ ，税收函数为 $T\left(t\right)=\epsilon Y\left(t-\tau \right)$ ；投资函数是关于供给和利率的函数 $I\left(Y,R\right)=\gamma Y\left(t\right)-{\beta }_{1}R\left(t\right)$ ；G是常数政府支出。

本文第二个方程中假设流动性需求函数 $L\left(Y\left(t\right),R\left(t\right)\right)=uY\left(t\right)-{\beta }_{2}R\left(t\right)$ ，供给函数 $L_s$ 是常数。则考虑的时滞IS-LM模型为

$\{\begin{array}{l}\stackrel{\dot{}}{Y}\left(t\right)=\alpha \left(\left(\sigma +\gamma -1\right)Y\left(t\right)-\sigma \epsilon Y\left(t-\tau \right)-{\beta }_{1}R\left(t\right)+G\right),\\ \stackrel{\dot{}}{R}\left(t\right)=\beta {}_0\left(uY\left(t\right)-{\beta }_{2}R\left(t\right)-L_s\right)\end{array}$ (1)

这里参数 $0<\alpha ,\sigma ,\gamma ,\epsilon ,{\beta }_0,{\beta }_1,{\beta }_2<1$ 。

命题1：在条件 $u{\beta }_1-{\beta }_2\left(\sigma \left(1-\epsilon \right)-1+\gamma \right)>0$ 和 $uG+\left(\sigma \left(1-\epsilon \right)-1+\gamma \right)L_s>0$ 都满足的情况下，标准时滞IS-LM模型(1)存在唯一的正均衡解

$\left(Y^{*},R^{*}\right)=\left(\frac{{\beta }_{2}G+{\beta }_1{L}_s}{u{\beta }_1-{\beta }_2\left(\sigma \left(1-\epsilon \right)-1+\gamma \right)},\frac{uG+\left(\sigma \left(1-\epsilon \right)-1+\gamma \right)L_s}{u{\beta }_1-{\beta }_2\left(\sigma \left(1-\epsilon \right)-1+\gamma \right)}\right)$ 。

令 $x\left(t\right)=Y\left(t\right)-Y^{*},x\left(t-\tau \right)=Y\left(t-\tau \right)-Y^{*}$ ， $y\left(t\right)=R\left(t\right)-R^{*}$ ，则模型(1)变成

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=\alpha \left(\left(\sigma +\gamma -1\right)x\left(t\right)-\sigma \epsilon x\left(t-\tau \right)-{\beta }_{1}y\left(t\right)\right),\\ \stackrel{\dot{}}{y}\left(t\right)={\beta }_0\left(ux\left(t\right)-{\beta }_{2}y\left(t\right)\right)\end{array}$ (2)

由于此系统受到随机扰动的影响，令变量 $x\left(t\right)$ ， $x\left(t-\tau \right)$ 和 $y\left(t\right)$ 对应的参数在两个方程中分别受到强度为 ${\sigma }_{1}w\left(t\right)$ ， ${\sigma }_{2}w\left(t\right)$ ， ${\sigma }_{3}w\left(t\right)$ ， ${\sigma }_{4}w\left(t\right)$ 和 ${\sigma }_{5}w\left(t\right)$ 的随机扰动。

即模型(2)变为

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=\left(\alpha \left(\sigma +\gamma -1\right)+{\sigma }_1\omega \left(t\right)\right)x\left(t\right)+\left(-\alpha \sigma \epsilon +{\sigma }_2\omega \left(t\right)\right)x\left(t-\tau \right)+\left(-\alpha {\beta }_1+{\sigma }_{3}w\left(t\right)\right)y\left(t\right),\\ \stackrel{\dot{}}{y}\left(t\right)=\left(\left({\beta }_{0}u+{\sigma }_4\omega \left(t\right)\right)x\left(t\right)+\left(-{\beta }_0{\beta }_2+{\sigma }_{5}w\left(t\right)\right)y\left(t\right)\right)\end{array}$ (3)

用矩阵表示为

$\text{d}X\left(t\right)=\left(AX\left(t\right)+BX\left(t-\tau \right)\right)dt+\left(CX\left(t\right)+DX\left(t-\tau \right)\right)\text{d}w\left(t\right)$ (4)

这里 $X\left(t\right)=\left(\begin{array}{l}x\left(t\right)\\ y\left(t\right)\end{array}\right),X\left(t-\tau \right)=\left(\begin{array}{l}x\left(t-\tau \right)\\ y\left(t-\tau \right)\end{array}\right),$

$A=\left(\begin{array}{cc}\alpha \left(\sigma +\gamma -1\right)& -\alpha {\beta }_1\\ {\beta }_{0}u& -{\beta }_0{\beta }_2\end{array}\right),B=\left(\begin{array}{cc}-\alpha \sigma \epsilon & 0\\ 0& 0\end{array}\right)$ ， $C=\left(\begin{array}{cc}{\sigma }_1& {\sigma }_3\\ {\sigma }_4& {\sigma }_5\end{array}\right),D=\left(\begin{array}{cc}{\sigma }_2& 0\\ 0& 0\end{array}\right)$ 。

考虑到上涨与紧缩周期性波动的影响，将二值离散状态的马尔可夫链 $r\left(t\right)$ 引入到模型(4)，则模型变成

$dX\left(t\right)=\left(A\left(r\left(t\right)\right)X\left(t\right)+B\left(r\left(t\right)\right)X\left(t-\tau \right)\right)\text{d}t+\left(C\left(r\left(t\right)\right)X\left(t\right)+D\left(r\left(t\right)\right)X\left(t-\tau \right)\right)\text{d}w\left(t\right)$ (5)

这里r(t) = 1，2。生成元 $\Gamma =\left(\begin{array}{cc}{\gamma }_{11}& {\gamma }_{12}\\ {\gamma }_{21}& {\gamma }_{22}\end{array}\right),{\gamma }_{11}=-{\gamma }_{12}<0,{\gamma }_{22}=-{\gamma }_{21}<0.$

3. 案例仿真

例令

${\alpha }_1=0.011,{\sigma }_1=0.018,{\gamma }_1=0.049,{\beta }_{01}=0.024,{\beta }_{11}=0.056,{\beta }_{21}=0.035,u_1=0.046$ ，

${\epsilon }_1=0.021,{\sigma }_{11}=0.048,{\sigma }_{21}=0.024,{\sigma }_{31}=0.015,{\sigma }_{41}=0.071,{\sigma }_{51}=0.171$ ，

${\alpha }_2=0.011,{\sigma }_2=0.016,{\gamma }_2=0.029,{\beta }_{02}=0.025,{\beta }_{12}=0.065,{\beta }_{22}=0.065,u_2=0.026$ ，

${\epsilon }_2=0.031,{\sigma }_{12}=0.058,{\sigma }_{22}=0.016,{\sigma }_{32}=0.025,{\sigma }_{42}=0.061,{\sigma }_{52}=0.271$

及生成元 $\Gamma =\left(\begin{array}{cc}-3& 3\\ 1& -1\end{array}\right)$ 。

命题1条件 $\{\begin{array}{l}{u}_1{\beta }_{11}-{\beta }_{21}\left({\sigma }_1\left(1-{\epsilon }_1\right)-1+{\gamma }_1\right)>0,\\ u_2{\beta }_{12}-{\beta }_{22}\left({\sigma }_2\left(1-{\epsilon }_{21}\right)-1+{\gamma }_2\right)>0\end{array}$ 满足，则模型(5)存在唯一平凡均衡解。

通过Euler-Maruyama方法 [4] 对案例进行MATLAB仿真，分别在两个不同随机样本和不同时滞( $\tau =0.2,0.8$ )下，得到图1，图2，图3，图4。

图1，图2描述的是在时滞为0.2的情况下两种不同随机样本下，模型(5)的解在时间t < 50以内波动性非常大。但随着时间t的增加(t > 150)模型(5)的解在不断地趋近于平凡均衡解。总体上呈现随机稳定。

图3，图4描述的是模型(5)在时滞为0.8的两种不同随机样本情况。图3中，模型(5)的解在时间t < 50以内波动性非常大，但随着时间t的增加(t > 150)在不断地趋近于均衡解。而图4中，模型(5)的解在时间t < 150以内波动性都非常大，而在t > 150后波动性较少且不断地趋近于均衡解。总体上都呈现随机稳定。

Figure 1. Solution x(t), y(t) of sde and r(t) (time lag = 0.2)

图1. 随机微分方程解x(t)，y(t)及马氏链r(t) (时滞 = 0.2)

Figure 2. Solution x(t), y(t) of sde and r(t) (time lag = 0.2)

图2. 随机微分方程解x(t)，y(t)及马氏链r(t) (时滞 = 0.2)

Figure 3. Solution x(t), y(t) of sde and r(t) (time lag = 0.8)

图3. 随机微分方程解x(t)，y(t)及马氏链r(t) (时滞 = 0.8)

Figure 4. Solution x(t), y(t) of sde and r(t) (time lag = 0.8)

图4. 随机微分方程解x(t)，y(t)及马氏链r(t) (时滞 = 0.8)

从图1~图4看出，r(t)是一个连续的阶梯函数，且在固定的子区间里是有限跳跃的。

基金项目

教育部人文社会科学研究一般项目(11YJCZH259)。

参考文献

[1]	Mao, X., Matasov, A. and Piunovskiy, A.B. (2000) Stochastic Differential Delay Equations with Markovian Switching. Bernoulli, 1, 73-90. https://doi.org/10.2307/3318634
[2]	Torre, V. (1977) The Existence of Limit Cycles and Control in Complete Keynesian Systems by Theory of Bifurcations. Econometrica, 45, 1457-1466. https://doi.org/10.2307/1912311
[3]	Turnovsky, S.J. (1996) Methods of Macroeconomic Dynamics. MIT Press, Boston.
[4]	Bing, L., Dingshi, L. and Daoyi, X. (2013) Stability Analysis for Impulsive Stochastic Delay Differential Equations with Markovian Switching. Journal of the Franklin Institute, 350, 1848-1864.

1. 引言

随机微分方程一般常用来描述经济系统的不确定性因素的影响，而如自然灾害的突然发生，或者经济政策的改变，或者企业生产技术的革新等会带来经济系统结构性的变化，因此具有离散状态特征的马氏链也会引入到随机微分方程中。

所以，具有马尔可夫链调制的随机微分方程与随机时滞微分方程的研究大量应用于经济、金融、控制等领域，其相应的稳定性理论、控制理论与数值仿真得到了极大的发展。

在标准的IS-LM经济周期模型应用方面，马氏链调制的随机IS-LM经济周期模型的研究尚少，尤其是模型中考虑到税收过程中存在时滞的情形，研究更少。本文考虑两状态的马尔可夫链到随机时滞IS-LM模型中，研究它的随机稳定性并给出数值仿真例子。

2. 预备知识及模型构建

Mao (2001)研究了如下的马尔可夫调制的随机微分方程 [1] ：

$\text{d}x\left(t\right)=f\left(x\left(t\right),x\left(t-\tau \right),t,r\left(t\right)\right)\text{d}t+g\left(x\left(t\right),x(t-\tau ),t,r\left(t\right)\right)\text{d}w(\; t\; )$

这里 $t\ge 0$ ，具有初始值 $x_0\in C_{F_0}^b\left(\Omega ;R^n\right),f:R^n\times R_{+}\times S\to R^n,g:R^n\times R_{+}\times S\to R^{n\times m}$ 。

假如A是一个向量，则用 $|\cdot |$ 表示Rⁿ中欧氏范数。假如A是一个矩阵(用A^T表示A的转置)，则定义 $|A|=\sqrt{trace\left(A^{\text{T}}A\right)}$ 为迹范数，定义 $\Vert A\Vert =\sup \left\{|Ax|:|x|=1\right\}$ 为算子范数。设 $r\left(t\right),t\ge 0,$ 是定义在完备概率空间 $\left(\Omega ,F,{\left\{F_t\right\}}_{t\ge 0},P\right)$ 上的右连续马尔可夫链，取两个值的状态空间 $S=\left\{1,2\right\}$ ，且具有生成元 $\Gamma ={\left({\gamma }_{ij}\right)}_{2\times 2}$ ，其中

$P\left\{r\left(t+\Delta \right)=j|r\left(t\right)=i\right\}=\{\begin{array}{l}{\gamma }_{ij}\Delta +o\left(\Delta \right),\\ 1+{\gamma }_{ij}\Delta +o\left(\Delta \right),\end{array}$ $\begin{array}{l}i\ne j\\ i=j\end{array}$

这里 $\Delta >0$ ， ${\gamma }_{ij}\ge 0$ 是由i状态转移到j状态的转移速率，且 ${\gamma }_{ii}=-{\displaystyle \underset{j}{\sum }{\gamma }_{ij}},i\ne j$ 。

马尔可夫链 $r\left(t\right)$ 与布朗运动 $w(•)$ 假设是相互独立的，则 $r\left(t\right)$ 的每一个样本轨迹都是一个右连续的阶梯函数，且在每一个有限子区间上只有有限个跳跃点。

V. Torre (1977)用分岔理论分析了非均衡市场的完全凯恩斯系统 [2]

$\{\begin{array}{l}\stackrel{\dot{}}{Y}\left(t\right)=\alpha \left(I\left(Y\left(t\right),R\left(t\right)\right)-S\left(Y\left(t\right),R\left(t\right)\right)\right),\\ \stackrel{\dot{}}{R}\left(t\right)={\beta }_0\left(L\left(Y\left(t\right),R\left(t\right)\right)-L_s\right),\end{array}$

其中 $\frac{\partial I}{\partial Y}>0,\frac{\partial I}{\partial R}<0,\frac{\partial S}{\partial Y}>0,\frac{\partial S}{\partial R}<0,\frac{\partial L}{\partial Y}>0,\frac{\partial L}{\partial R}<0$ ；Y，R，I(Y,R)，S(Y,R)，L(Y,R)，Ls， $\alpha >0$ 和 $\beta >0$ 分别表示产品供给函数、利率、投资(产出和利率的函数)、储蓄(产出和利率的函数)、货币流动需求量、常数货币流动供给量、产品市场和货币市场调整参数。第一个方程表示产品市场的非均衡调整，第二个表示货币市场的非均衡调整。

本文中的第一个方程利用萨谬尔森的动态乘数-加速数模型 [3] 的第一个方程即 $\stackrel{\dot{}}{Y}\left(t\right)=\alpha \left(D\left(t\right)-Y\left(t\right)\right)$ ，其中需求函数 $D\left(t\right)=C\left(t\right)+I\left(t\right)+G$ 表示成消费函数、投资函数及政府支出的和。而消费函数 $C\left(t\right)=\sigma Y^d\left(t\right)=\sigma \left(Y\left(t\right)-T\left(t\right)\right)$ 是关于可支配收入的函数。税收函数是产品供给函数的正常数倍，且税收过程中存在时间延迟现象，因此引入时滞 $\tau$ ，税收函数为 $T\left(t\right)=\epsilon Y\left(t-\tau \right)$ ；投资函数是关于供给和利率的函数 $I\left(Y,R\right)=\gamma Y\left(t\right)-{\beta }_{1}R\left(t\right)$ ；G是常数政府支出。

本文第二个方程中假设流动性需求函数 $L\left(Y\left(t\right),R\left(t\right)\right)=uY\left(t\right)-{\beta }_{2}R\left(t\right)$ ，供给函数 $L_s$ 是常数。则考虑的时滞IS-LM模型为

$\{\begin{array}{l}\stackrel{\dot{}}{Y}\left(t\right)=\alpha \left(\left(\sigma +\gamma -1\right)Y\left(t\right)-\sigma \epsilon Y\left(t-\tau \right)-{\beta }_{1}R\left(t\right)+G\right),\\ \stackrel{\dot{}}{R}\left(t\right)=\beta {}_0\left(uY\left(t\right)-{\beta }_{2}R\left(t\right)-L_s\right)\end{array}$ (1)

这里参数 $0<\alpha ,\sigma ,\gamma ,\epsilon ,{\beta }_0,{\beta }_1,{\beta }_2<1$ 。

命题1：在条件 $u{\beta }_1-{\beta }_2\left(\sigma \left(1-\epsilon \right)-1+\gamma \right)>0$ 和 $uG+\left(\sigma \left(1-\epsilon \right)-1+\gamma \right)L_s>0$ 都满足的情况下，标准时滞IS-LM模型(1)存在唯一的正均衡解

$\left(Y^{*},R^{*}\right)=\left(\frac{{\beta }_{2}G+{\beta }_1{L}_s}{u{\beta }_1-{\beta }_2\left(\sigma \left(1-\epsilon \right)-1+\gamma \right)},\frac{uG+\left(\sigma \left(1-\epsilon \right)-1+\gamma \right)L_s}{u{\beta }_1-{\beta }_2\left(\sigma \left(1-\epsilon \right)-1+\gamma \right)}\right)$ 。

令 $x\left(t\right)=Y\left(t\right)-Y^{*},x\left(t-\tau \right)=Y\left(t-\tau \right)-Y^{*}$ ， $y\left(t\right)=R\left(t\right)-R^{*}$ ，则模型(1)变成

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=\alpha \left(\left(\sigma +\gamma -1\right)x\left(t\right)-\sigma \epsilon x\left(t-\tau \right)-{\beta }_{1}y\left(t\right)\right),\\ \stackrel{\dot{}}{y}\left(t\right)={\beta }_0\left(ux\left(t\right)-{\beta }_{2}y\left(t\right)\right)\end{array}$ (2)

由于此系统受到随机扰动的影响，令变量 $x\left(t\right)$ ， $x\left(t-\tau \right)$ 和 $y\left(t\right)$ 对应的参数在两个方程中分别受到强度为 ${\sigma }_{1}w\left(t\right)$ ， ${\sigma }_{2}w\left(t\right)$ ， ${\sigma }_{3}w\left(t\right)$ ， ${\sigma }_{4}w\left(t\right)$ 和 ${\sigma }_{5}w\left(t\right)$ 的随机扰动。

即模型(2)变为

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=\left(\alpha \left(\sigma +\gamma -1\right)+{\sigma }_1\omega \left(t\right)\right)x\left(t\right)+\left(-\alpha \sigma \epsilon +{\sigma }_2\omega \left(t\right)\right)x\left(t-\tau \right)+\left(-\alpha {\beta }_1+{\sigma }_{3}w\left(t\right)\right)y\left(t\right),\\ \stackrel{\dot{}}{y}\left(t\right)=\left(\left({\beta }_{0}u+{\sigma }_4\omega \left(t\right)\right)x\left(t\right)+\left(-{\beta }_0{\beta }_2+{\sigma }_{5}w\left(t\right)\right)y\left(t\right)\right)\end{array}$ (3)

用矩阵表示为

$\text{d}X\left(t\right)=\left(AX\left(t\right)+BX\left(t-\tau \right)\right)dt+\left(CX\left(t\right)+DX\left(t-\tau \right)\right)\text{d}w\left(t\right)$ (4)

这里 $X\left(t\right)=\left(\begin{array}{l}x\left(t\right)\\ y\left(t\right)\end{array}\right),X\left(t-\tau \right)=\left(\begin{array}{l}x\left(t-\tau \right)\\ y\left(t-\tau \right)\end{array}\right),$

$A=\left(\begin{array}{cc}\alpha \left(\sigma +\gamma -1\right)& -\alpha {\beta }_1\\ {\beta }_{0}u& -{\beta }_0{\beta }_2\end{array}\right),B=\left(\begin{array}{cc}-\alpha \sigma \epsilon & 0\\ 0& 0\end{array}\right)$ ， $C=\left(\begin{array}{cc}{\sigma }_1& {\sigma }_3\\ {\sigma }_4& {\sigma }_5\end{array}\right),D=\left(\begin{array}{cc}{\sigma }_2& 0\\ 0& 0\end{array}\right)$ 。

考虑到上涨与紧缩周期性波动的影响，将二值离散状态的马尔可夫链 $r\left(t\right)$ 引入到模型(4)，则模型变成

$dX\left(t\right)=\left(A\left(r\left(t\right)\right)X\left(t\right)+B\left(r\left(t\right)\right)X\left(t-\tau \right)\right)\text{d}t+\left(C\left(r\left(t\right)\right)X\left(t\right)+D\left(r\left(t\right)\right)X\left(t-\tau \right)\right)\text{d}w\left(t\right)$ (5)

这里r(t) = 1，2。生成元 $\Gamma =\left(\begin{array}{cc}{\gamma }_{11}& {\gamma }_{12}\\ {\gamma }_{21}& {\gamma }_{22}\end{array}\right),{\gamma }_{11}=-{\gamma }_{12}<0,{\gamma }_{22}=-{\gamma }_{21}<0.$

3. 案例仿真

例令

${\alpha }_1=0.011,{\sigma }_1=0.018,{\gamma }_1=0.049,{\beta }_{01}=0.024,{\beta }_{11}=0.056,{\beta }_{21}=0.035,u_1=0.046$ ，

${\epsilon }_1=0.021,{\sigma }_{11}=0.048,{\sigma }_{21}=0.024,{\sigma }_{31}=0.015,{\sigma }_{41}=0.071,{\sigma }_{51}=0.171$ ，

${\alpha }_2=0.011,{\sigma }_2=0.016,{\gamma }_2=0.029,{\beta }_{02}=0.025,{\beta }_{12}=0.065,{\beta }_{22}=0.065,u_2=0.026$ ，

${\epsilon }_2=0.031,{\sigma }_{12}=0.058,{\sigma }_{22}=0.016,{\sigma }_{32}=0.025,{\sigma }_{42}=0.061,{\sigma }_{52}=0.271$

及生成元 $\Gamma =\left(\begin{array}{cc}-3& 3\\ 1& -1\end{array}\right)$ 。

命题1条件 $\{\begin{array}{l}{u}_1{\beta }_{11}-{\beta }_{21}\left({\sigma }_1\left(1-{\epsilon }_1\right)-1+{\gamma }_1\right)>0,\\ u_2{\beta }_{12}-{\beta }_{22}\left({\sigma }_2\left(1-{\epsilon }_{21}\right)-1+{\gamma }_2\right)>0\end{array}$ 满足，则模型(5)存在唯一平凡均衡解。

通过Euler-Maruyama方法 [4] 对案例进行MATLAB仿真，分别在两个不同随机样本和不同时滞( $\tau =0.2,0.8$ )下，得到图1，图2，图3，图4。

图1，图2描述的是在时滞为0.2的情况下两种不同随机样本下，模型(5)的解在时间t < 50以内波动性非常大。但随着时间t的增加(t > 150)模型(5)的解在不断地趋近于平凡均衡解。总体上呈现随机稳定。

图3，图4描述的是模型(5)在时滞为0.8的两种不同随机样本情况。图3中，模型(5)的解在时间t < 50以内波动性非常大，但随着时间t的增加(t > 150)在不断地趋近于均衡解。而图4中，模型(5)的解在时间t < 150以内波动性都非常大，而在t > 150后波动性较少且不断地趋近于均衡解。总体上都呈现随机稳定。

Figure 1. Solution x(t), y(t) of sde and r(t) (time lag = 0.2)

图1. 随机微分方程解x(t)，y(t)及马氏链r(t) (时滞 = 0.2)

Figure 2. Solution x(t), y(t) of sde and r(t) (time lag = 0.2)

图2. 随机微分方程解x(t)，y(t)及马氏链r(t) (时滞 = 0.2)

Figure 3. Solution x(t), y(t) of sde and r(t) (time lag = 0.8)

图3. 随机微分方程解x(t)，y(t)及马氏链r(t) (时滞 = 0.8)

Figure 4. Solution x(t), y(t) of sde and r(t) (time lag = 0.8)

图4. 随机微分方程解x(t)，y(t)及马氏链r(t) (时滞 = 0.8)

从图1~图4看出，r(t)是一个连续的阶梯函数，且在固定的子区间里是有限跳跃的。

基金项目

教育部人文社会科学研究一般项目(11YJCZH259)。