0 引言

在宏观经济问题研究中, 影响经济运行的因素很多, 且各因素联系紧密.而模型和实际问题间存在误差，这些误差在模型中主要以两种形式存在:一是作为模型的随机扰动项, 二是隐含在参数中, 此时参数就是随机参数.因此, 研究随机参数对经济周期模型响应的影响有重要的理论意义和应用价值.

目前, 研究含有界随机参数模型的方法主要有Monte Carlo方法^[1-2], 随机有限元法^[3-4]和正交多项式逼近法^[5-13].Monte Carlo方法简单但费时, 随机有限元法只能解决随机参数是一个小量的情况, 正交多项式逼近法则适用性更强.正交多项式逼近法应用较多的是第二类Chebyshev多项式.但并不是任何情况下都可以使用第二类Chebyshev多项式对模型进行近似逼近, 应根据实际情况选择最合适的正交多项式.另外, 在现有的研究成果中, 很少见到对多项式近似逼近误差的讨论.没有考虑近似误差的多项式逼近有可能是无效的, 因此, 本文根据模型参数的实际意义, 选取Legendre多项式作为正交多项式, 以非线性经济周期模型为例, 在均方误差最小的准则下讨论了参数取值对Legendre多项式近似逼近误差的影响.

1 含有界随机参数的非线性经济周期模型

根据Puu的投资函数^[14]和Goodwin的消费函数^[15], 建立如下经济周期模型

$ \ddot x + \upsilon {{\dot x}^3} + u\dot x + \left( {1-\alpha } \right)x = f\cos \omega t. $

(1)

其中：x表示国民收入增长率, α(0≤α≤1) 是边际消费率, 表示资本市场的供需关系; v是边际投资率, 满足v=1/(1-a), u=2α-1+1/(1-α); f与ω均为无量纲参数, 分别表示周期噪声的强度和频率.由于不同时间的边际消费率不同, 因此α是一个有界随机参数.

设随机参数α=α+σξ, 其中, α是α的均值, σ是α的标准差, ξ是个有界的随机变量, 且|ξ| < 1.

2 正交多项式基的选取(Legendre多项式)

由于ξ的分布没有先验信息, 一般可以将ξ看作服从(-1, 1) 上均匀分布的随机变量, 即ξ的概率密度函数为

$ \rho \left( \xi \right) = \left\{ \begin{array}{l} 1/2, \;\;\;\;\left| \xi \right| < 1, \\ \;\;0, \;\;\;\;\;\left| \xi \right| \ge 1. \end{array} \right. $

(2)

基于上述的概率密度函数, 本文选取Legendre多项式为正交多项式基.这类多项式的一般表达式为

$ {P_\lambda }\left( \xi \right) = \sum\limits_{\kappa = 0}^{\left[{\lambda /2} \right]} {{{\left( { -1} \right)}^\kappa }} \frac{{\left( {2\lambda -2\kappa } \right)!}}{{{2^\lambda }\kappa !\left( {\lambda -\kappa } \right)!\left( {\lambda - 2\kappa } \right)!}}{\xi ^{\lambda - 2\kappa }}. $

(3)

由此可以得到Legendre多项式的具体表达式，即

$ {P_0}\left( \xi \right) = 1, {P_1}\left( \xi \right) = \xi, {P_2}\left( \xi \right) = \left( {3{\xi ^2}-1} \right)/2, {P_3}\left( \xi \right) = \left( {5{\xi ^3}-3\xi } \right)/2, \cdots $

(4)

由Legendre多项式三项式的递推公式, 即

$ \left( {l + 1} \right){P_{l + 1}}\left( \xi \right) = \left( {2l + 1} \right)\xi {P_l}\left( \xi \right)-l{P_{l-1}}\left( \xi \right), \left( {l \ge 1} \right), $

可以得到

$ \xi {P_l}\left( \xi \right) = \left[{\left( {l + 1} \right){P_{l + 1}}\left( \xi \right) + l{P_{l-1}}\left( \xi \right)} \right]/\left( {2l + 1} \right), \left( {l \ge 1} \right). $

(5)

另外, Legendre多项式的正交性可表示为

$ \int_{-1}^1 {{P_k}\left( \xi \right){P_l}\left( \xi \right)/2{\rm{d}}\xi = \left\{ \begin{array}{l} 1/\left( {2l + 1} \right), \;\;\;\;k = l, \\ 0, \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;k \ne l. \end{array} \right.} $

(6)

3 Legendre多项式逼近

为方便后续推导, 给式(1) 两边同时乘以(1-α), 得到

$ \left( {1-\alpha } \right)\ddot x + {{\dot x}^3} + \left( {-2{\alpha ^2} + 3\alpha } \right)\dot x + {\left( {1-\alpha } \right)^2}x = \left( {1 - \alpha } \right)f\cos \omega t. $

(7)

模型(7) 的响应是关于时间t和ξ的函数

$ x = x\left( {t, \xi } \right). $

(8)

由随机函数的正交分解法, 模型(7) 的响应可表示为

$ x\left( {t, \xi } \right) = \sum\limits_{l = 0}^\infty {{x_l}\left( t \right){P_l}\left( \xi \right)} . $

其中, P_l(ξ) 是第l个Legendre多项式.在实际计算中, 选取满足计算精度的有限项近似, 即

$ x\left( t,\xi \right)\approx \sum\limits_{l=0}^{N}{{{x}_{l}}\left( t \right){{P}_{l}}\left( \xi \right)}. $

(9)

以N=4为例给出详细计算过程, N等于其他值的情况可以做类似推导.当N=4时, 将式(9) 代入式(7), 得到

$ \begin{array}{l} \left( {1-\alpha } \right)\sum\limits_{l = 0}^4 {{{\ddot x}_l}\left( t \right){P_l}\left( \xi \right)} + {\left( {\sum\limits_{l = 0}^4 {{{\dot x}_l}\left( t \right){P_l}\left( \xi \right)} } \right)^3} + \left( {-2{\alpha ^2} + 3\alpha } \right)\sum\limits_{l = 0}^4 {{{\dot x}_l}\left( t \right){P_l}\left( \xi \right)} + \\ {\left( {1-\alpha } \right)^2}\sum\limits_{l = 0}^4 {{{\dot x}_l}\left( t \right){P_l}\left( \xi \right)} = \left( {1 - \alpha } \right)f\cos \omega t. \end{array} $

(10)

根据公式(5), 有

$ \begin{array}{l} \left( {1- \alpha } \right)\sum\limits_{l = 0}^4 {{{\ddot x}_l}\left( t \right){P_l}\left( \xi \right)} = \sum\limits_{l = 0}^4 {{{\ddot x}_l}\left( t \right)} \left[{-\frac{{\sigma \left( {l + 1} \right)}}{{2l + 1}}{P_{l + 1}}\left( \xi \right) + \left( {1-\bar \alpha } \right){P_l}\left( \xi \right)-\frac{{\sigma l}}{{2l + 1}}{P_{l - 1}}\left( \xi \right)} \right], \\ \left( { - 2{\alpha ^2} + 3\alpha } \right)\sum\limits_{l = 0}^4 {{{\dot x}_l}\left( t \right){P_l}\left( \xi \right)} = \sum\limits_{l = 0}^4 {{{\dot x}_l}\left( t \right)} \left\{ { - \frac{{2{\sigma ^2}{{\left( {l + 1} \right)}^2}}}{{{{\left( {2l + 1} \right)}^2}}}{P_{l + 2}}\left( \xi \right) + \frac{{\left( { - 4\bar \alpha + 3} \right)\sigma \left( {l + 1} \right)}}{{2l + 1}}{P_{l + 1}}\left( \xi \right)} \right. + \\ \left. {\left[{\left( {-2{{\bar \alpha }^2} + 3\bar \alpha } \right)-\frac{{4{\sigma ^2}\left( {l + 1} \right)l}}{{{{\left( {2l + 1} \right)}^2}}}} \right]{P_l}\left( \xi \right) + \frac{{\left( { - 4\bar \alpha + 3} \right)\sigma l}}{{2l + 1}}{P_{l - 1}}\left( \xi \right) - \frac{{2{\sigma ^2}{l^2}}}{{{{\left( {2l + 1} \right)}^2}}}{P_{l - 2}}\left( \xi \right)} \right\}, \\ {\left( {1 - \alpha } \right)^2}\sum\limits_{l = 0}^4 {{x_l}\left( t \right){P_l}\left( \xi \right)} = \sum\limits_{l = 0}^4 {{x_l}\left( t \right)} \left\{ {\frac{{{\sigma ^2}{{\left( {l + 1} \right)}^2}}}{{{{\left( {2l + 1} \right)}^2}}}{P_{l + 2}}\left( \xi \right) - \frac{{2\sigma \left( {1 - \bar \alpha } \right)\left( {l + 1} \right)}}{{2l + 1}}} \right.{P_{l + 1}}\left( \xi \right) + \\ \left. {\left[{\left( {1-2\bar \alpha + {{\bar \alpha }^2}} \right)-\frac{{2{\sigma ^2}\left( {l + 1} \right)l}}{{{{\left( {2l + 1} \right)}^2}}}} \right]{P_l}\left( \xi \right) -\frac{{2\sigma \left( {1 -\bar \alpha } \right)l}}{{2l + 1}}{P_{l -1}}\left( \xi \right) + \frac{{{\sigma ^2}{l^2}}}{{{{\left( {2l + 1} \right)}^2}}}{P_{l - 2}}\left( \xi \right)} \right\}, \\ \left( {1 - \alpha } \right)f\cos \omega t = \left( {1 - \bar \alpha } \right)f\cos \omega t - \sigma \xi f\cos \omega t. \end{array} $

利用公式(5) 可以得到非线性项$ {{\left( \sum\limits_{l=0}^{4}{{{{\dot{x}}}_{l}}\left( t \right){{P}_{l}}\left( \xi \right)} \right)}^{3}}\approx \sum\limits_{l=0}^{12}{{{X}_{l}}{{P}_{l}}\left( \xi \right)}$, 其中X_l为p_l(ξ) 的系数.

式(10) 两端同乘以p₀(ξ), p₁(ξ), p₂(ξ), p₃(ξ), p₄(ξ), 关于ξ求期望得模型(7) 的近似确定性系统(11), 即

$ \begin{array}{*{20}{c}} {\left( {1-\bar \alpha } \right){{\ddot x}_0}\left( t \right)-\frac{\sigma }{3}{{\ddot x}_1}\left( t \right) + \left( {-2{{\bar \alpha }^2} + 3\bar \alpha } \right){{\dot x}_0}\left( t \right) + \frac{{\left( { - 4\bar \alpha + 3} \right)\sigma }}{3}{{\dot x}_1}\left( t \right) - \frac{8}{{25}}{\sigma ^2}{{\dot x}_2}\left( t \right) + }\\ {{{\left( {1 - \bar \alpha } \right)}^2}{x_0}\left( t \right) - \frac{{2\sigma \left( {1 - \bar \alpha } \right)}}{3}{x_1}\left( t \right) + \frac{4}{{25}}{\sigma ^2}{x_2}\left( t \right) + {X_0} = \left( {1 - \bar \alpha } \right)f\cos \omega t, } \end{array} $

(11-1)

$ \begin{array}{l} - \sigma {{\ddot x}_0}\left( t \right) + \left( {1- \bar \alpha } \right){{\ddot x}_1}\left( t \right)- \frac{{2\sigma }}{5}{{\ddot x}_2}\left( t \right) + \left( { - 4\bar \alpha + 3} \right)\sigma {{\dot x}_0}\left( t \right) + \\ \;\;\;\;\;\;\left[{\left( {-2{{\bar \alpha }^2} + 3\bar \alpha } \right)-8{\sigma ^2}/9} \right]{{\dot x}_1}\left( t \right) + \frac{{\left( { - 4\bar \alpha + 3} \right)2\sigma }}{5}{{\dot x}_2}\left( t \right) - \frac{{18}}{{49}}{\sigma ^2}{{\dot x}_3}\\ - 2\sigma \left( {1 - \bar \alpha } \right){x_0}\left( t \right) + \left[{{{\left( {1-\bar \alpha } \right)}^2} + \frac{4}{9}{\sigma ^2}} \right]{x_1}\left( t \right) -\frac{4}{5}\left( {1 -\bar \alpha } \right){x_2}\left( t \right) + \\ \;\;\;\;\;\;\;\frac{9}{{49}}{\sigma ^2}{x_3}\left( t \right) + {X_1} = \sigma f\cos \omega t, \end{array} $

(11-2)

$ \begin{array}{l} - \frac{{2\sigma }}{3}{{\ddot x}_1}\left( t \right) + \left( {1- \bar \alpha } \right){{\ddot x}_2}\left( t \right)- \frac{{3\sigma }}{7}{{\ddot x}_3}\left( t \right) - 2{\sigma ^2}{{\dot x}_0}\left( t \right) + \frac{{2\left( { - 4\bar \alpha + 3} \right)\sigma }}{3}{{\dot x}_1}\left( t \right) + \\ \left[{\left( {-2{{\bar \alpha }^2} + 3\bar \alpha } \right)-\frac{{24{\sigma ^2}}}{{25}}} \right]{{\ddot x}_2}\left( t \right) + \frac{3}{7}\left( { - 4\bar \alpha + 3} \right)\sigma {{\dot x}_3}\left( t \right)\\ + \frac{{32{\sigma ^2}}}{{81}}{{\dot x}_4}\left( t \right) + {\sigma ^2}{x_0}\left( t \right) - \frac{4}{3}\sigma \left( {1 - \bar \alpha } \right){x_1}\left( t \right) + \\ \left[{{{\left( {1-\bar \alpha } \right)}^2} + \frac{{12}}{{25}}{\sigma ^2}} \right]{x_2}\left( t \right) -\frac{6}{7}\sigma \left( {1 -\bar \alpha } \right){x_3}\left( t \right) + \frac{{16}}{{81}}{\sigma ^2}{x_4}\left( t \right) + {X_2} = 0, \end{array} $

(11-3)

$ \begin{array}{l} - \frac{{3\sigma }}{5}{{\ddot x}_2}\left( t \right) + \left( {1- \bar \alpha } \right){{\ddot x}_3}\left( t \right)- \frac{{8\sigma }}{9}{{\ddot x}_4}\left( t \right) - \frac{4}{9}{\sigma ^2}{{\dot x}_1}\left( t \right) + \frac{{3\left( { - 4\bar \alpha + 3} \right)\sigma }}{5}{{\dot x}_2}\left( t \right) + \\ \left[{\left( {-2{{\bar \alpha }^2} + 3\bar \alpha } \right)-\frac{{48{\sigma ^2}}}{{49}}} \right]{{\dot x}_3}\left( t \right) + \frac{4}{9}\left( { - 4\bar \alpha + 3} \right)\sigma {{\dot x}_4}\left( t \right)\\ + \frac{4}{9}{\sigma ^2}{x_1}\left( t \right) - \frac{6}{5}\left( {1 - \bar \alpha } \right)\sigma {x_2}\left( t \right) + \left[{{{\left( {1-\bar \alpha } \right)}^2} + \frac{{24}}{{49}}{\sigma ^2}} \right]{x_3}\left( t \right) -\frac{8}{9}\sigma \left( {1 -\bar \alpha } \right){x_4}\left( t \right) + {X_3} = 0, \end{array} $

(11-4)

$ \begin{array}{l} - \frac{{4\sigma }}{7}{{\ddot x}_3}\left( t \right) + \left( {1- \bar \alpha } \right){{\ddot x}_4}\left( t \right)- \frac{{18}}{{25}}{\sigma ^2}{{\dot x}_2}\left( t \right) + \frac{{4\left( { - 4\bar \alpha + 3} \right)\sigma }}{7}{{\dot x}_3}\left( t \right) + \\ \left[{\left( {-2{{\bar \alpha }^2} + 3\bar \alpha } \right)-\frac{{80{\sigma ^2}}}{{81}}} \right]{{\dot x}_4}\left( t \right) + \frac{9}{{25}}{\sigma ^2}{x_2}\left( t \right) - \frac{8}{7}\left( {1 - \bar \alpha } \right)\sigma {x_3}\left( t \right) + \\ \left[{{{\left( {1-\bar \alpha } \right)}^2} + \frac{{40}}{{81}}{\sigma ^2}} \right]{x_4}\left( t \right) + {X_4} = 0. \end{array} $

(11-5)

4 Legendre多项式逼近的均方误差

一个好的近似逼近应该有较小的近似误差.由于原模型和近似模型的解析解难以得到, 因此本文用数值解分析参数取值对近似误差的影响.表 1~5列出了N=2, 3, 4, 5时近似确定性模型(11) 与原模型(1) 数值解的均方误差.

由表 1可以看出, 当模型(1) 没有周期性噪声扰动, 即f=0时, 不论Legendre近似逼近项数取多少, 在其他参数不变的情况下, 随着σ的增大, Legendre近似逼近的均方误差迅速增大.当σ增大10倍时, 均方误差增加了10⁴倍.由此看见σ的取值显著影响着Legendre近似逼近的精度, 因此, 在实际应用中应该选取较小的σ, 以提高近似逼近精度.

由表 2可以看出, f取值相同时, 当其他参数都不变, 不同项数的近似逼近均方差异很小, 最小的差异在10⁵倍, 最大的差异在1.04倍.但当f增大5倍时, 相同项数的近似逼近均方误差最大差异超过22倍, 最小差异也大于6倍.也有个别例外, 如表 3~4中, α=α=0.2时, f从1.5增大到2.5时, 均方误差没有增大, 反而有较小程度的变小.总体上, f的变化对均方误差的影响显著.但由于f是模型的周期扰动强度, 在实际问题分析中, 不能随意选取, 因此, 应该考虑在允许的范围内选取最小的f值.

由表 5可以看出, 在其他参数不变的情况下, 周期扰动的频率ω对均方误差的影响比较小.综合比较表 1~5可以看出:(1) 均方误差随着参数的变化而变化; (2) 不论模型是否受到周期噪声的扰动, 在其他参数不变的情况下, σ越大均方误差越大, 因此应选取尽可能小的σ, 以提高近似逼近精度；(3) 在其他参数不变的情况下, α越大误差越小.由于在实际问题中α是个确定值, 不能通过调整α取值提高近似逼近精度.且在参数α取值不同的情况下, 均方误差最小的近似逼近阶数不同.例如在表 5中, 当α=α=0.8, σ=0.2, f=2.5时, ω=1.43时N=5的近似逼近均方误差最小, ω=3.43时N=2的近似逼近均方误差最小.但均方误差的差异小于0.03, 因此在近似逼近对精度要求不是很高的情况下, 尽量选择项数较少的近似逼近, 因为项数每增加一阶, 近似逼近方程组的方程个数将增加1个, 每个方程所含的项数也会随之增加, 导致近似逼近的计算量增加.

5 结束语

正交多项式逼近法使用简单而且应用广泛, 但在实际应用中需要注意选择合适的正交多项式和近似项数.本文在均方误差最小的准则下讨论了参数取值对Legendre多项式近似逼近项数选取的影响.数值分析结果发现，模型参数取值对近似逼近误差有显著影响.随机参数的方差越大, 均方误差越大, 模型的周期扰动强度越大, 均方误差越大.因此, 在实际应用中, 应根据参数的取值情况选取均方误差最小的近似逼近项数, 且在满足精度要求的情况下, 为了后续计算简便, 尽可能选取项数较少的近似逼近.