预测多元时间序列

这个例子展示了如何在捕食者拥挤的情况下，对从捕食者和被捕食者种群中测量的数据进行多变量时间序列预测。利用线性和非线性时间序列模型建立了捕食者-被捕食者种群变化动力学模型。比较了这些模型的预测性能。

数据描述

数据是一个双变量时间序列，包括1-捕食者1-猎物种群(以千为单位)，每年收集10次，持续20年。有关数据的更多信息，请参见三种生态种群系统:MATLAB和C MEX-File时间序列建模．

加载时间序列数据。

负载PredPreyCrowdingDataz = iddata (y, [], 0.1,“TimeUnit”，“年”，“Tstart”, 0);

z是一个iddata对象，它包含两个输出信号，日元而且y2，分别指捕食者和被捕食者的种群。的OutputData的属性z以201 × 2矩阵的形式包含人口数据，这样z.OutputData (: 1)捕食者的数量和z.OutputData (: 2)是猎物的数量。

图数据。

情节(z)标题(“捕食者-猎物人口数据”) ylabel (“人口(千)”）

图中包含2个轴对象。标题为y1的Axes对象1包含一个类型为line的对象。axis对象2的标题为y2，包含一个类型为line的对象。这个对象表示z。

数据显示由于拥挤，捕食者数量下降。

用前半部分作为估计数据来识别时间序列模型。

泽= z (1:120);

利用剩余数据搜索模型订单，并对预测结果进行验证。

zv z =(121:结束);

线性模型的估计

将时间序列建模为线性自回归过程。可以使用以下命令以多项式形式或状态空间形式创建线性模型基于“增大化现实”技术(仅针对标量时间序列)，arx，armax，n4sid而且党卫军．由于线性模型不捕获数据偏移量(非零条件均值)，首先对数据进行趋势化处理。

[zed, Tze] =趋势(ze, 0);[zvd, Tzv] =趋势(zv, 0);

识别需要型号订单的说明。对于多项式模型，您可以使用arxstruc命令。自arxstruc只适用于单输出模型，对每个输出分别执行模型顺序搜索。

na_list = (1:10) ';V1 = arxstruc (zed(: 1:)、zvd (: 1:), na_list);na1 = selstruc (V1 0);V2 = arxstruc (zed(:, 2:)、zvd (:, 2:), na_list);钠= selstruc (V2, 0);

的arxstrucCommand分别提出了7和8阶的自回归模型。

使用这些模型顺序来估计一个交叉项被任意选择的多方差ARMA模型。

Na = [na1 na1-1;na2-1钠);数控= [na1;钠);sysARMA = armax(zed，[na nc])

sysARMA =离散ARMA模型:模型输出“日元”:一个(z) y_1 (t) = - ai (z) y_i (t) + C (z) e_1 (t) (z) = 1 - 0.885 z ^ 1 - 0.1493 z ^ 2 + 0.8089 z ^ 3 - 0.2661 z ^ 4 - 0.9487 z ^ 5 + 0.8719 z ^ 6 - 0.2896 z ^ 7 A₂(z) = 0.3433 z ^ 1 - 0.2802 z ^ 2 - 0.04949 z ^ 3 + 0.1018 z ^ 4 - 0.02683 z ^ 5 - 0.2416 z ^ 6 C (z) = 1 - 0.4534 z ^ 1 - 0.4127 z ^ 2 + 0.7874 z ^ 3 + 0.298 z ^ 4 - 0.8684 z z ^ 6 ^ 5 + 0.6106 + 0.3616 z ^ 7模型输出“y2”:(z) y_2 (t) = - ai (z) y_i (t) + C (z) e_2 (t) (z) = 1 - 0.5826 z ^ 1 - 0.4688 z ^ 2 - 0.5949 z ^ 3 - 0.0547 z ^ 4 + 0.5062 z ^ 5 + 0.4024 z ^ 6 - 0.01544 z ^ 7 - 0.1766 z ^ 8 A_1 (z) = 0.2386 z ^ 1 + 0.1564 z ^ 2 - 0.2249 z ^ 3 - 0.2638 z ^ 4 - 0.1019 z ^ 5 - 0.07821 z ^ 6 + 0.2982 z ^ 7 C (z) = 1 - 0.1717 z ^ 1 - 0.09877 z ^ 2 - 0.5289 z ^ 3 - 0.24 z ^ 4 + 0.06555 z ^ 5 + 0.2217 z ^ 6 - 0.05765 z ^ 7 - 0.1824 z ^ 8样时间:0.1年参数化:多项式订单:na=[7 6;7 8] nc=[7;8]自由系数数量:43使用"polydata"， "getpvec"， "getcov"表示参数及其不确定性。状态:使用ARMAX对时域数据“zed”进行估计。拟合估计数据:[89.85;90.97]%(预测焦点)FPE: 3.814e-05, MSE: 0.007533

计算提前10步(1年)的预测输出，以在估计数据的时间跨度内验证模型。由于数据是为了估计而去趋势化的，您需要为有意义的预测指定那些偏移量。

predOpt = predictOptions (“OutputOffset”Tze.OutputOffset ');yhat1 = predict(sysARMA,ze,10, predOpt);

的预测Command预测被测数据在时间跨度内的响应，是验证估计模型质量的工具。响应时间t用t = 0时的测量值计算，…, t-10。

绘制预测响应和实测数据图。

情节(泽,yhat1)标题(“与测量数据相比，10步预测反应”）

$图中包含2个轴对象。标题为y1的axis对象1包含两个类型为line的对象。这些对象表示ze, ze\ _expected。标题为y2的Axes对象2包含两个类型为line的对象。这些对象表示ze, ze\ _expected。$

注意，预测响应的生成和用测量数据绘制它，可以使用比较命令。

compareOpt = compareOptions (“OutputOffset”Tze.OutputOffset ');比较(泽sysARMA 10 compareOpt)

图中包含2个轴对象。坐标轴对象1包含2个line类型的对象。这些对象表示验证数据(y1)， sysARMA: 75.49%。坐标轴对象2包含两个line类型的对象。这些对象表示验证数据(y2)， sysARMA: 76.21%。

使用比较还以百分比形式显示了归一化均方根(NRMSE)的拟合优度度量。

在验证数据后，预测模型的输出sysARMA100步(10年)以外的估计数据，并计算输出标准差。

forecastOpt = forecastOptions (“OutputOffset”Tze.OutputOffset ');[yf1,x01,sysf1,ysd1] = forecast(sysARMA, ze, 100, forecastOpt);

yf1预测的响应是否返回为iddata对象。yf1。OutputData包含预测值。

sysf1系统类似于sysARMA但它是状态空间形式。模拟sysf1使用sim卡在初始条件下，x01，再现预测响应，yf1．

ysd1是标准差矩阵。它测量了由于数据中附加扰动的影响而导致的预测的不确定性sysARMA。NoiseVariance)，参数不确定度(如getcov (sysARMA))以及与将过去数据映射到预测所需初始条件的过程相关的不确定性(参见data2state)．

绘制模型的测量、预测和预测输出sysARMA．

t = yf1.SamplingInstants;te = ze.SamplingInstants;t0 = z.SamplingInstants;次要情节(1、2、1);情节(t0 z.y (: 1),.．.te, yhat1.y (: 1),.．.t, yf1.y (: 1),“米”，.．.t, yf1.y (: 1) + ysd1 (: 1),“k——”，.．.t, yf1.y (: 1) -ysd1 (: 1),“k——”)包含(的时间(年)）;ylabel (“捕食者的数量，成千上万”）;次要情节(1、2、2);情节(t0 z.y (:, 2),.．.te, yhat1.y (:, 2),.．.t, yf1.y (:, 2),“米”，.．.t, yf1.y (:, 2) + ysd1 (:, 2),“k——”，.．.t, yf1.y (:, 2) -ysd1 (:, 2),“k——”）把数字画大一点。无花果= gcf;p = fig.Position;fig.Position = [p (1), (2) - p (4) * 0.2, p (3) * 1.4, p (4) * 1.2);包含(的时间(年)）;ylabel (“猎物数量，以数千计”）;传奇({“测量”，“预测”，“预测”，“预测不确定性(1sd)”}，.．.“位置”，“最佳”）

图中包含2个轴对象。axis对象1包含5个类型为line的对象。Axes对象2包含5个line类型的对象。这些对象表示测量、预测、预测、预测不确定度(1 sd)。

这些图表明，使用线性ARMA模型(增加了对偏移量的处理)进行预测取得了一定效果，结果显示，与12-20年的实际人口相比，其不确定性很高。这表明种群变化动态可能是非线性的。

估计非线性黑箱模型

在本节中，我们通过在模型中添加非线性元素来扩展上一节的黑箱识别方法。这是使用非线性ARX建模框架进行的。非线性ARX模型在概念上类似于线性模型(如sysARMA)，因为他们使用两步过程计算预测的输出:

将测量信号映射到有限维回归空间。
使用可以是线性或非线性的静态函数将回归量映射到预测输出。

2021 - 03 - 30 - _14 29 - 33. png

非线性可以在两个地方添加——要么在回归函数集中(步骤1)，要么在静态映射函数中(步骤2)。在本节中，我们将探索这两个选项。首先，我们创建了非线性ARX模型的线性回归形式，其中所有的非线性都仅限于回归定义;将回归量转换为预测输出的静态映射函数是线性的(确切地说，是仿射的)。然后，我们建立了一个模型，该模型使用高斯过程函数作为静态映射函数，但其回归量保持线性。

多项式AR模型

创建一个有两个输出，没有输入的非线性ARX模型。

L = linearRegressor (ze.OutputName {1});sysNLARX = idnlarx(泽。OutputName,泽。InputName, L, []);sysNLARX。t = 0.1;sysNLARX。TimeUnit =“年”；

sysNLARX为不使用非线性函数的一阶非线性ARX模型;它使用一阶回归量的加权和来预测模型响应。

getreg (sysNLARX)

ans =2 x1细胞{“y1 (t - 1)”}{}“y2 (t - 1)”

为了引入非线性函数，在模型中加入多项式回归量。

创建二阶回归量并包含交叉项(上面列出的标准回归量的乘积)。2022世界杯八强谁会赢？将这些回归量添加到模型中。

P = polynomialRegressor (ze.OutputName,{1}, 2,假的,真的);sysNLARX.Regressors结束(+ 1)= P;getreg (sysNLARX)

ans =5 x1细胞{y1 (t - 1)的}{y2 (t - 1)的}{y1 (t - 1) ^ 2的}{y2 (t - 1) ^ 2的}{“y1 (t - 1) * y2 (t - 1)}

使用估计数据估计模型的系数(回归系数和偏移量)，泽．

sysNLARX sysNLARX = nlarx(泽)

sysNLARX =具有两个输出的非线性多元时间序列模型输出:y1, y2变量y1 y2的线性回归量。变量y1, y2中的二阶回归量输出函数:输出1:线性带偏移输出2:线性带偏移采样时间:0.1年状态:终止条件:接近(局部)最小值，(范数(g) < tol ..迭代次数:0，函数计算次数:1使用NLARX对时域数据“ze”进行估计。拟合估计数据:[88.34;88.91]%(预测焦点)FPE: 3.265e-05, MSE: 0.01048模型“报告”属性中的更多信息。

计算一个提前10步的预测输出来验证模型。

yhat2 =预测(sysNLARX泽10);

预测模型的输出超过估计数据100步。

yf2 =预测(sysNLARX,泽,100);

非线性ARX模型没有计算预测响应的标准差。这些数据是不可用的，因为在估计这些模型时没有计算参数协方差信息。

绘制测量的、预测的和预测的输出。

t = yf2.SamplingInstants;次要情节(1、2、1);情节(t0 z.y (: 1),.．.te, yhat2.y (: 1),.．.t, yf2.y (: 1),“米”)包含(的时间(年)）;ylabel (“捕食者数量(千)”）;标题(“多项式回归ARX模型的预测结果”)次要情节(1、2、2);情节(t0 z.y (:, 2),.．.te, yhat2.y (:, 2),.．.t, yf2.y (:, 2),“米”)传说(“测量”，“预测”，“预测”)包含(的时间(年)）;ylabel (“猎物种群(千)”）;

图中包含2个轴对象。标题为“ARX模型的多项式回归预测结果”的坐标轴对象1包含3个类型为直线的对象。Axes对象2包含3个line类型的对象。这些对象代表测量，预测，预测。

结果表明，非线性ARX模型的预测效果优于线性模型。非线性黑箱建模不需要关于控制数据的方程的先验知识。

注意，要减少回归量的数量，可以使用主成分分析(参见主成分分析)或特性选择(参见sequentialfs)中的统计和机器学习工具箱™。

估计一个高斯过程模型

模型模型sysNLARX使用线性映射函数;非线性只包含在回归量中。更灵活的方法是选择一个非平凡的映射函数，例如高斯过程函数(GP)。创建一个非线性ARX模型，只使用线性回归量，但添加一个高斯过程函数将回归量映射到输出。

复制sysNLARX，这样我们就不必再创建模型结构了sysGP = sysNLARX;%移除多项式回归集sysGP.Regressors (2) = [];将线性映射函数替换为GP函数。GP默认使用平方指数核。GP = idGaussianProcess;sysGP。OutputFcn =全科医生;disp (sysGP)

[2x2 double] InputName: {0x1 cell} InputUnit: {0x1 cell} InputGroup: [1x1 struct]注释:[0x1 string] UserData: [] Name: " Ts: 0.1000 TimeUnit: 'years' Report: [1x1 idresults.nlarxhw]

现在，您有了一个基于线性回归器和GP映射函数的模板模型sysGP。它的参数是它的平方指数核的两个系数。使用nlarx来估计这些参数。

sysGP = nlarx(ze, sysGP)

sysGP =具有2个输出的非线性多元时间序列模型输出:y1, y2回归量:变量y1, y2中的线性回归量输出函数:输出1:使用SquaredExponential核的高斯过程函数输出2:使用SquaredExponential核的高斯过程函数采样时间:0.1年状态:使用时域数据“ze”上的NLARX估计。拟合估计数据:[88.07;88.84]%(预测焦点)FPE: 3.369e-05, MSE: 0.01083模型“报告”属性中的更多信息。

getpvec (sysGP)

ans =10×1-0.6671 0.0196 -0.0000 -0.0571 -3.6447 -0.0221 -0.5752 0.0000 0.1725 -3.1068

和前面一样，计算预测和预测的响应并绘制结果图。

计算一个提前10步的预测输出来验证模型。

yhat3 =预测(sysGP泽10);

预测模型的输出超过估计数据100步。

yf3 =预测(sysGP,泽,100);

绘制测量的、预测的和预测的输出。

t = yf3.SamplingInstants;次要情节(1、2、1);情节(t0 z.y (: 1),.．.te, yhat3.y (: 1),.．.t, yf3.y (: 1),“米”)包含(的时间(年)）;ylabel (“捕食者数量(千)”）;标题(“高斯过程模型的预测结果”)次要情节(1、2、2);情节(t0 z.y (:, 2),.．.te, yhat3.y (:, 2),.．.t, yf3.y (:, 2),“米”)传说(“测量”，“预测”，“预测”)包含(的时间(年)）;ylabel (“猎物种群(千)”）;

图中包含2个轴对象。标题为“高斯过程模型预测结果”的坐标轴对象1包含3个类型为line的对象。Axes对象2包含3个line类型的对象。这些对象代表测量，预测，预测。

结果表明，与线性模型相比，该模型提高了精度。使用基于GP的模型的一个优点是，与其他形式如小波网络(idWaveletNetwork)或Sigmoid Networks (idSigmoidNetwork)．这使得它们易于训练，且参数估计的不确定性较低。

非线性灰盒模型的估计

如果您有系统动力学的先验知识，您可以使用非线性灰盒模型拟合估计数据。总的来说，了解动态的本质有助于提高模型的质量，从而提高预测的准确性。对于捕食者-被捕食者动态，捕食者的变化(日元)和猎物(y2)可表示为:

${y_{}^{˙}}_{1} （ t ）＝ p_{1} ＊ y_{1} （ t ） + p_{2} ＊ y_{2} （ t ）＊ y_{1} （ t ）$

${y_{}^{˙}}_{2} （ t ）＝ p_{3.} ＊ y_{2} （ t ） - p_{4} ＊ y_{1} （ t ）＊ y_{2} （ t ） - p_{5} ＊ y_{2} （ t ）^{2}$

有关这些方程的更多信息，请参见三种生态种群系统:MATLAB和C MEX-File时间序列建模．

根据这些方程建立非线性灰盒模型。

指定一个描述捕食者-猎物系统模型结构的文件。该文件将状态导数和模型输出指定为时间、状态、输入和模型参数的函数。选择两个输出(捕食者和猎物种群)作为状态，推导出动态的非线性状态空间描述。

文件名=“predprey2_m”；

指定模型顺序(输出、输入和状态的数量)

Order = [2 0 2];

指定参数的初始值 $p_{1}$ ， $p_{2}$ ， $p_{3.}$ ， $p_{4}$ , $p_{5}$ ，表示所有参数都要估计。注意，在估计黑盒模型时，不存在指定参数初始猜测的要求sysARMA而且sysNLARX．

参数=结构(“名字”, {“生存因素,捕食者”“死亡因素,捕食者”.．.“生存因素,猎物”“死亡因素,猎物”.．.“拥挤因素,猎物”}，.．.“单位”, {1 /年的1 /年的1 /年的1 /年的1 /年的}，.．.“价值”，{-1.1 0.9 1.1 0.9 0.2}，.．.“最低”，{-Inf -Inf -Inf -Inf}，.．.“最大”，{Inf Inf Inf Inf Inf}，.．.“固定”，{假假假假假假假});

类似地，指定模型的初始状态，并指出要估计两个初始状态。

InitialStates =结构(“名字”, {捕食者种群的“猎物人口”}，.．.“单位”, {的大小(千)的大小(千)}，.．.“价值”{1.8 - 1.8},.．.“最低”, {0 0},.．.“最大”{正正无穷},.．.“固定”,{假假});

指定模型为连续时间系统。

t = 0;

创建一个具有指定结构、参数和状态的非线性灰盒模型。

sysGrey = idnlgrey(文件名、秩序、参数、InitialStates Ts,“TimeUnit”，“年”）;

估计模型参数。

sysGrey = nlgreyest(泽,sysGrey);礼物(sysGrey)

sysGreyGrey =由'predprey2_m'定义的连续时间非线性灰盒模型(MATLAB文件):dx/dt = F(t, x(t)， p1，…， p5) y(t) = H(t, x(t)， p1，…，p5）+e(t) with 0 input(s), 2 state(s), 2 output(s), and 5 free parameter(s) (out of 5). States: Initial value x(1) Predator population(t) [Size (thou..] xinit@exp1 2.01325 (estimated) in [0, Inf] x(2) Prey population(t) [Size (thou..] xinit@exp1 1.99687 (estimated) in [0, Inf] Outputs: y(1) y1(t) y(2) y2(t) Parameters: Value Standard Deviation p1 Survival factor, predators [1/year] -0.995895 0.0125269 (estimated) in [-Inf, Inf] p2 Death factor, predators [1/year] 1.00441 0.0129368 (estimated) in [-Inf, Inf] p3 Survival factor, preys [1/year] 1.01234 0.0135413 (estimated) in [-Inf, Inf] p4 Death factor, preys [1/year] 1.01909 0.0121026 (estimated) in [-Inf, Inf] p5 Crowding factor, preys [1/year] 0.103244 0.0039285 (estimated) in [-Inf, Inf] Status: Termination condition: Change in cost was less than the specified tolerance.. Number of iterations: 6, Number of function evaluations: 7 Estimated using Solver: ode45; Search: lsqnonlin on time domain data "ze". Fit to estimation data: [91.21;92.07]% FPE: 8.613e-06, MSE: 0.005713 More information in model's "Report" property. Model Properties

计算一个提前10步的预测输出来验证模型。

yhat4 =预测(sysGrey泽10);

预测模型的输出比估计数据高出100步，并计算输出标准差。

[yf4, x04 sysf4 ysd4] =预测(sysGrey,泽,100);

绘制测量的、预测的和预测的输出。

t = yf4.SamplingInstants;次要情节(1、2、1);情节(t0 z.y (: 1),.．.te, yhat4.y (: 1),.．.t, yf4.y (: 1),“米”，.．.t, yf4.y (: 1) + ysd4 (: 1),“k——”，.．.t, yf4.y (: 1) -ysd4 (: 1),“k——”)包含(的时间(年)）;ylabel (“捕食者数量(千)”）;标题(“灰盒模型的预测结果”)次要情节(1、2、2);情节(t0 z.y (:, 2),.．.te, yhat4.y (:, 2),.．.t, yf4.y (:, 2),“米”，.．.t, yf4.y (:, 2) + ysd4 (:, 2),“k——”，.．.t, yf4.y (:, 2) -ysd4 (:, 2),“k——”)传说(“测量”，“预测”，“预测”，“预测不确定性(1sd)”)包含(的时间(年)）;ylabel (“猎物种群(千)”）;

图中包含2个轴对象。标题为“灰盒模型预测结果”的坐标轴对象1包含5个类型为line的对象。Axes对象2包含5个line类型的对象。这些对象表示测量、预测、预测、预测不确定度(1 sd)。

结果表明，采用非线性灰盒模型进行预测具有较好的预测效果和较低的预测输出不确定性。

比较预测性能

比较从识别的模型得到的预测响应，sysARMA，sysNLARX,sysGrey．前两个是离散时间模型和sysGrey是一个连续时间模型。

clf情节(z, yf1, yf2 yf3, yf4)传说({“测量”，“ARMA”，“多项式AR”，“全科医生”，“方框”})标题(“预测反应”23) xlim ([7])

图中包含2个轴对象。标题为y1的axis对象1包含5个类型为line的对象。这些对象表示度量，ARMA，多项式AR, GP，灰盒。标题为y2的Axes对象2包含5个类型为line的对象。这些对象表示度量，ARMA，多项式AR, GP，灰盒。

非线性ARX模型的预测效果优于线性模型的预测效果。在非线性灰盒模型中加入动力学知识，进一步提高了模型的可靠性，从而提高了预测精度。

注意，灰盒建模中使用的方程与非线性ARX模型中使用的多项式回归量密切相关。如果你用一阶差分近似控制方程中的导数，你会得到类似于:

$y_{1} （ t ）＝ {年代}_{1} ＊ y_{1} （ t - 1 ） + {年代}_{2} ＊ y_{1} （ t - 1 ）＊ y_{2} （ t - 1 ）$

$y_{2} （ t ）＝ {年代}_{3.} ＊ y_{2} （ t - 1 ） - {年代}_{4} ＊ y_{1} （ t - 1 ）＊ y_{2} （ t - 1 ） - {年代}_{5} ＊ y_{2} （ t - 1 ）^{2}$

在那里, ${年代}_{1} ，．．．， {年代}_{5}$ 是否有些参数与原参数有关 $p_{1} ，．．．， p_{5}$ 以及用于差分的样本时间。这些方程表明，在构建非线性ARX模型时，第一个输出只需要2个回归量(y1(t-1)和*y1(t-1)*y2(t-1))，第二个输出需要3个回归量。即使在缺乏这些先验知识的情况下，使用多项式回归量的线性回归模型结构在实践中仍然是一个流行的选择。

利用非线性灰盒模型预测200年的数值。

[yf5，~，~，ysd5] =预报(sysGrey, ze, 2000);

绘制数据的后一部分(显示1sd的不确定度)

t = yf5.SamplingInstants;N = 700:2000;次要情节(1、2、1);情节(t (N) yf5.y (N, 1),“米”，.．.t (N), yf5.y (N - 1) + ysd5 (N, 1),“k——”，.．.t (N), yf5.y (N, 1) -ysd5 (N, 1),“k——”)包含(的时间(年)）;ylabel (“捕食者数量(千)”）;标题(“200年后的人口预测”) ax = gca;斧子。YLim = [0.8 1];斧子。XLim = [82 212];次要情节(1、2、2);情节(t (N) yf5.y (N, 2),“米”，.．.t (N), yf5.y (N, 2) + ysd5 (N, 2),“k——”，.．.t (N), yf5.y (N, 2) -ysd5 (N, 2),“k——”)传说(人口预测的，“预测不确定性(1sd)”)包含(的时间(年)）;ylabel (“猎物种群(千)”）;甘氨胆酸ax =;斧子。YLim = [0.9 1.1];斧子。XLim = [82 212];

图中包含2个轴对象。标题为“200年后预测人口”的坐标轴对象1包含3个类型为line的对象。Axes对象2包含3个line类型的对象。这些对象代表预测人口，预测不确定性(1 sd)。

图中显示，捕食者的数量预计将达到大约890的稳定状态，而被捕食者的数量预计将达到990。

预测多元时间序列

数据描述

线性模型的估计

估计非线性黑箱模型

多项式AR模型

估计一个高斯过程模型

非线性灰盒模型的估计

比较预测性能

相关的话题