多元时间序列预测
这个例子展示了如何执行多元时间序列预测的数据测量从捕食者和猎物拥挤的场景。利用线性和非线性时间序列模型对捕食者-食饵种群变化动态进行了建模。比较了这些模型的预测效果。
数据描述
数据是一个二元时间序列，由1-捕食者1-被捕食者(以千计)组成，每年采集10次，持续20年。有关数据的更多信息，请参见三个生态种群系统:MATLAB和C mex -文件时间序列建模．
加载时间序列数据。
负载预排序数据z = iddata (y, [], 0.1,“TimeUnit”，“年”，“Tstart”,0);
z是一个iddata对象包含两个输出信号，日元和y2，分别指捕食者和被捕食者的数量。的OutputData的属性z以201×2矩阵形式包含总体数据，因此z、 输出数据（：，1）是捕食者的种群和数量z、 输出数据（：，2）是猎物的数量。
绘制数据。
图（z）标题(“捕食者-猎物人口数据”)伊拉贝尔(‘人口（千）’）
数据显示，由于拥挤，捕食者数量下降。
将前半部分用作识别时间序列模型的估计数据。
ze=z（1:120）；
使用剩余的数据来搜索模型订单，并验证预测结果。
zv z =(121:结束);
估计线性模型
将时间序列建模为线性的、自回归的过程。线性模型可以以多项式形式或状态空间形式创建，使用命令如基于“增大化现实”技术(仅对标量时间序列)，阿克斯，阿玛克斯，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);钠= selstruc (V2, 0);
的arxstrucCommand建议分别使用7和8阶的自回归模型。
在任意选择交叉项的情况下，利用这些模型阶数估计多元方差ARMA模型。
na=[na1na1-1；na2-1na2]；nc=[na1；na2]；sysARMA=armax（zed，[na-nc]）
sysARMA =离散时间ARMA模型:输出“y1”模型:(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 ^ ^ 5 + 4 + 0.06555 z参数化:多项式阶数: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 (“OutputOffset”Tze.OutputOffset ');yhat1 =预测(sysARMA,ze,10, predOpt);
的预测Command预测在测量数据的时间跨度内的响应，是验证估计模型质量的工具。响应时间t是用t = 0时的实测值计算的，…, t-10。
绘制预测响应和测量数据。
情节（ze，yhat1）标题(“与实测数据相比的10步预测响应”）
注意，预测响应的生成和测量数据的绘图，可以自动使用比较命令。
compareOpt = compareOptions (“OutputOffset”Tze.OutputOffset ');比较(泽sysARMA 10 compareOpt)
使用比较也以百分比形式显示了归一化均方根(NRMSE)的拟合优度度量。
验证数据后，预测模型的输出西萨玛超出估计数据100步(10年)，计算输出标准差。
forecastOpt=forecastOptions(“OutputOffset”Tze.OutputOffset ');[yf1,x01,sysf1,ysd1] = forecast(sysARMA, ze, 100, forecastOpt);
yf1预测的响应是否返回为iddata对象yf1。OutputData包含预测值。
sysf1系统是否类似于西萨玛而是处于状态空间的形式。模拟sysf1使用sim卡在初始条件下，x01，复制预测的响应，yf1．
ysd1是标准差的矩阵。它测量的不确定性是预测由于附加干扰在数据的影响(如测量sysARMA。NoiseVariance)，参数不确定性(由getcov (sysARMA))，以及将过去的数据映射到预测所需的初始条件过程中的不确定性(见data2state)．
绘制模型的测量、预测和预测输出西萨玛．
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——”)xlabel(的时间(年))；伊莱贝尔(“捕食者种群，以千为单位”)子地块（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);包含(的时间(年))；伊莱贝尔(“猎物数量，以千计”）;传奇({“测量”，“预测”，“预测”，“预测的不确定性(1 sd)”}，...“位置”，“最佳”）
图中显示，使用线性ARMA模型(添加了补偿处理)进行预测的效果还不错，与12-20年时间跨度内的实际人口相比，结果显示出较高的不确定性。这表明种群变化动态可能是非线性的。
估计一个非线性黑箱模型
拟合一个非线性黑箱模型到估计数据。您不需要关于控制估计数据的方程的先验知识。一个线性回归形式的非线性ARX模型将被估计。
创建一个具有2个输出而无输入的非线性ARX模型。
sysNLARX = idnlarx([1 1;1 1]，[]，“t”, 0.1,“TimeUnit”，“年”）;
sysNLARX为不使用非线性函数的一阶非线性ARX模型;它预测模型响应使用其一阶回归量的加权和。
getreg (sysNLARX)
输出1:y1(t-1) y2(t-1) y2(t-1)
要引入非线性函数，请在模型中添加多项式回归。
创建幂次为2的回归量，并包括交叉项(上面列出的标准回归量的乘积)。下载188bet金宝搏将这些回归器作为自定义回归器添加到模型中。
R = polyreg (sysNLARX“MaxPower”，2,“CrossTerm”，“上”); sysNLARX.CustomRegressors=R；getreg（sysNLARX）
回归器：用于输出1:y1（t-1）y2（t-1）y1（t-1）。^2 y1（t-1）。*y2（t-1）y2（t-1）。^2用于输出2:y1（t-1）y2（t-1）y1（t-1）。^2 y1（t-1）。*y2（t-1）y2（t-1）。^2
使用估计数据估计模型的系数（回归权重和偏移量），泽．
sysNLARX sysNLARX = nlarx(泽)
sysNLARX=具有2个输出和0个输入的非线性ARX模型输入：输出：y1，y2对应于顺序的标准回归器：na=[11；11]nb=[]nk=[]自定义回归器：用于输出1:y1（t-1）。^2 y1（t-1）。*y2（t-1）y2（t-1）。^2用于输出2:y1（t-1）。^2 y1（t-1）。*y2（t-1）y2（t-1）.2非线性回归：对于输出1:none对于输出2:none模型输出在其回归中是线性的。样本时间：0.1年状态：使用NLARX对时域数据“ze”进行估计。适合估计数据：[88.34；88.91]（预测焦点）FPE:3.265e-05，MSE:0.01048
计算提前10步的预测输出，以验证模型。
yhat2 =预测(sysNLARX泽10);
预测超出估计数据100步的模型输出。
yf2=预测（sysNLARX，ze，100）；
对于非线性ARX模型，未计算预测响应的标准差。这个数据是不可用的，因为在估计这些模型时没有计算参数协方差信息。
绘制测量的、预测的和预测的输出。
t = yf2.SamplingInstants;次要情节(1、2、1);情节(t0 z.y (: 1),...te, yhat2.y (: 1),...t, yf2.y (: 1),“米”)xlabel(的时间(年))；伊莱贝尔(“捕食者数量(千)”)子地块（1,2,2）；地块（t0，z.y.：，2），...te, yhat2.y (:, 2),...t, yf2.y (:, 2),“米”)传说(“测量”，“预测”，“预测”)xlabel(的时间(年))；伊莱贝尔(“猎物种群(千)”）;
从图中可以看出，使用非线性ARX模型的预测结果比使用线性模型的预测结果要好。非线性黑箱建模不需要关于控制数据的方程的先验知识。
请注意，为了减少回归器的数量，可以使用主成分分析选择（转换）变量的最佳子集（参见主成分分析)或特征选择(见sequentialfs)在统计和机器学习工具箱中™.
如果您有系统动力学的先验知识，您可以使用非线性灰箱模型拟合估计数据。
估计一个非线性灰箱模型
了解动态的本质有助于提高模型的质量，从而提高预测的准确性。对于捕食者-被捕食者动力学，捕食者(日元)及猎物(y2)人口可表示为:
有关方程式的详细信息，请参见三个生态种群系统:MATLAB和C mex -文件时间序列建模．
基于这些方程建立了一个非线性灰箱模型。
为捕食者-猎物系统指定一个描述模型结构的文件。该文件将状态导数和模型输出指定为时间、状态、输入和模型参数的函数。选择两个输出(捕食者和食饵种群)作为状态，导出动力学的非线性状态空间描述。
文件名=“predprey2_m”；
指定模型顺序（输出、输入和状态的数量）
订单= [2 0 2];
指定参数的初始值
，
，
，
和
，表示所有参数均需估计。注意，在估计黑盒模型时，不存在为参数指定初始猜测的需求西萨玛和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 Inf}，...“固定”，{假假假假});
类似地，指定模型的初始状态，并指出两个初始状态都要估计。
InitialStates=struct(“名字”,{“捕食者种群”“猎物人口”}，...“单位”,{‘大小（千）’‘大小（千）’}，...“价值”{1.8 - 1.8},...“最低限度”, {0 0},...“最大”{正正无穷},...“固定”,{假假});
将模型指定为连续时间系统。
t = 0;
创建具有指定结构、参数和状态的非线性灰箱模型。
sysGrey = idnlgrey(文件名、秩序、参数、InitialStates Ts,“TimeUnit”，“年”）;
估计模型参数。
Sysgree=nlgreyest（ze，Sysgree）；present（Sysgree）
sysGrey = '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.
计算提前10步的预测输出，以验证模型。
yhat3=预测（系统灰色，ze，10）；
预测超出估计数据100步的模型输出，并计算输出标准偏差。
[yf3, x03 sysf3 ysd3] =预测(sysGrey,泽,100);
绘制测量的、预测的和预测的输出。
t = yf3.SamplingInstants;次要情节(1、2、1);情节(t0 z.y (: 1),...te, yhat3.y (: 1),...t、 yf3.y（：，1），“米”，...t, yf3.y (: 1) + ysd3 (: 1),“k——”，...t, yf3.y (: 1) -ysd3 (: 1),“k——”)xlabel(的时间(年))；伊莱贝尔(“捕食者数量(千)”)子地块（1,2,2）；地块（t0，z.y.：，2），...te, yhat3.y (:, 2),...t, yf3.y (:, 2),“米”，...t, yf3.y (:, 2) + ysd3 (:, 2),“k——”，...t, yf3.y (:, 2) -ysd3 (:, 2),“k——”)传说(“测量”，“预测”，“预测”，“预测的不确定性(1 sd)”)xlabel(的时间(年))；伊莱贝尔(“猎物种群(千)”）;
结果表明，采用非线性灰箱模型进行预测具有较好的预测效果和较低的预测输出不确定性。
比较预测性能
比较从确定的模型中获得的预测响应，西萨玛，sysNLARX和sysGrey.前两个是离散时间模型和sysGrey是一个连续时间模型。
clf情节(z, yf1, yf2 yf3)传说({“测量”，的线性自回归滑动平均，“非线性AR”，“非线性灰盒”})标题(“预测反应”）
使用非线性ARX模型进行预测比使用线性模型进行预测效果更好。在非线性灰箱模型中加入动力学知识，进一步提高了模型的可靠性，从而提高了预测精度。
注意，灰箱模型中使用的方程与非线性ARX模型中使用的多项式回归器密切相关。如果你用一阶差分近似控制方程的导数，你会得到类似的方程:
哪里 有些参数与原始参数有关吗 以及用于差分的采样时间。这些方程表明，在构建非线性ARX模型时，第一个输出只需要2个回归器(y1(t-1)和*y1(t-1)*y2(t-1))，第二个输出只需要3个回归器。即使在没有这些先验知识的情况下，线性回归或使用多项式回归的模型结构在实践中仍然是一个流行的选择。
使用非线性灰箱模型预测200年来的数值。
[yf4，~，~，ysd4] = forecast(sysGrey, ze, 2000);
绘制数据的后一部分(显示1标准差的不确定性)
t = yf4.SamplingInstants;N = 700:2000;次要情节(1、2、1);情节(t (N) yf4.y (N, 1),“米”，...t (N), yf4.y (N - 1) + ysd4 (N, 1),“k——”，...t (N), yf4.y (N, 1) -ysd4 (N, 1),“k——”)xlabel(的时间(年))；伊莱贝尔(“捕食者数量(千)”); ax=gca；ax.YLim=[0.81]；ax.XLim=[82 212]；子批次（1,2,2）；图（t（N），yf4.y（N，2），“米”，...t (N), yf4.y (N, 2) + ysd4 (N, 2),“k——”，...t (N), yf4.y (N, 2) -ysd4 (N, 2),“k——”)传说(“预测人口”，“预测的不确定性(1 sd)”)xlabel(的时间(年))；伊莱贝尔(“猎物种群(千)”）;甘氨胆酸ax =;斧子。YLim = [0.9 1.1];斧子。XLim = [82 212];
该图显示，捕食者种群预计将达到约890的稳定状态，猎物种群预计将达到990。
