0 引言

为解决“小样本、贫信息”的不确定系统分析、模拟、预测与决策等问题, 我国著名学者邓聚龙教授提出的灰色系统理论^[1], 已成功应用到经济社会活动中的多个领域.多变量灰色预测模型是灰色系统理论中的一个重要预测模型, 该模型通过引入驱动因素与灰作用量共同作用机制, 建立微分方程模拟和预测系统行为序列的发展趋势.鉴于该模型在实际生活中具有较广泛的应用背景, 能够较好地解决多因素影响下主系统行为的预测问题, 许多学者对模型的数乘变换、鲁棒性等性质^[2-3]、参数优化^[4-5]、模型改进与拓展^[6-7]进行了研究, 从而有效地提高了该模型的建模性能.

从GM(1, N)模型的建模过程来看, 其本质上是因子分析模型, 能够对含有多因子的系统进行整体的、全局的、动态的分析, 但因其不具有全信息而限制了其预测能力, 造成外推困难.事实上, 对于改善模型预测精度, 提高GM(1, N)模型的可预测性, 已经取得了一些有价值的成果.丁松等^[8]基于驱动因素变化趋势提出了TGM(1, N)模型; Tien基于卷积积分技术提出GMC(1, N)模型^[9], 并在此基础上提出了一系列改进模型^[10-13]. Hsu^[14]采用遗传算法优化背景值插值系数, 并在台湾地区集成电路产业的预测问题上取得了较好的应用效果.因此, 目前的研究主要集中在系统行为序列和影响因素同期变化的GM(1, N)模型优化与改进方面.然而, 涉及多变量之间影响大小、影响极性及协调程度时, 影响效果总会呈现出一定的延迟现象, 如从投资到产出, 从交通拥堵到城市环境污染形成均具有系统时滞特征.黄继^[15]针对输入和输出存在时滞特征的系统, 提出了GM(1, N, τ)模型, 并成功应用到武汉科技投入与产出的预测当中.在此基础上, 毛树华等^[16]引入分数阶累加算子, 将时滞模型扩展为分数阶累加GM(1, N, τ)模型; 王正新^[17]针对相关变量累加序列变化较大的时滞系统, 提出了时滞GM(1, N)的基本模型, 并推导出一种派生模型; 张可等^[18]通过引入时滞控制因子构建多变量离散模型, 并结合灰色扩维识别方法探索时滞参数的测算, 在农村水环境与区域发展滞后效应测算中取得了较好的应用效果.

虽然现有的研究在一定程度上提升了多变量预测模型的应用范围, 但他们主要侧重于研究同期及滞后τ期的主系统行为的预测问题, 在实际建模中仍然会产生误差较大的情况.分析其原因, 主要是采用离散形式估计参数与利用微分方程进行预测之间存在跳跃性误差.另外, 驱动项时滞作用机制的复杂性和不确定性, 表现为往期相关因素对于当期系统行为序列间的影响存在时滞累积效应.例如: GDP增长与RD经费投入之间的关系, 由于科技投入效益长期性的特点, GDP增长不仅与当年资金投入水平有关, 还与往期资金投入累计效果有着密切关联, 呈现出多因素能量累计产生的综合效力, 这一实际应用背景在以往的时滞模型中没有得到体现.因此, 本文从该应用背景出发, 对传统多变量预测模型建模过程及其缺陷进行详细分析, 充分考虑了驱动因素往期数据对于当期系统行为序列影响的延迟性和时滞累积效果, 提出TDDGM(1, N)模型, 给出模型的参数识别, 并从信息的完备性角度利用经验分析法和粒子群优化算法求解滞后系数, 最后利用本文模型解决了我国高新技术企业产值预测问题, 同时利用DGM(1, 1)、DGM(1, N)模型进行建模.通过结果对比发现, 本文模型能够更合理地呈现系统行为特征, 并对系统行为进行精确模拟和预测.

1 传统多变量预测模型建模及其缺陷分析

定义 1 设原始数据序列为X_i⁽⁰⁾, 其一次累加生成序列为X_i⁽¹⁾, 称

(1)

为GM(1, N)模型.其中: a为系统发展系数, b_ix_i⁽¹⁾(k)为驱动项, b_i为驱动系数, 为模型的参数列.

定义 2 设, 则称

(2)

为GM(1, N)模型的白化方程,

GM(1, N)模型能够依靠自身动态特性, 反映相关因素序列对系统行为变量的同步作用效果, 在已知影响因子变量变化趋势信息的条件下, 能够对系统行为变量进行预测.但是, 该模型在利用差分形式估计模型参数后, 转而通过微分方程进行模拟预测, 难免会产生模型转换导致的计算误差和模型误差.另外, 在经济社会与工程技术领域的实际应用中, 相关因素序列自身以及不同行为序列之间具有复杂的非线性关系, 可能会出现X_i⁽¹⁾(i=2, 3, ..., N)具有较大变化幅度的情况, 此时驱动作用项不能作为灰常量进行模拟和预测, 从而该模型将不能进行实际预测, 因此构建离散多变量灰色模型显得尤为重要.

定义 3  设原始数据序列为X_i⁽⁰⁾, 其一次累加生成序列为X_i⁽¹⁾, 称

(3)

为多变量离散灰色模型, 记为DGM(1, N), 其中β₁, β₂, ..., β_N+1为DGM(1, N)模型的参数.

DGM(1, N)模型实现了差分形式与微分形式的统一, 有效避免了二者间转换的跳跃误差, 并且能够实现多变量模型的可预测性.然而, 该模型未能充分考虑影响因子对于主系统行为序列的滞后作用机制, 默认主系统行为序列与影响因子序列之间是同期作用关系, 不存在滞后效应.但是, 现实生活中各系统不断发展演化, 系统行为与影响因子之间存在错综复杂的作用关系, 各因素因子在不同时期对当期系统行为可能存在一定的滞后影响, 而且其对系统作用强度也不尽相同.因此, 本文在传统多变量离散模型的基础上, 引入滞后效应参数, 构建基于时滞效应的多变量离散灰色预测模型.

另外, 本文对于时滞的理解与文献[15-18]有一定的差别.以往时滞模型的构建, 是基于系统行为序列X₁(k)当期值仅受到相隔周期p的单期影响因子X_i(k-p)(i=2, 3, ..., N)的作用, 更侧重于对时滞周期的求解; 而本文考虑时滞效应的长期性, 即系统行为序列X₁(k)当期值受到当期及前期所有影响因子序列X_i(1)... X_i(k)(i=2, 3, ..., N)的影响.

需要指出的是, 对于以上两种多变量模型以及本文模型的构建, 由于变量较多并存在差异较大的数量级, 与单变量灰色模型相比, 更容易出现矩阵漂移的现象, 进而影响模型的建模精度.因此, 需要在建立模型之前最好对数量级较大的原始序列通过数乘变换进行一定的处理.

2 基于时滞效应的TDDGM(1, N)模型构建及参数求解 2.1 TDDGM(1, N)模型的定义及参数估计

定义 4 设原始序列为X_i⁽⁰⁾={x_i⁽⁰⁾(1), x_i⁽⁰⁾(2), ..., x_i⁽⁰⁾(n)}, i=1, 2, ..., N, 其中X₁⁽⁰⁾为系统行为序列, X₂⁽⁰⁾, X₃⁽⁰⁾, ..., X_N⁽⁰⁾为影响因素序列.一次累加生成序列为X_i⁽¹⁾ = {x_i⁽¹⁾(1), x_i⁽¹⁾(2), ..., x_i⁽¹⁾(n)}, i=1, 2, ..., N, 其中.称

(4)

为基于时滞效应的多变量灰色离散模型, 记为TDDGM(1, N)).其中λ_i+1(0＜λ_i+1＜1) 称为第i+1(i=1, 2, ..., N-1) 个变量的滞后效应参数, 称为滞后效应驱动项, 反映了此时此刻第i+1个驱动项对系统滞后影响效果的累积.滞后系数λ_i^k-j随着时间由远及近呈现递减趋势, 并且其递减速率与λ_i取值相关, 取值越大, 递减越慢, 反之越快. 图 1所示为取k=10, 在λ的不同取值下滞后系数的变化情况描述.该系数一定程度上揭示了相关因素对于系统行为的影响大小随着时间的往后推移在不断变化, 比较贴近实际.当取j=k时, 此时该模型即退化为传统模型, 说明在系统运行过程中, 所有影响因素只对当期系统行为产生强度相当的影响, 不存在滞后因素, 即TDDGM(1, N)与传统DGM(1, N)模型是等价的, 两种模型会产生相同的建模结果.因此可以将DGM(1, N)看成是本文模型的特例.

另外, 本文模型还能运算滞后驱动因素与非滞后因素同时对主系统行为产生作用的情况, 这是传统离散模型所不能实现的.如对农作物产量预测时, 其投入因素农作物播种面积和人力、资本要素投入, 播种面积是非滞后因素, 而人力、资本投入是具有滞后累积效果的因素, 此时构建农作物产量预测时应包含这两类影响因素, 才能更好地体现模型构建的实际应用意义.因此, 在实际考察时滞影响关系时, 需要针对具体应用案例加以分析其影响因素是否具有滞后累积效果, 构建合适的模型进行模拟预测.

定理 1  X_i⁽⁰⁾、X_i⁽¹⁾、λ_i+1(0＜λ_i+1＜1) 如定义4所示, 参数列, 且

由最小二乘法的矩阵运算可得:

1) n=N+2时, ;

2) 当n>N+2时, , |B^TB|≠ 0;

3) 当n＜N+2时, , | B^TB|≠ 0.

证明 将k=2, 3, ..., n代入TDDGM(1, N)模型(4), 可得方程组

即由最小二乘法可得Y=Bβ.

1) 当n=N+2且|B|≠ 0时, B存在逆矩阵, 方程组有唯一解, 即可得到.

2) 当n>N+2且B为列满秩时, 有B的满秩分解为B=DC, 进而可以得到B的广义逆矩阵

由于B为满秩矩阵, C可取为单位矩阵, B=DI_N, B=D, 即

3) 当n＜N+2且B为行满秩矩阵时, D可取为单位矩阵, B=I_n-1C, B=C, 则

定理 2 设B, Y, 如定理1所描述, =[β₁, β₂, ..., β_N+1]^T, 记, 可得:

1) 取初值x₁⁽¹⁾(1)=x₁⁽⁰⁾(1)(k=1, 2, ..., n-1) 时模型的模拟值为

(5)

2) 模型的还原值为

(6)

证明  1) 利用数学归纳法可证明.当k= 1时, 有

结论成立.假设当k=m时结论成立, 即可得

由定义4中的式(4) 可以得到

(7)

将k=m的值代入式(8) 可得到

因此, 当k=m+1时结论也成立.

2) 由累加生成的逆过程可知, , k=1, 2, ..., n-1.

定理 3 当j=k时, 滞后效应驱动项退化为T_i+1(k)=x_i+1⁽¹⁾(k), TDDGM(1, N)模型退化为传统DGM(1, N)模型, 则可得:

1) 取初值x₁⁽¹⁾(1)=x₁⁽⁰⁾(1)(k=1, 2, ..., n-1) 时模型的模拟值为

(8)

2) 模型的还原值为

(9)

证明同定理2, 此处略.

定理1通过最小二乘法构建了TDDGM(1, N)的模型参数估计方法, 给出了不同样本特征情况下参数估计的求解过程; 定理2以系统行为序列的初始值作为模型建模的初始点, 采用迭代的方法建立了时滞系统的主行为序列的模拟和预测体系.当影响因素序列的滞后信息已知时, 依据定理1和定理2可以进行建模模拟和预测.定理3讨论了本文模型与传统多变量离散模型之间的关系, 本文模型是传统模型的延伸和拓展, 在传统多变量离散模型的基础上, 丰富了其实际应用意义和内涵.另外, 当驱动因素白化信息不断丰富或者需要进行中长期预测时, 可以考虑引入新陈代谢控制机制, 去掉最早期的一个数据, 保持数列等维, 然后建立TDDGM(1, N)模型预测下一个数据, 不断重复, 依次递补, 直到完成预测目标为止.该做法不仅符合驱动因素累积时滞效果的时效性特征, 还能简化计算过程, 提升模型精度.

2.2 TDDGM(1, N)时滞参数的识别方法

在基于时滞效应控制的多变量离散灰色预测模型建模过程中, 最核心的部分是其时滞参数的识别与计算, 直接影响模型建模和预测精度.由于多变量系统中结合信息存在灰性, 本文拟从白化信息充分和白化信息匮乏两个角度识别滞后参数的计算方法.

1) 当系统结构信息充分, 驱动因素对于主系统行为的作用方式及效果明确, 或者不存在明显的滞后效应时, 根据经验分析可以选定λ_i的值, 进而代入TDDGM(1, N)估计模型结构参数, 求解模型时间响应函数进行模拟和预测.

2) 当系统行为状态或者影响因素对于系统的作用方式发生变化或者缺乏充分白化信息时, 这也是实际建模过程中经常存在的情况, 此时对于反映各个变量对于主系统影响递减程度的λ_i未知.本文拟从控制模型精度的角度, 利用最小一乘法, 通过模型参数之间的约束关系来构建非线性优化模型, 从而确定各影响因素对主系统延迟的作用机制.

文献[15-18]均以残差平方和最小对参数进行寻优, 与模型的检验标准平均相对误差最小不一致, 两个准则的不一致性可能会造成模型优化效果不理想, 导致模型平均相对误差变劣^[19].残差平方和最小准则侧重于等精度序列, 表现出大致相同的绝对误差; 而相对误差最小准则侧重于非等精度序列, 表现为大致相同的相对误差, 两者内涵存在差异性.鉴于模型精度一般以平均相对误差作为衡量模型拟合效果的指标, 即原始数据越大, 所能接受的拟合绝对误差越大.本文以还原序列的相对误差绝对值之和的平均数作为最优目标函数, 即最小一乘法, 实现最优目标函数与误差检验标准原则的统一, 以期进一步提升模型的精度.

设λ=[λ₂, λ₃, ..., λ_N]分别为各驱动因素的驱动因素的滞后参数, 为TDDGM(1, N)模型的模拟值, x₁⁽⁰⁾(k)为建模的真实值, 记

则时滞参数的求解可以转化为相对误差绝对值和绝对值最小的非线性优化模型.对于TDDGM(1, N)模型, 建立如下所示的非线性优化模型:

(10)

式(10) 可以通过优化软件LINGO或者采用智能优化算法(如粒子群优化算法、遗传算法等)进行求解, 一旦完成滞后参数的最优取值, 然后代入时间响应函数即可计算和预测结果.

2.3 TDDGM(1, N)模型建模有效性检验

除了常用的平均相对误差作为模型建模效果检验标准以外, 本文将通过关联模型从拟合序列相似程度的角度做检验, 对比拟合序列与原始序列的形状相似程度, 以此作为可否进行预测的依据, 保证建模得到的拟合序列与原始序列保持趋势的一致性.采用绝对关联度进行拟合序列和原始序列形状相似性的判断依据, 计算公式如下(X₀⁽⁰⁾和X_j⁽⁰⁾分别为序列X₀和X_j的始点零化像)

则关联度计算公式为

2.4 TDDGM(1, N)模型建模步骤

综上所述, TDDGM(1, N)模型建模与预测结果可以分为如下几个步骤:

Step1:利用获得的白化信息, 分析提炼影响主系统行为的序列, 确定驱动影响因子X_i⁽⁰⁾, i=2, 3, ..., N, 并初步识别这些影响因子是否包含滞后信息, 再根据定义4构建TDDGM(1, N)模型.

Step 2:若已知系统行为信息及滞后参数信息, 则得出滞后参数λ_i, 根据定理1估计模型参数, 再跳转Step 4.

Step 3:若缺乏系统白化信息和滞后效果信息, 则依据式(10) 构造优化模型, 利用作为中间参数, 通过粒子群智能算法获得滞后参数估计值, 求出滞后参数以后, 根据定理1估计出模型参数数.

Step 4:依据定理2计算主系统行为序列的时间响应函数, 并进行模拟和预测分析, 检验模型精度.若模型误差超出容忍范围, 则转回Step1, 重新对影响因子及滞后参数信息进行分析.

Step 5:采用GM(1, 1) 分别预测各影响因素X_i⁽⁰⁾ (i=2, 3, ..., N)的预测值.

Step 6:根据影响因子的预测值, 利用TDDGM(1, N)预测主系统行为序列X₁⁽⁰⁾, 并对预测结果的合理性和有效性进行分析.

3 我国高新技术企业产值预测分析

高新技术企业是提升国家创新能力的重要载体, 是国家之间经济、科技竞争的制高点, 代表了国家未来技术与产业的发展方向.从高新技术企业本身来看, 高新技术企业的产出人力投入和自己投入是创新产出的重要影响因素.在人力资本投入方面, R & D人才是高新技术企业持续发展的基本动力, 拥有高质量的R & D人才和高规模的人力投入, 才能保证企业研发创新活动的高效持续进行, 进而提升企业创新产出.在资金投入方面, R & D经费投入能够反映高新技术企业研发创新的基础和能力及潜在产出, 并为高新技术企业研发创新活动的顺利完成提供基础保障.因此, 增加和合理利用R & D人员及R & D经费, 会直接促进专利授权量, 进而推动产品产值的增加, 提升高新技术企业的创新产出.

为了能验证本文模型在实际应用中的有效性, 本文拟通过多种模型的对比分析, 对我国高新技术企业产值进行预测分析.高新技术企业十分重视技术创新的投入, 这部分主要包括R & D人员投入和R & D经费投入.根据Step1, 首先进行变量的选择, 本文以我国高新技术企业产值为产出指标(系统行为变量), 以R & D人员折合全时当量(万人年)和R & D经费内部支出(亿元)为投入指标(相关因素变量), 这3个指标2005～2014年的数据如表 1所示, 以现价计算, 3个指标在2006～2014年的增长率如图 2所示.

由表 1和图 2可以看出, 2005～2014年间, 我国高新技术产业发展总体趋势较好, 呈现波动上升趋势.另外, 可计算高新技术产业产值、R & D人员和R & D经费的年均增长率分别为15.84 %、14.21 %和20.36 %, 这一方面表明了我国研究与发展人员和经费投入力度较大, 同时也说明了投入要素的增长和产出要素的增长并不同步, 存在一定的差异性.究其原因, 主要是这两种要素的投入到产生效益存在明显的时滞累积特征, 而且, 往期要素投入对当期产值水平均会产生一定的影响.因此, 对于这种少信息多变量且投入与产出具有滞后性的控制系统, 适合用TDDGM(1, N)模型进行建模分析.

由于系统行为序列和驱动因素滞后信息比较缺乏, 根据建模Step 3, 本文采用粒子群算法对时滞参数进行优化, 明确驱动项时滞效应的作用机制.为此, 选用2005～2012年数据作为建模数据, 选用2013～2014年数据作为样本外推预测.由建模结果可知, 模型的各参数结果如表 2所示.

根据Step 4可以建立TDDGM(1, N)模型为

(11)

由表 2可以看出, R & D人员投入和R & D经费投入的滞后系数分别为0.999 9和0.608 3, 表明R & D人员和经费投入对企业产值的促进不能一蹴而就, 存在滞后效应, 其滞后效果变化情况如图 3所示. R & D人员投入对我国高新技术企业产值的影响非常大, 且对企业的影响具有较强的滞后效应, 能够长期影响企业未来发展.而R & D经费投入的滞后影响也较强, 其累计投入效果持续时间相对较短, 在某期经费投入10年时, 经费滞后效果几乎消失.因此, R & D人才和经费投入对于高新技术企业的发展至关重要, 坚持推进R & D活动, 从长远角度看, 企业终将受益.为了对比模型的性能, 分别建立TDDGM(1, N)、DGM(1, 1)、DGM(1, N)模型, 具体如下.

DGM(1, 1) 模型为

(12)

DGM(1, N)模型为

(13)

由表 3可以得出, DGM(1, 1) 模型虽然避免了差分方程到微分方程间的跳跃性误差, 但其为利用驱动因素作用信息, 仅适用于增长率恒定的近似指数规律的系统建模与预测, 对于存在波动性、滞后性特点的序列建模性能较差, 因此模拟和预测平均误差分别为16.02 %和7.79 %. DGM(1, N)模型虽然采用了与本文同样的驱动因素变量, 但其并未进一步分析R & D人员投入和R & D经费投入与高新技术产业产值之间的滞后效应, 更未讨论往期驱动因素对当期高新技术产值还存在明显的累积影响效果, 因此模型存在较大的模拟和预测误差, 分别为31.71 %和38.49 %, 已失去建模意义.本文模型通过控制驱动变量时滞累积效应的作用机制, 不仅建模误差远低于DGM(1, 1) 和DGM(1, N)模型, 而且预测误差也远低于后两类模型, 避免了后两类模型在2013～2014年间出现预测误差跳跃现象, 两项误差分别为5.28%和0.48%. 图 4反映了3种模型的模拟和预测相对误差波动情况.

另外, 从灰色关联方法的有效性检验角度来看, 本文模型与原始序列的绝对关联度为0.943 1, 而传统灰色DGM(1, 1) 模型、DGM(1, N)模型得到的拟合序列与原始序列之间的绝对关联度分别为0.874 9和0.667 1.由此可知, 本文模型相对于其他两种模型对序列特征的描述更为精准, 用于预测的可行性更高.综合采用两种检验方法对比各种模型的建模效果来看, 本文的离散灰色模型通过在驱动项中加入驱动因素及时滞累积效应的作用机制等白化信息, 能够持续降低建模和预测误差, 适用于我国高新技术产业产值等具有滞后累积效果的问题预测.

4 结论

灰色多变量预测模型的本质就是通过挖掘驱动变量与主系统行为变量之间的白化系统信息, 发现系统行为规律, 然后建模并预测.本文基于灰色系统建模方法, 提出了基于时滞效应的离散多变量灰色预测模型, 并给出了其与传统多变量模型间的关系, 探讨了参数估计方法; 然后利用最小一乘法, 结合粒子群寻优算法能够减少平均相对误差变劣现象, 进一步提升模型的建模精度, 为解决具有时滞累计效应特点的小样本系统建模问题提供了有效的分析工具. TDDGM(1, N)模型不仅克服了传统模型的跳跃性误差, 还明确了驱动项及其时滞累计效应作用机制, 能够持续降低模拟和预测误差.实例分析表明, TDDGM(1, N)模型能够有效描述高新技术产业产值与R & D人员和R & D经费间的时滞关系, 并能够显著地提高对高新技术产值的模拟和预测精度.