预测多变量时间序列

打开直播脚本

此示例显示了如何在猎物拥挤场景中执行从捕食者和猎物群体测量的数据的多变量时间序列预测。使用线性和非线性时间序列模型建模捕食者 - 猎物人口变化动态。比较预测这些模型的性能。

数据描述

数据是一款由1-predator 1-prey群体组成的双变量时间序列（成千上万），每年收集10次。有关数据的更多信息，请参阅三种生态人口系统：MATLAB和CMEX文件建模时序。

加载时间序列数据。

加载predpreycrowdingdata.z = iddata（y，[]，0.1，'不时'那'年'那'tstart'，0）;

Z.是一个iddata.包含两个输出信号的对象，y1.和Y2.，分别指的是捕食者和猎物群体。这outputData.财产Z.包含人口数据作为201-BY-2矩阵，使得这一点z.outputdata（：，1）是捕食者人口和z.outputdata（:,2）是猎物人口。

绘制数据。

绘制（Z）标题（'捕食者 - 猎物人口数据'）ylabel（'人口（成千上万）'）

图包含2个轴。带标题Y1的轴1包含类型线的对象。该对象表示z。标题Y2的轴2包含类型线的对象。该对象表示z。

由于拥挤，数据表现出捕食者人口的下降。

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

ZE = Z（1：120）;

使用剩余的数据来搜索模型订单，并验证预测结果。

zv = z（121：结束）;

估计线性模型

将时间序列模拟为线性，自回归过程。线性模型可以使用多项式形式或使用命令（如）以多项式形式或状态空间形式创建AR.（仅适用于标量序列），ARX.那armax.那n4sid.和SSEST.。由于线性模型不会捕获数据偏移（非零条件平均值），因此首先拒绝数据。

[zed，tze] = detrend（ze，0）;[zvd，tzv] = detrend（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）;na2 = selstruc（v2,0）;

这ARXSTRUC命令旨在分别建议订单7和8的自动增加模型。

使用这些模型订单来估计一个多方方arma模型，其中横向术语已被任意选择。

na = [na1 na1-1;na2-1 na2];nc = [na1;na2];sysarma = ARMAX（ZED，[NA NC]）

sysarma =离散时间arma模型：输出“y1”模型：a（z）y_1（t）=  -  a_i（z）y_i（t）+ c（z）e_1（z）a（z）= 1  -  0.885Z ^ -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_2（z）= 0.3433 z ^ -1-  0.2802 Z ^ -2  -  0.04949 Z ^ -3 + 0.018 Z ^ -4  -  0.02683 Z ^ -5  -  0.2416 Z ^ -6 C（Z）= 1  -  0.4534 Z ^ -1  -  0.4127 Z ^ -2 + 0.7874Z ^ -3 + 0.298 Z ^ -4  -  0.8684 Z ^ -5 + 0.6106 Z ^ -6 + 0.3616 Z ^ -7输出“y2”的型号：a（z）y_2（t）=  -  a_i（z）y_i（T）+ C（Z）E_2（T）A（Z）= 1  -  0.5826 Z ^ -1  -  0.4688 Z ^ -2  -  0.5949 Z ^ -3  -  0.0547 Z ^ -4 + 0.5062 Z ^ -5 + 0.4024Z ^ -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.2217Z ^ -6  -  0.05765 z ^ -7  -  0.1824 z ^ -8采样时间：0.1年参数化：多项式订单：na = [7 6; 7 8] nc = [7; 8]自由系数的数量：43使用“polydata“，”getpvec“，”getcov“参数及其不正常散流。状态：在时域数据“zed”时估计使用ARMAX。适合估计数据：[89.85; 90.97]％（预测焦点）FPE：3.814E-05，MSE：0.007533

计算10步前（1年）预测输出以在估计数据的时间跨度上验证模型。由于数据被解除为估计，因此您需要指定有意义预测的偏移量。

predopt = predictoptions（'Outporfoffset'，tze.outputoffset'）;Yhat1 =预测（Sysarma，ZE，10，PredOpt）;

这预测命令预测测量数据时间跨度的响应，并且是用于验证估计模型的质量的工具。时间的响应T.在时间t = 0，...，T-10使用测量值计算。

绘制预测的响应和测量数据。

绘制（ZE，YHAT1）标题（与测量数据相比，'10 -Step预测响应'）

$图包含2个轴。具有标题Y1的轴1包含2个类型的线。这些对象代表ze，ze \ _preedette。带标题Y2的轴2包含2个类型的2个物体。这些对象代表ze，ze \ _preedette。$

注意，可以使用测量数据来实现预测响应和绘制响应的生成，可以自动使用相比命令。

compareopt = compareOptions（'Outporfoffset'，tze.outputoffset'）;比较（ze，sysarma，10，compareopt）

图包含2个轴。轴1包含2个类型的2个物体。这些物体代表ZE（Y1），SYSARMA：75.49％。轴2包含2个类型的2个物体。这些物体代表ZE（Y2），SYSARMA：76.21％。

使用的曲线相比还显示了赋予百分比的归一化均方根（Nrmse）效率的良好。

验证数据后，预测模型的输出sysarma.超出估计数据的100步（10年），并计算输出标准偏差。

Forecastopt =预测选项（'Outporfoffset'，tze.outputoffset'）;[YF1，X01，SYSF1，YSD1] =预测（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），'M'那......t，yf1.y（：，1）+ ysd1（：1），'k-'那......t，yf1.y（：，1）-ysd1（：，1），'k-'）Xlabel（'时间（年）'）;ylabel（'捕食者人口，成千上万）;子图（1,2,2）;绘图（T0，Z.Y（:,2），......Te，Yhat1.y（：2），......T，YF1.Y（:,2），'M'那......T，YF1.Y（:,2）+ YSD1（:,2），'k-'那......t，yf1.y（:,2）-ysd1（:,2），'k-'）％使图更大。图= GCF;p = fig.position;图。页位= [P（1），P（2）-P（4）* 0.2，P（3）* 1.4，P（4）* 1.2];Xlabel（'时间（年）'）;ylabel（'猎物人口，数千人'）;传奇（{'衡量'那'预料到的'那'预测'那'预测不确定性（1 SD）'}，......'地点'那'最好'）

图包含2个轴。轴1包含5个类型的线。轴2包含5型对象的类型。这些对象表示测量，预测，预测，预测不确定性（1 SD）。

该图表明，使用线性ARMA模型（额外处理偏移处理）的预测有所作为，结果表现出高度不确定性，与实际群体相比，在12 - 20年的时间范围内。这表明人口变更动态可能是非线性的。

估计非线性黑匣子型号

将非线性黑盒模型适合估计数据。您不需要先前了解有关估算数据的方程式的知识。将估计线性反回形式的非线性ARX模型。

创建具有2个输出和没有输入的非线性ARX模型。

sysnlarx = idnlarx（[1 1; 1 1]，[]，'ts'，0.1，'不时'那'年'）;

Sysnlarx.是一个使用没有非线性功能的第一阶非线性ARX模型;它使用其一阶回归量的加权和预测模型响应。

getreg（sysnlarx）

ans =.2x1细胞{'y1（t-1）'} {'y2（t-1）'}

要引入非线性功能，请向模型添加多项式回归。

创建回归给电源2，并包括交叉术语（上面列出的标准回归器的产品）。下载188bet金宝搏将这些回归作为自定义回归器添加到模型中。

r = polyreg（sysnlarx，'最大功率'2，'克罗斯特勒'那'在'）;sysnlarx.customreogressors = r;getreg（sysnlarx）

ans =.5x1细胞{'y1（t-1）'} {'y2（t-1）'} {'y1（t-1）。^ 2'} {'y1（t-1）。* y2（t-1）'} {'y2（t-1）。^ 2'}

使用估计数据估算模型的系数（回归权重和偏移量），泽。

sysnlarx = nlarx（ze，sysnlarx）

sysnlarx =非线性多变量时间序列模型，具有2个输出输出：y1，y2回归量：1。变量中的线性回归y1，y2 2.自定义回归线：y1（t-1）。^ 2 3.自定义回归：Y1（T-1）。* y2（t-1）4。自定义回归线：y2（t-1）。^ 2所有回归器的列表所有模型输出都是线性的回归器。采样时间：0.1年状态：使用nlarx估计时域数据“ze”。适合估计数据：[88.34; 88.91]％（预测焦点）FPE：3.265E-05，MSE：0.01048

计算10步前预测的输出以验证模型。

Yhat2 =预测（Sysnlarx，Ze，10）;

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

YF2 =预测（Sysnlarx，ZE，100）;

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

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

t = yf2.samplinginstants;子图（1,2,1）;绘图（T0，Z.Y（：1），......Te，Yhat2.y（：1），......t，yf2.y（：1），'M'）Xlabel（'时间（年）'）;ylabel（'捕食者人口（成千上万）'）;子图（1,2,2）;绘图（T0，Z.Y（:,2），......Te，Yhat2.y（:,2），......T，YF2.Y（：2），'M'） 传奇（'衡量'那'预料到的'那'预测'）Xlabel（'时间（年）'）;ylabel（'猎物人口（成千上万）'）;

图包含2个轴。轴1包含3个类型的线。轴2包含3个类型的线。这些对象代表测量，预测预测。

该图表明，使用非线性ARX模型的预测比使用线性模型提供了更好的预测结果。非线性黑匣子盒建模不需要先前了解管理数据的方程。

请注意，要减少回归数量的数量，可以使用主成分分析选择（转换）变量的最佳子集（参见PCA.）或特征选择（见序列）在统计和机器学习工具箱™中。

如果您已经先前了解系统动态，则可以使用非线性灰度框模型拟合估计数据。

估计非线性灰度盒模型

关于动态性质的知识可以帮助提高模型质量，从而有助于提高预测精度。对于捕食者 - 猎物动态，捕食者的变化（y1.）猎物（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和CMEX文件建模时序。

基于这些方程构造非线性灰度盒模型。

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

filename =.'predpley2_m';

指定模型订单（输出，输入和状态）

订单= [2 0 2];

指定参数的初始值 ${P.}_{1}$ 那 ${P.}_{2}$ 那 ${P.}_{3.}$ 那 ${P.}_{4.}$ ，和 ${P.}_{5.}$ ，并表明估计所有参数。请注意，在估计黑匣子型号时，要求对参数指定初始猜测的要求sysarma.和Sysnlarx.。

参数= struct（'名称'，{'生存因子，捕食者''死亡因素，捕食者'......“生存因子，拍摄'死亡因素，拍摄......'拥挤因素，拍摄}，......'单元'，{'1年''1年''1年''1年''1年'}，......'价值'，{ -  1.1 0.9 1.1 0.9 0.2}，......'最低限度'，{ -  inf -inf -inf -inf -inf}，......'最大'，{inf inf inf inf}，......'固定的'，{false value false false};

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

initialstates = struct（'名称'，{'捕食者人口''猎物人口'}，......'单元'，{'大小（成千上）''大小（成千上）'}，......'价值'，{1.8 1.8}，......'最低限度'，{0 0}，......'最大'，{inf inf}，......'固定的'，{false false}）;

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

ts = 0;

创建具有指定结构，参数和状态的非线性灰度框模型。

sysgrey = idnlgrey（文件名，订单，参数，initialstates，ts，'不时'那'年'）;

估计模型参数。

sysgrey = nlgreyest（ze，sysgrey）;目前（SYSGEY）

sysgrey ='predprey2_m'（matlab文件）定义的连续时间非线性灰盒模型：dx / dt = f（t，x（t），p1，...，p5）y（t）= h（t，X（t），p1，...，p5）+ e（t），0输入，2个状态，2个输出和5个免费参数（为5）。状态：初始值x（1）捕食者群（t）[size（thou] xinit @ exp1 2.01325（估计）在[0，inf] x（2）猎物群（t）[size（thoul.c] xinit中@ Exp1 1.99687（估计）在[0，INF]输出：Y（1）Y1（T）Y（2）Y2（T）参数：值标准偏差P1存活因子，捕食者[1 /年] -0.995895 0.0125269（估计）在[-inf，inf] p2死亡因子，捕食者[1 /年] 1.00441 0.0129368（估计）在[-inf，inf] p3存活因子，Preys [1 /年] 1.01234 0.0135413（估计）在[-inf，inf] p4死亡因子，射精[1 /年] 1.01909 0.0121026（估计）在[-inf，inf] p5挤压系数，术前[1 /年] 0.103244 0.0039285（估计）在[-inf，inf]状态：终止条件：成本的变化小于指定的公差。迭代数量：6，功能评估数：7使用求解器估计：oDE45;搜索：LSQNONLIN在时域数据“ZE”。适合估计数据：[91.21; 92.07]％FPE：8.613E-06，MSE：0.005713 Model“报告”属性中的更多信息。

计算10步前预测的输出以验证模型。

Yhat3 =预测（Sysgrey，Ze，10）;

预测型号100步的输出超出估计数据，并计算输出标准偏差。

[YF3，X03，SYSF3，YSD3] =预测（SYSGREY，ZE，100）;

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

t = yf3.samplinginstants;子图（1,2,1）;绘图（T0，Z.Y（：1），......Te，Yhat3.y（：1），......T，YF3.Y（：1），'M'那......t，yf3.y（：，1）+ ysd3（：1），'k-'那......t，yf3.y（：，1）-ysd3（：1），'k-'）Xlabel（'时间（年）'）;ylabel（'捕食者人口（成千上万）'）;子图（1,2,2）;绘图（T0，Z.Y（:,2），......Te，Yhat3.y（:,2），......t，yf3.y（:,2），'M'那......t，yf3.y（:,2）+ ysd3（:,2），'k-'那......t，yf3.y（:,2）-ysd3（:,2），'k-'） 传奇（'衡量'那'预料到的'那'预测'那'预测不确定性（1 SD）'）Xlabel（'时间（年）'）;ylabel（'猎物人口（成千上万）'）;

图包含2个轴。轴1包含5个类型的线。轴2包含5型对象的类型。这些对象表示测量，预测，预测，预测不确定性（1 SD）。

该图表明，使用非线性灰度盒模型的预测得到了良好的预测结果和低预测输出不确定性。

比较预测性能

比较从已识别的模型获得的预测响应，sysarma.那Sysnlarx.，和sysgrey。前两个是离散时间模型sysgrey是一个连续时间模型。

CLF图（Z，YF1，YF2，YF3）传奇（{'衡量'那'linear arma'那'非线性AR'那'非线性灰盒'}） 标题（'预测回复'）

图包含2个轴。带标题Y1的轴1包含4个类型的线。这些物体表示测量，线性ARMA，非线性AR，非线性灰度盒。带标题Y2的轴2包含4型对象。这些物体表示测量，线性ARMA，非线性AR，非线性灰度盒。

与非线性ARX模型的预测比用线性模型预测结果比预测更好。包含非线性灰度箱模型中动力学的知识进一步提高了模型的可靠性，从而提高了预测精度。

请注意，灰度盒建模中使用的方程与非线性ARX模型使用的多项式回归密切相关。如果通过一阶差异近似于管理方程式的衍生物，则将获得类似于：的等式

$y_{1} （ T. ） = {S.}_{1} * y_{1} （ T. - 1 ） + {S.}_{2} * y_{1} （ T. - 1 ） * y_{2} （ T. - 1 ）$

$y_{2} （ T. ） = {S.}_{3.} * y_{2} （ T. - 1 ） - {S.}_{4.} * y_{1} （ T. - 1 ） * y_{2} （ T. - 1 ） - {S.}_{5.} * y_{2} （ T. - 1 ）^{2}$

在哪里， ${S.}_{1} 那。。。那 {S.}_{5.}$ 是与原始参数相关的一些参数 ${P.}_{1} 那。。。那 {P.}_{5.}$ 以及用于差异的采样时间。这些等式表明，在构造非线性ARX模型时，第一输出（Y1（T-1）和* Y1（T-1）* Y2（T-1））和3的第二输出仅需要2个回归量。即使在没有这样的先前知识，也采用多项式回归的线性反回的模型结构仍然是实践中的流行选择。

预测使用200年超过200年的非线性灰度框模型的值。

[YF4，〜，〜，YSD4] =预测（Sysgrey，Ze，2000）;

绘制数据的后半部分（显示1个SD不确定性）

t = yf4.samplinginstants;n = 700：2000;子图（1,2,1）;绘图（t（n），yf4.y（n，1），'M'那......t（n），yf4.y（n，1）+ ysd4（n，1），'k-'那......t（n），yf4.y（n，1）-ysd4（n，1），'k-'）Xlabel（'时间（年）'）;ylabel（'捕食者人口（成千上万）'）;斧头= GCA;ax.ylim = [0.8 1];ax.xlim = [82 212];子图（1,2,2）;绘图（t（n），yf4.y（n，2），'M'那......t（n），yf4.y（n，2）+ ysd4（n，2），'k-'那......t（n），yf4.y（n，2）-ysd4（n，2），'k-'） 传奇（'预测人口'那'预测不确定性（1 SD）'）Xlabel（'时间（年）'）;ylabel（'猎物人口（成千上万）'）;斧头= GCA;ax.ylim = [0.9 1.1];ax.xlim = [82 212];