第6章+随机型时间序列预测法.ppt

1、第6章 随机型时间序列预测法,6.1 随机型时间序列模型 6.2 ARMA模型的相关分析 6.3 模型的识别 6.4 ARMA序列的参数估计 6.5 模型的检验与预报 6.6 ARIMA模型,利用随机型时间序列预测法建立预测模型的过程：,（1） 根据建模的目的和理论分析，确定模型的基本形式；,（2） 进行模型识别，即从一大类模型中选出一类试验模型；,（3） 将所选的模型应用于所取得的历史数据，求得模型的参数；,（4） 检验得到的模型是否合适。若合适，则可用于预测或控制；若不适合，则返回到第二阶段重新选择模型。,6.1 随机型时间序列模型,时间序列 取自某一个随机过程，则称：,一、平稳时间序列,

2、过程是平稳的随机过程的随机特征不随时间变化而变化,过程是非平稳的随机过程的随机特征随时间变化而变化,1、宽平稳时间序列的定义,设时间序列,，对于任意的t,k和m，满足：,则称 宽平稳。,2、ARMA模型三种基本形式： 自回归模型（AR：Auto-regressive） 移动平均模型（MA：Moving-Average） 自回归移动平均模型（ARMA） 求和自回归移动平均模型（ARI MA）,如果时间序列 满足,二、自回归模型AR(p),其中 是独立同分布的随机变量序列,且满足：,引入向后推移算子B：,则称时间序列 服从p阶自回归模型,AR(p),则称自回归模型是平稳的。,记,则AR模型：,称多

3、项式方程,为AR模型的特征方程,称特征方程p个根：,为AR模型的特征根,如果p个特征根都在单位圆外，即：,平稳性条件,如：AR（1）模型：,其特征方程,特征根,且有,平稳性条件,实际上,因此，对平稳的AR（1）模型，,可由过去各期的误差线性表示。,推广,对平稳的AR（p）模型，,可由过去各期的误差线性表示。,AR（2）模型平稳性条件：,参数 2， 1的取值变化分三种情形讨论： （1）当 12 + 4 2 = 0 时，有相等特征根。2, 1取值在图中的抛物线上，称为临界阻尼状态。,(2) 当 12 + 4 2 0 时，有不等特征根。2, 1的值位于过阻尼区（自相关函数呈指数衰减）,(3)当 12

4、 + 4 2 0 时，为共轭复根。 2, 1的值位于欠阻尼区（自相关函数呈正弦震荡衰减）。,平稳域：是一个三角形区域，见图阴影部分,例 考察四个模型的平稳性,平稳序列时序图,例 非平稳序列时序图,AR(1)模型平稳条件,特征根,平稳域,AR(2)模型,例 平稳性判别,如果时间序列 满足,三、移动平均模型MA(q),则称时间序列 服从q阶移动平均模型。,或者记为,称多项式方程,为MA模型的特征方程,称特征方程q个根：,为MA模型的特征根,如果 q个特征根都在单位圆外，,则称此模型是可逆的。,如：MA（1）模型,其特征方程：,特征根,可逆性条件,由,得,因此，对可逆的MA（1）模型，,可由过去各期

5、 的误差线性表示。,因此，AR（p）模型与 MA（q）模型是 互相对偶的两个模型,四、自回归移动平均模型ARMA(p,q),如果时间序列,满足：,则称时间序列 服从(p,q)阶自回归移动平均 模型ARMA(p,q) 。,或者记为：,q=0，模型即为AR(p)； p=0，模型即为MA(q)。,ARMA(p,q)模型特殊情况：,五、求和自回归移动平均模型ARI MA(p,d,q),用来处理非平稳序列，转换成平稳序列,齐次非平稳序列：通过一次或多次差分， 可转换成平稳序列。 差分的次数称为齐次化的阶数,定义差分：,引入差分算子,n阶差分为：,二阶差分为：,(p,d,q)阶求和自回归移动平均模型为,亦

6、即 是ARMA(p,q)序列,p,q分别是 序列的自回归和移动平均的阶数，d为求和阶数。,假设,前d个随即变量是均值为0，方差有限，且与 不相关，,这样序列 可由它的初值 及平稳序列 表达。,六、季节性模型,对含有季节性周期的时间序列，也可用季节差分的方法将之化成平稳序列。,如：对季度波动，可用4次差分,如：对月度波动，可用12次差分,鲍克斯季节模型,如取,则有,若随机干扰项也是与季节相关，则模型为,例如,描述了一个既有线性发展趋势，又含月度周期变动的随机型时间序列模型。,6.2 ARMA模型的相关分析,自相关分析法是进行时间序列分析的有效方 法，它简单易行, 较为直观，根据绘制的自 相关分析

7、图和偏自相关分析图，我们可以初 步地识别平稳序列的模型类型和模型阶数。 利用自相关分析法可以测定时间序列的随机性 和平稳性，以及时间序列的季节性。,一、自相关分析,（1）自相关函数的定义,滞后期为k的自协方差函数为：,则自相关函数为：,其中,当序列平稳时，自相关函数可写为：,（2）样本自相关函数,其中,样本自相关函数可以说明不同时期的数 据之间的相关程度，其取值范围在-1到 1之间，值越接近于1，说明时间序列的 自相关程度越高。,如果 满足AR(p）模型：则称其为AR(p)序列,二、 AR(p)序列的自相关函数,用 乘上式两端，再取均值，则有,两边同除以 ,得AR(p)自相关函数方程,算子形式

8、：,利用AR(p)模型参数求解自相关函数的方法：,若 的根 互不相同，则通解,其中，,为常系数，其求解方法为：,因为,故,若 有重根，不妨为 ，,其重数分别为,则方程通解为,不管 有无重根，总可证明存在正常数,亦即 随k的增加按指数形式衰减，呈“拖尾”状。,利用统计样本估计自相关函数，再求模型参数的方法：,将样本估计值 代入方程，,可求得模型参数 的估计值.,Yule-Walker方程,模型中白噪声 的方差的方法：,而,所以,故,常用AR模型自相关系数递推公式,AR(1)模型 AR(2)模型,例:考察如下AR模型的自相关图,自相关系数按复指数单调收敛到零,自相关系数呈现出“伪周期”性,自相关系

9、数不规则衰减,如果 满足MA(q) 序列时，即,三、 MA(q)序列的自相关函数,则,表明：当 与 的相距步数 时,与 不相关，即MA(q)序列的自相关函数 从 以后全部为0。称这一性质为“截尾”的。,反之，若一个平稳时间序列的自相关（自协方差）函数 截尾，那么它必定是 MA(q)序列,常用MA模型的自相关系数,MA(1)模型,MA(2)模型,MA模型的自相关系数截尾,MA模型的自相关系数截尾,四、 ARMA(p,q)序列的自相关函数,结论：,（1）当,若 的根 互不相同，则通解,它具有拖尾性。,（2）若已知,得关于 线性方程组,可解出自回归参数,再求移动平均参数,因,则,令,即 为MA(q)

10、序列，,其自协方差函数,可通过 其自协方差函数 求出：,进而可以求出移动平均参数,例:ARMA模型的相关函数拖尾性,模型ARMA(1,1): 样本自相关图,五、AR偏相关函数,最优系数,使下式达到最小,考虑用 对 作最小方差估计，即考虑回归方程,只需求偏导：,即有,求得此方程即得,称序列 为 的偏相关函数,对AR(p)序列 ，由 定义知，对任意,若取,则 达到最小值,由此可见AR(p)序列的偏相关函数在 后等于0，即截尾的，它是特有标志。,对于AR(1)过程： 当k = 1时， 11 0；当k 1时，kk = 0， 所以AR(1)过程的偏自相关函数特征是在k = 1出现峰值（11 = 1），然

11、后截尾。,对于AR(2)过程， 当k 2时，kk 0，；当k 2时，kk = 0 偏自相关函数在滞后期2以后有截尾特性。,对于AR(p)过程： 当k p时，kk 0；当k p时，kk = 0。 偏自相关函数在滞后期p以后有截尾特性，因此可用此特征识别AR(p)过程的阶数。,偏自相关函数由下式中的粗项组成： xt = 11 xt-1 + ut xt = 21 xt-1 + 22 xt-2 + ut xt = k 1 xt-1 + k 2 xt-2 + + kk xt-k + ut,因偏自相关函数中每一个回归系数 kk 恰好表示xt 与xt-k在 排除了其中间变量xt-1, xt-2, , xt-

12、k +1 影响之后的相关系数， xt - k 1 xt-1 - k 2 xt-2 - - kk-1 xt-k +1 = kk xt-k + ut 所以偏自相关函数由此得名。,偏相关函数定义：,kk,样本的偏自相关函数公式,其中，,偏自相关函数计算流程,例:考察如下AR模型的偏自相关图,理论偏自相关系数,样本偏自相关图,理论偏自相关系数,样本偏自相关图,理论偏自相关系数,样本偏自相关图,理论偏自相关系数,样本偏自相关系数图,六、MA(q)偏自相关函数,因为任何一个可逆的MA(q) 过程都可以转换成一个无限阶的系数按几何递减的AR过程，所以MA(q) 过程的偏自相关函数呈缓慢衰减特征。,MA(1)

13、 过程的偏自相关函数呈指数衰减特征。 若1 0, 偏自相关函数呈交替改变符号式 指数衰减； 若1 0，偏自相关函数呈负数的指数衰减。,MA模型的偏自相关系数拖尾,MA模型的偏自相关系数拖尾,ARMA( p, q) 过程的偏自相关函数也是无限延长的，其表现形式与MA(q)过程的偏自相关函数相类似。根据模型中移动平均部分的阶数q以及参数i的不同，偏自相关函数呈指数衰减和（或）正弦衰减混合形式。,七、RAMA(p, q)偏自相关函数,ARMA模型偏自相关系数拖尾性,样本偏自相关图,模型ARMA(1,1),偏相关图,ARMA模型相关性特征,6.3 模型的识别,一、相关图与偏相关图的分析法,模型的识别主

14、要依赖于对相关图与偏相关图的分析。在对经济时间序列进行分析之前，首先应对样本数据取对数，目的是消除数据中可能存在的异方差，然后分析其相关图。 主要识别模型的阶数,识别的第1步是判断随机过程是否平稳 如果一个随机过程是平稳的，其特征方程的根都应在单位圆之外。若特征根接近单位圆，自相关函数将衰减的很慢。所以在分析相关图时，如果发现其衰减很慢，即可认为该时间序列是非平稳的。 这时应对该时间序列进行差分，同时分析差分序列的相关图以判断差分序列的平稳性，直至得到一个平稳的序列。对于经济时间序列，差分次数，即模型中的参数d通常只取0，1或2。,实际中也要防止过度差分。 缺点：（1）序列的样本容量减小；（2

15、）方差变大。 对于一个序列，差分后若数据的极差变大，说明差分过度。,第2步是在平稳时间序列基础上识别ARMA模型阶数p, q。 下面表给出了不同ARMA模型的自相关函数和偏自相关函数。当然一个过程的自相关函数和偏自相关函数通常是未知的。用样本得到的只是估计的自相关函数和偏自相关函数，即相关图和偏相关图。 建立ARMA模型，时间序列的相关图与偏相关图可为识别模型参数p, q提供信息。相关图和偏相关图（估计的自相关系数和偏自相关系数）通常比真实的自相关系数和偏自相关系数的方差要大，并表现为更高的自相关。实际中相关图，偏相关图的特征不会像自相关函数与偏自相关函数那样“规范”，所以应该善于从相关图，偏

16、相关图中识别出模型的真实参数p, q。 另外，估计的模型形式不是唯一的，所以在模型识别阶段应多选择几种模型形式，以供进一步选择。,二、模型定阶经验方法,95的置信区间 模型定阶的经验方法 如果样本(偏)自相关系数在最初的d阶明显大于两倍标准差范围，而后几乎95的自相关系数都落在2倍标准差的范围以内，而且通常由非零自相关系数衰减为小值波动的过程非常突然。这时，通常视为(偏)自相关系数截尾。截尾阶数为d。,时间序列的随机性，是指时间序列各项之间没有相关关系的特征。使用自相关分析图判断时间序列的随机性，一般给出如下准则：,若时间序列的自相关函数基本上都落入 置信区间，则该时间序列具有随机性； 若较多

17、自相关函数落在置信区间之外， 则认为该时间序列不具有随机性。,判断随机性准则,判断时间序列是否平稳，是一项很重要的工作。运用自相关分析图判定时间序列平稳性的准则是：,若时间序列的自相关函数在k3时都落入置 信区间，且逐渐趋于零，则该时间序列具有平稳性； 若时间序列的自相关函数更多地落在置信区 间外面，则该时间序列就不具有平稳性。,判断平稳性准则,运用自相关分析图判定时间序列季节性的准则是：,若某一时间序列在k2或3以后的自相关函数值，存在着显著不为零的值，则该时间序列具有季节性；用差分方法剔除长期趋势，再算自相关函数，用自相关分析图确定季节性。,判断季节性准则,三、基于自相关函数和偏相关函数的

18、定阶方法,对于ARMA(p,q)模型，可以利用其样本的自相关函数和样本偏自相关函数的截尾性判定模型的阶数。,具体方法如下：,对于每一个q,计算,.,（M 取,为 或者 ），考察其中满足,或者,的个数是否占M个的68.3%或者95.5%。,均近似于零,达到其相应的比例，则可以近似地判定,是 步截尾，平稳时间序列,为,。,如果,， 都明显地异于零，而,并且满足上述不等式之一的 的个数,类似，我们可通过计算序列,满足,，考察,或者,的个数是否占M个的68.3%或者95.5%。即可以近似,地判定,是 步截尾，平稳时间序列,为,如果对于序列,和,截尾，即不存在上述的,来说，均不,和,判定平稳时间序列,，

19、则可以,为ARMA模型。,X=6241 7275 8619 11859 22116 20191 26390 19057 . 24498 25021 69643 7859 7053 7106 9764 13818 . 22998 21081 27043 20703 15363 27251 69500 . 12536 7449 6859 9527 10206 18625 20122 27100 . 17865 17097 29653 79071 10644 6971 7696 9182 8543 . 17293 19578 29357 26060 21897 25179 69196 12145 8

20、152 . 7526 10552 12719 17721 21049 28859 22147 23639 24313 . 64842 8487; t=1:60; plot(t,X) for n=1:48 Y(n)=X(n+12)-X(n); end Y figure(2) plot(1:48,Y) figure(3) autocorr(X,59) a=autocorr(X,59) figure(4) plot(1:60,a) figure(5) parcorr(X,59) b=parcorr(X,59) figure(6) plot(1:60,b) figure(7) autocorr(Y,4

21、7) c=autocorr(Y,47) figure(8) plot(1:48,c) figure(9) parcorr(Y,47) d=parcorr(Y,47) figure(10) plot(1:48,d) Z=diff(Y) figure(11) autocorr(Z,46) e=autocorr(Z,46) figure(12) plot(1:47,e) figure(13) parcorr(Z,46) f=parcorr(Z,46) figure(14) plot(1:47,f),6.4 ARMA序列参数估计,待估参数 常用估计方法 矩估计 极大似然估计 最小二乘估计,一、矩估计方

22、法,原理 用样本自相关系数估计总体自相关系数,1、AR(p)模型参数矩估计,Yule-Walker方程,用 代替 ，并求出,则称 为 Yule-Walker估计,AR(1)模型:,Yule-Walker方程,矩估计（Yule-Walker方程的解）,AR(2)模型:,Yule-Walker方程,矩估计（Yule-Walker方程的解）,2、 MA(q)模型参数估计,用 代替 ，并求出 即得矩估计,求解方法：直接法；线性迭代法,MA(1)模型,矩估计,直接法,线性迭代法,改写成,给出一组初值,计算第一次迭代值,继续,直到对于m步与m-1步两组值相差不大，停止迭代。,则解为,注意：模型必须可逆，否

23、则改变初始值，重新迭代,可编程利用MATLAB计算,由于模型结构的复杂性，比较困难，可选用其它方法。,3、ARMA(p,q)模型的参数估计,ARMA(1,1)模型,方程,矩估计,二、最小二乘估计,原理 使残差平方和达到最小的那组参数值即为最小二乘估计值,如设样本序列,定义参差平方和,使 达到极小 的即为最小二乘估计,三、极大似然估计,ARMA(p,q)模型参数的精估计，一般采用极大似然估计，由于模型结构的复杂性，无法直接给出参数的极大似然估计，只能通过迭代方法来完成，这时，迭代初值常常利用矩估计得到的值。,原理 在极大似然准则下，认为样本来自使该样本出现概率最大的总体。因此未知参数的极大似然估

24、计就是使得似然函数（即联合密度函数）达到最大的参数值,例题分析,设,为一AR(2)序列，,其中,。,求,的自协方差函数,。, 例 1,解答：,Yule-Walker方程为：,即：,且：,联合上面三个方程，解出：, 例 2,考虑如下AR(2) 序列：,若已知观测值,（1）试预报,（2）给出（1）预报的置信度为95%的预报区间,解答：,（1）,（2）,预报的置信度为95%的预报区间分别为：,MATLAB命令,函数名：corrcoef 语法：M=corrcoef(x) M=corrcoef(x1,x2),说明： M=corrcoef(x) 计算时间序列x的自相关系数， 结果为1； M=corrcoe

25、f(x1,x2)，计算两时间序列x1,x2的 相关系数M，其值在0，1之间,1、计算两时间序列的相关系数,2、计算并描绘时间序列的自相关函数,函数名：autocorr 语法：autocorr(series,nLags,M,nSTDs) ACF,Lags,Bounds=autocorr(),说明：series：时间序列 nLags：延迟，当nLags= 或缺省时，计算ACF时在延迟 点 0，1，T，而T=min(20,length(series)-1; M：表示Lags的非负整数，当M = 或缺省时，且autocorr 假设为高斯白噪声。 nSTDs：表示计算出的相关函数ACF估计误差的标准差；

26、 ACF：相关函数； Lags：对应于ACF的延迟； Bounds：置信区间的近似上下限，假设序列是MA（M)模型,运行结果： Bounds= 0.0899 -0.0899 autocorr(y,2) %如图所示,例：randn(state,0) x=randn(1000,1);%生成1000点的Gaussian白噪声 y=filter(1 1 1,1,x);%生成MA（2）过程 ACF,Lags,Bounds= autocorr(y,2);%计算95置信度下的 相关系数 Bounds,3、计算并描绘时间序列的偏相关函数,函数名：parcorr 语法：a,b=parcorr(series) P

27、artialACF,Lags,Bounds=Parcorr(series,nLags,R, nSTDs ),说明：series：时间序列;a 阶数，b与阶数对应的偏相关函数; nLags：延迟，当nLags= 或缺省时，计算ACF时在延迟 点 0，1，T，而T=min(20,length(series)-1; R：表示Lags的非负整数，当R = 或缺省时，且parcorr 假设为AR(R)过程。 nSTDs：表示计算出的相关函数ACF估计误差的标准差； ACF：相关函数； Lags：对应于ACF的延迟； Bounds：置信区间的近似上下限，假设序列是AR(R)过程,运行结果： Bounds=

28、 0.0633 -0.0633 parcorr(y,2) %如图所示,例：randn(state,0) x=randn(1000,1); %生成1000点的Gaussian白噪声 y=filter(1,1 0.6 0.08,x); %生成AR(2）过程 PartialACF,Lags,Bounds= Parcorr(y,2); %计算95置信度下的相关系数 Bounds,时间序列模型估计,评价时间序列模型的FPE准则： 最终预报误差（Final Prediction Error）的定阶准则称为FPE准则,AR模型定阶的步骤：,N是样本的个数，P作为AR模型最大的阶,第一步：任意选取正整数P：,

29、第二步：依次计算FPE(h)：h=1,2,P 使得FPE最大的正整数P即为AR模型的阶。,也可用：偏相关系数parcorr(y) ，自相关系数autocorr(y)确定AR模型、MA模型的滞后阶数,AR模型参数估计,函数名：ar （1）m=ar(y,n) (2) m,refl=ar(y,n,approach, window),输入参数：y 是数据结构由iddata函数得到；y= iddata(y)， 后面y是给定的时间序列， n 是AR阶次; approach ：估计时采用的方法，有下面几种： Approachfb:前向后； ls:最小二乘法； yw:Yule-Walker方法; Burg：基

30、于Burg谱估计方法 Window：处理Y中缺失值的方法 Window now：表示观察值中没有缺失值； Window yw：表示Yule-Walker方法处理缺失值 输入参数：m AR模型的文字形式 refl %AR 模型的系数,例：已知正弦含噪声时间序列sinusoid，比较Burg谱估计和前向后无窗方法预测AR模型参数,y=sin(1:300)+0.5*randn(300,1);%原始数据 y=iddata(y); y的数据结构 mb=ar(y,4,burg); Burg谱估计 mfb =ar(y,4); 无窗估计 figure(1); bode(mb,mfb) bode频域相应图 fi

31、gure(2); plot(1:length(y),y)时序图,Discrete-time IDPOLY model: A(q)y(t) = e(t) A(q) = 1 - 0.2864 q-1 + 0.2024 q-2 + 0.4591 q-3 + 0.1087 q-4 Estimated using AR (burg/now) from data set y Loss function 0.35826 and FPE 0.367943 Sampling interval: 1,Discrete-time IDPOLY model: A(q)y(t) = e(t) A(q) = 1 - 0.

32、2841 q-1 + 0.2045 q-2 + 0.458 q-3 + 0.1087 q-4 Estimated using AR (fb/now) from data set y Loss function 0.358254 and FPE 0.367937 Sampling interval: 1,时序图,lpc命令: 功能：计算AR模型的系数 格式：a,g=lpc(x,p) %使用全极点模型自互相关法计算,例：已知某一地区19801999年的肿瘤引起的死亡率（）如下： 10.010 11.260 9.000 9.090 9.440 9.090 8.730 8.680 9.040 9.04

33、5 10.050 7.330 6.190 5.680 5.860 5.630 5.560 5.640 5.700 6.360 利用此数据建立AR模型。,x=10.010 11.260 9.000 9.090 9.440 9.090 8.730. 8.680 9.040 9.045 10.050 7.330 6.190 5.680. 5.860 5.630 5.560 5.640 5.700 6.360; a=lpc(x,3) %模型参数 estx=filter(0-a(2:end),1,x); 估计时间序列 e=x-estx; 预测误差 acs,lags=xcorr(e,coeff); ACS

34、的预测误差 subplot(121); plot(1:20,x,1:20,estx,-.); title(原始信号); xlabel(采样点);ylabel(幅度);grid; legend(原始信号,LPC估计) subplot(122); plot(lags,acs); title(预测误差的自相关函数); xlabel(延迟);ylabel(归一化值);grid;,a = 1.0000 -1.0153 0.1922 -0.1179,ARMAX模型参数估计 考虑外部解释变量X ARMAX model ： A(q) y(t) = B(q) u(t-nk) + C(q) e(t),na,nb,

35、nc是滞后多项式的阶数， nk为延迟,函数名：armax m=armax(data,orders) 输入参数：data：输入输出时间序列 Ordersna,nb,nc,nk Z = IDDATA (y),y原始序列 M = ARMAX(Z,na nb nc nk) M = ARMAX(Z,na,na,nb,nb,nc,nc,nk,nk) na,nb,nc是滞后多项式的阶数， nk为延迟,例：研究深发展收益率（000001）与上证指数收益率之间的关系。如选用的时间段为2005年10月21日到2006年9月29日。深发展的收益率保存在变量y中，上证指数的收益率保存在变量u中， 收益率为算术收益率

36、（t时刻的收益率 ） 用ARMAX模型进行估计。,z=iddata(y,u); m=armax(z,2,2,2,1);,A(q) = 1 1.273 q-1 +0.9398 q-2 B(q) = -0.305 q-1 + 0.2697 q-2 C(q) = 1 1.184 q-1 +0.8833 q-2 损失函数值0.000590723,FPE=0.000626166,用ARMAX模型可对MA模型进行估计,z=iddata(y); m=armax(z,nc,5);,y(t) = C(q) e(t) C(q) = 1+ 0.01409 q-1 +0.05519 q-2- 0.004185 q-3

37、+ 0.09333 q-4+ 0.1789 q-5 损失函数值0.000145294,FPE=0.000152486,A(q) y(t) = B(q) u(t-nk) + C(q) e(t)令A(q) 1， B(q) 0，可得MA模型,用ARMAX模型可对ARMA模型进行估计,na=2;nc=3 orders=na nc z=iddata(y); m=armax(z,orders);,A(q) y(t) = C(q) e(t)A(q) = 1 + 0.1337 q-1 +0.9058 q-2 C(q) = 1+ 0.1239 q-1 +0.9026 q-2- 0.06353 q-3 损失函数值

38、0.000611174,FPE=0.00064143,A(q) y(t) = B(q) u(t-nk) + C(q) e(t)令B(q) 0，可得ARMA模型,ARX模型参数估计 考虑外部解释变量X A(q) y(t) = B(q) u(t-nk) + e(t) M = ARX(DATA,ORDERS) ORDERS = na nb nk m = ARX(DATA,na,na,nb,nb,nk,nk),例：研究深发展收益率（000001）与上证指数收益率之间的关系。如选用的时间段为2005年10月21日到2006年9月29日。深发展的收益率保存在变量y中，上证指数的收益率保存在变量u中， 收益

39、率为算术收益率 （t时刻的收益率 ） 用ARX模型进行估计。,z=iddata(y,u); m=arx(z,na,2,nb,3);,A(q) y(t) = B(q) u(t-nk) + e(t) A(q) = 1 - 0.1118 q-1 - 0.09373 q-2 B(q) = -0.2222 q-1 - 0.2341 q-2 - 0.07306 q-3 损失函数0.000557467,FPE=0.000584268,例：已知系统实测输入与输出数据如表所示，试根据这些数据辨识出系统模型,第一步，利用AIC辨识模型的阶次； 第二步，利用ARX确定具体的模型,u=1.4601 0.8849 1.

40、1854 1.0887 1.413 1.3096 1.0651 0.7184. 1.3571 1.0557 1.1923 1.3335 1.4374 1.2905 0.841 1.0245. 1.4483 1.4335 1.0282 1.4149 0.7463 0.9822 1.3503 0.7078. 0.8111 0.8622 0.8589 1.183 0.9177 0.859 0.7122 1.2974. 1.056 1.4454 1.0727 1.0349 1.3769 1.1201 0.8621 1.2377. 1.3704 0.7157 1.245 1.0035 1.3654 1.

41、1022 1.2675 1.0431;,U=iddata(y,u,0.1); for n=1:7, for m=1:7 T1=arx(U,n,m,0); TAic0(n,m)=aic(T1) T2=arx(U,n,m,1);,y=0 0 8.7606 13.1939 17.41 17.6361 18.7627 18.5296 17.0414. 13.4154 14.4539 14.59 16.1104 17.6853 19.4981 19.5935. 16.4106 14.3359 15.7463 18.1179 17.784 18.8104 15.3086. 13.7004 14.8178

42、13.2354 12.2993 11.6001 11.6074 13.7662. 14.195 13.763 11.8713 13.8566 14.6944 17.8659 17.6543. 16.6386 17.1071 16.5373 14.643 15.0862 16.8058 14.7641. 15.4976 14.679 16.6552 16.6301;,TAic1(n,m)=aic(T2) T3=arx(U,n,m,2); TAic2(n,m)=aic(T3) end, end,不同阶次组合下的AIC准则值 （延迟步数分别为0，1，2）,TAic0 = 1.3484 1.3737

43、-0.2335 -0.6313 -1.0041 -1.5317 -2.6039 1.2379 1.1946 -2.0987 -2.3506 -4.8816 -5.2171 -7.2810 1.0426 1.0427 -2.8730 -3.4498 -5.4326 -5.5666 -7.5037 1.0221 1.0345 -7.8012 -9.9035 -11.4173 -11.4332 -11.6348 1.0079 1.0287 -9.5807 -11.2235 -11.3944 -11.4389 -11.7450 1.0292 1.0575 -11.1253 -11.3542 -11.5

44、753 -11.5470 -11.8350 0.9850 1.0261 -11.3697 -11.3424 -11.5572 -11.6429 -11.8065 合适的阶次是选（4，5，0）,TAic1 = 1.4841 -0.2543 -0.6614 -1.0458 -1.5673 -2.6351 -3.3914 1.3459 -2.1256 -2.3677 -4.9090 -5.2106 -7.3222 -7.4634 1.0658 -2.8874 -3.4736 -5.4459 -5.5889 -7.5454 -7.6691 1.0330 -7.8342 -9.9322 -11.4589

45、 -11.4709 -11.6686 -11.7934 1.0043 -9.5929 -11.2505 -11.4341 -11.4767 -11.7747 -11.7948 1.0230 -11.1669 -11.3952 -11.6094 -11.5794 -11.8711 -11.8424 0.9909 -11.4108 -11.3833 -11.5890 -11.6760 -11.8429 -11.8016 合适的阶次是选（4，4，1）,TAic2 = -0.2911 -0.7030 -1.0810 -1.6028 -2.6764 -3.3990 -3.5630 -2.1664 -2.

46、4094 -4.9496 -5.2515 -7.3297 -7.5019 -9.6802 -2.9280 -3.5094 -5.4840 -5.6178 -7.5797 -7.7073 -10.0269 -7.8648 -9.9713 -11.4714 -11.4780 -11.6853 -11.7937 -11.7925 -9.6345 -11.2570 -11.4468 -11.4824 -11.7949 -11.8010 -11.8077 -11.2008 -11.4010 -11.6174 -11.5929 -11.8887 -11.8538 -11.9053 -11.4221 -11

47、.3945 -11.6050 -11.7030 -11.8537 -11.8124 -11.8673 合适的阶次是选（4，3，2）,第二步：只需取一个合适的阶次即可： U=iddata(y,u,0.1); T1=arx(U,4,5,0)，T2=arx(U,4,4,1)，T3=arx(U,4,3,2),Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t) A(q) = 1 - 1.004 q-1 + 0.2498 q-2 + 0.2528 q-3 - 0.1269 q-4 B(q) = 8.694e-005 - 0.001947 q-1 +

48、 6.004 q-2 - 0.6219 q-3 - 0.1433 q-4 Estimated using ARX from data set U Loss function 7.56226e-006 and FPE 1.10525e-005 Sampling interval: 0.1,Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t) A(q) = 1 - 1.004 q-1 + 0.2498 q-2 + 0.2528 q-3 - 0.1269 q-4 B(q) = -0.001934 q-1 + 6.004 q-2 - 0.622

49、q-3 - 0.1434 q-4 Estimated using ARX from data set U Loss function 7.56269e-006 and FPE 1.05878e-005 Sampling interval: 0.1,Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t) A(q) = 1 - 1.003 q-1 + 0.2496 q-2 + 0.2528 q-3 - 0.1268 q-4 B(q) = 6.004 q-2 - 0.6199 q-3 - 0.143 q-4 Estimated using ARX from data set U Loss function 7.78671e-006 and FPE 1.04456e-005 Sampling interval: 0.1,6.5 模型