多元时间序列预测

打开生活的脚本

这个例子展示了如何执行多元时间序列预测的数据测量从猎物的捕食者和猎物种群拥挤场景。捕食者-猎物人口变化动力学的建模使用线性和非线性时间序列模型。这些模型的预测性能进行比较。

数据描述

数据组成的二元时间序列1-predator 1-prey数量(千)收集20年来每年的10倍。关于数据的更多信息,请参阅三个人口生态系统: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轴对象。坐标轴对象1标题y1包含一个类型的对象。这个对象表示z轴对象2标题y2包含一个类型的对象。该对象代表z。

数据展示捕食者种群下降由于拥挤。

利用上半年估计数据识别时间序列模型。

泽= z (1:120);

用剩下的订单数据搜索模型,并验证预测结果。

zv z =(121:结束);

估计线性模型

时间序列的线性模型,自回归过程。线性模型可以在多项式形式使用命令或状态方程形式等基于“增大化现实”技术(仅供标量时间序列),arx,armax,n4sid和党卫军。由于线性模型不捕获数据偏移量(零条件的意思),第一次去趋势数据。

[zed,老子]=去趋势(泽,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);

的arxstruc命令显示订单7和8的自回归模型,分别。

使用这些模型的订单估计multi-variance 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]数控=(7、8)数量的免费系数:43使用“polydata”、“getpvec”、“getcov”参数及其不确定性。状态:估计ARMAX使用时域数据“z”。适合估算数据:[89.85,90.97]%(预测聚焦)消防工程:3.814 e-05, MSE: 0.007533

计算10-step-ahead(1年)输出验证模型预测的时间跨度估计数据。由于数据是去趋势估计,您需要指定偏移量的有意义的预测。

predOpt = predictOptions (“OutputOffset”Tze.OutputOffset ');yhat1 =预测(sysARMA泽10,predOpt);

的预测命令响应预测的时间跨度测量数据和验证的工具的质量估计模型。响应时间t有时使用测量值计算t = 0时,…,t-10。

预测的响应和测量数据的阴谋。

情节(泽,yhat1)标题(所述的预测响应相比,测量数据的)

$图包含2轴对象。坐标轴对象1标题y1包含2线类型的对象。这些对象代表泽泽\ _Predicted。坐标轴对象2与标题y2包含2线类型的对象。这些对象代表泽泽\ _Predicted。$

注意,生成预测响应和策划的测量数据,可以自动使用比较命令。

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

图包含2轴对象。坐标轴对象1包含2线类型的对象。这些对象代表泽(日元),sysARMA: 75.49%。坐标轴对象包含2线类型的对象。这些对象代表泽(y2), sysARMA: 76.21%。

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

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

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

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 (人口的猎物,数以千计的);传奇({“测量”,“预测”,“预测”,“预测不确定性(1 sd)”},…“位置”,“最佳”)

图包含2轴对象。坐标轴对象1包含5线类型的对象。坐标轴对象2包含5线类型的对象。这些对象代表测量,预测,预测,预报不确定性(1 sd)。

情节表明,使用一个线性ARMA模型预测(偏移量添加处理)工作,结果显示高不确定性比实际人口12-20年时间跨度。这表明人口变化可能是非线性动力学方程。

估计非线性黑箱模型

在本节中,我们扩展黑箱辨识方法将非线性元素添加到前一节的模型。这是使用非线性ARX建模框架。非线性ARX模型在概念上类似于线性的(如sysARMA以上),他们计算出预测输出使用一个两步的过程:

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

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

多项式AR模型

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

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)”

引入一个非线性函数,多项式解释变量添加到模型中。

创建订单2的解释变量,包括交叉项(上面列出的产品的标准解释变量)。下载188bet金宝搏这些解释变量添加到模型中。

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 = 2输出输出非线性多变量时间序列模型:y1, y2解释:1。线性解释变量y1, y2 2。订单2解释变量y1, y2的所有解释变量的列表输出功能:输出1:没有输出2:没有一个样品时间:0.1年状态:估计使用NLARX时域数据“泽”。适合估算数据:[88.34,88.91]%(预测聚焦)消防工程:3.265 e-05, MSE: 0.01048

计算一个10-step-ahead预测输出验证模型。

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轴对象。坐标轴对象1标题预测结果ARX模型与多项式解释包含3线类型的对象。坐标轴对象2包含3线类型的对象。这些对象代表测量,预测,预测。

情节表明,使用非线性ARX模型预测了比使用线性模型更好的预测结果。非线性黑箱建模不需要先验知识的方程管理数据。

注意,减少解释变量的数量,你可以选择一个最优的子集(转换)变量使用主成分分析(参见主成分分析)或特征选择(见sequentialfs在统计和机器学习的工具箱™。

估计高斯过程模型

该模型模型sysNLARX使用一个线性映射函数;中包含的非线性是解释变量。更灵活的方法是选择一个简单的映射函数,如高斯过程(GP)的函数。创建一个只使用线性的非线性ARX模型解释变量,但添加一个高斯过程函数的解释变量映射到输出。

%复制sysNLARX所以我们没有再次创建模型结构sysGP = sysNLARX;%去除多项式回归量sysGP.Regressors (2) = [];%替代线性映射函数的高斯过程(GP)的功能。%的医生使用一个平方指数默认内核。GP = idGaussianProcess;sysGP。OutputFcn =全科医生;disp (sysGP)

非线性多变量时间序列模型2输出输出:y1, y2解释变量:线性解释变量y1, y2输出功能:输出1:高斯过程使用SquaredExponential内核函数。输出2:高斯过程使用SquaredExponential内核函数。样品时间:0.1年状态:估计用NLARX预测的焦点

您现在拥有了一个模板模型sysGP基于线性解释变量和GP映射函数。它的参数是两个系数的平方指数内核。使用nlarx估计这些参数。

sysGP sysGP = nlarx(泽)

sysGP = 2输出输出非线性多变量时间序列模型:y1, y2解释变量:线性解释变量y1, y2的所有解释变量的列表输出功能:输出1:高斯过程使用SquaredExponential内核函数。输出2:高斯过程使用SquaredExponential内核函数。样品时间:0.1年状态:估计使用NLARX时域数据“泽”。适合估算数据:[88.07,88.84]%(预测聚焦)消防工程:3.369 e-05, MSE: 0.01083

getpvec (sysGP)

ans =10×1-0.6671 0.0196 0.9558 0.9445 0.0261 -0.0221 -0.5752 0.8479 1.1883 0.0447

和之前一样,计算预测和预测反应和阴谋的结果。

计算一个10-step-ahead预测输出验证模型。

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线类型的对象。坐标轴对象2包含3线类型的对象。这些对象代表测量,预测,预测。

结果表明改进的线性模型的准确性。使用基于GP模型的一个优势是,他们使用很少的参数,而不是其他形式,如小波网络(idWaveletNetwork)或乙状结肠网络(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”;

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

订单= (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},…“最低”,{负负负负无穷},…“最大”,{正正正正正},…“固定”,{假假假假假});

同样,指定模型的初始状态,表明两个初始状态估计。

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

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

t = 0;

创建一个指定的非线性灰色矩形模型结构、参数和状态。

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

估计模型参数。

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

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-step-ahead预测输出验证模型。

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——”)传说(“测量”,“预测”,“预测”,“预测不确定性(1 sd)”)包含(的时间(年));ylabel (“猎物种群(千)”);

图包含2轴对象。坐标轴对象1标题预测结果灰色矩形模型包含5线类型的对象。坐标轴对象2包含5线类型的对象。这些对象代表测量,预测,预测,预报不确定性(1 sd)。

情节表明,使用非线性灰色矩形模型预测了良好的预测效果和较低的预测输出不确定性。

比较预测性能

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

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

图包含2轴对象。坐标轴对象1标题y1包含5线类型的对象。这些对象代表测量,ARMA,多项式AR, GP,灰色矩形。坐标轴对象2与标题y2包含5线类型的对象。这些对象代表测量,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}$ 以及样品的时间用于差分。这些方程表明,只有两个解释变量所需的第一个输出(y1 (t - 1)和* (t - 1) * y2 (t - 1))和3第二当构造非线性ARX模型的输出。即使在缺乏先验知识,linear-in-regressor模型结构采用多项式解释变量在实践中仍然是一个受欢迎的选择。

使用非线性灰色矩形模型预测的值超过200年了。

[yf5, ~, ~, ysd5] =预测(sysGrey,泽,2000);

图数据的后期(显示1 sd不确定性)

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 =;斧子。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——”)传说(人口预测的,“预测不确定性(1 sd)”)包含(的时间(年));ylabel (“猎物种群(千)”);甘氨胆酸ax =;斧子。YLim = (0.9 - 1.1);斧子。XLim = (82 - 212);