时间序列模型的相关概念以及Matlab拟合ARIMA(p,d,q)模型

时间序列模型

统计量：

通常，我们用 X t = { X ∣ t ∈ T } X_t=\{X|t∈T\} Xt​={X∣t∈T}表示一个时间序列,$x_1,x_2,\dots ,x_n$是观察值，时间序列有下面几个常用的统计量，这些统计量，为我们拟合时间序列模型做出准备.这几个统计量分别是均值，方差，协方差，相关系数，自相关系数，偏自相关系数

均值：

理论的均值函数
$${\mu }_t=E(X_t)$$
均值的估计值
$$\hat{\mu }=\frac{1}{n}\sum _{i=1}^n{x}_i$$

方差：

理论的方差函数
$${\sigma }_t^2=D(X_t)=E(X_t^2)-(E(X_t){)}^2$$
方差的估计值
$${\hat{\sigma }}^2=\frac{1}{n-1}\sum _{i=1}^n(x_i-\hat{\mu }{)}^2$$

差分

对序列1阶差分
$$\Delta x_t=x_t-x_{t-1}$$
2阶差分
$${\Delta }^2{x}_t=\Delta x_t-\Delta x_{t-1}=x_t-2x_{t-1}+x_{t-2}$$
P阶差分
$${\Delta }^p{x}_t={\Delta }^{p-1}{x}_t-{\Delta }^{p-1}{x}_{t-1}$$
P步差分
$${\Delta }_p{x}_t=x_t-x_{t-p}$$
延迟算子
$$\begin{gathered} L{x}_t=x_{t-1} \\ {L}^i{x}_t=x_{t-i} \end{gathered}$$
常见的时间序列模型：
$$x_t={\varphi }_1{x}_{(}t-1)+{\varphi }_2{x}_{(}t-2)+\cdots +{\varphi }_p{x}_{(}t-p)+{\epsilon }_t$$
用延迟算子表示是
$$(1-{\varphi }_{1}L-{\varphi }_2{L}^2-\dots -{\varphi }_p{L}^p)x_t={\epsilon }_t$$

平稳时间序列

白噪声序列

一个时间序列，如果完全符合正态分布，说明这个序列是白噪声序列，一旦时间序列通过了白噪声检验，说明这个序列已经没有了分析的价值，我们应该停止分析

白噪声序列满足条件
$$E(X_t)=\mu$$

$$\gamma (s,t)=\{\begin{array}{l}{\sigma }^2,s\ne t\\ 0,s=t\end{array}$$
记为 $X_t\sim WN(\mu ,{\sigma }^2)$

白噪声检验

我们不加证明的给出预备理论：

1.如果一个序列是白噪声序列，那么这个序列的以非零延迟的自相关系数${\hat{\rho }}_k$近似服从方差为观察期数的倒数，均值为0的正态分布
$${\hat{\rho }}_k\sim N(0,\frac{1}{n}),\forall k\ne 0$$
2.若 $Y_i\sim N(0,1),i=1,2,\dots ,n$ ,则
$$\sum _{i=1}^n{Y}_i^2\sim {\chi }^2(n)$$
N个标准正态分布的平方和服从自由度n的卡方分布

原假设：延迟期数小于m的序列完全不相关

备择假设：延迟期数小于m的序列存在相关性

我们构造原假设 $H_0:{\rho }_0={\rho }_1=\cdots ={\rho }_m$

备择假设 $H_1:{\rho }_0,{\rho }_1,\dots ,{\rho }_m$ 至少存在非零值

白噪声检验常用的统计量有两个

1. Q统计量

$$Q=n\sum _{i=1}^m{\hat{\rho }}_k^2$$
n是序列观测期数，m是指定的延期数

由于$${\hat{\rho }}_k\sim N(0,\frac{1}{n}),\forall k\ne 0$$，所以 $$\sqrt{n}{\hat{\rho }}_k\sim N(0,1)$$ ,所以
$$Q=n\sum _{i=1}^m{\hat{\rho }}_k^2\sim {\chi }^2(m)$$
若 $Q\ge {\chi }_{1-\alpha }^2(m)$的分位点，或者Q的p值小于 $\alpha$ ,可以在$1-\alpha$的置信水平下拒绝原假设，认为原序列不是白噪声序列，存在相关性，可以继续分析。
Q统计量适合大样本的情形

2. LB统计量

   $$LB=n(n+2)\sum _{k=1}^m(\frac{{\hat{\rho }}_k^2}{n-k})$$
   n是序列观测期数，m是指定的延期数
   我们不加证明的告诉大家LB近似服从自由度为m的卡方分布
   $$LB\sim {\chi }^2(m)$$

自协方差相关系数(ACF):

$${\rho }_{s,t}=Corr(X_s,X_t)=\frac{Cov(X_s,X_t)}{D(X_s)D(X_t)}=\frac{{\gamma }_{s,t}}{{\sigma }_s{\sigma }_t}$$

回顾一下协方差和相关系数的定义

$$\begin{gathered} Cov(X,Y)=E(X-E(X))(Y-E(Y)) \\ Corr(X,Y)=Cov(X,Y)/D(X)D(Y) \\ {\rho }_{s,s}=1 \end{gathered}$$

上面的定义，意味着在 时刻要大量的数据，才可以计算ACF，很明显这是理想情况，这是不现实的。为此，我们引入平稳时间序列的概念

平稳序列满足三个条件

(1) $\forall t\in T,E(X_t)<+\infty$

(2) $\forall t\in T,E(X_t)=\mu ,\mu$, 为常数

(3) $\forall r,s,t,{\gamma }_{s,t}={\gamma }_{s+r,t+r}$

可以看出，平稳时间序列中，任意两个时间点的ACF,只和这两个点之间的距离,即滞后 $|t-s|$有关，即${\gamma }_{s,t}={\gamma }_{0,t-s}$
这样的话，如果想要计算ACF，我们需要计算的就是
偏自相关系数比较复杂，在此从略计算，大家会代码即可
$${\rho }_{lag}=\frac{{\gamma }_{lag}}{{\gamma }_0},lag=0,1,2,\dots$$
，

只需要知道，AR§的PACF是p阶截尾，MA(q)的ACF是q阶截尾即可

定性理论

差分方程通解问题

$$(1+a_{1}L+a_2{L}^2+\cdots +a_p{L}^p)z_t=h(t)$$
若$h(t)=0$ ,则上面的方程是齐次差分方程

齐次差分方程通解问题

上面的差分方程的特征方程为${\lambda }^p+a_1{\lambda }^{p-1}+a_2{\lambda }^{p-2}+\cdots +a_p=0$

根据代数学基本定理，这个p次多项式必定有p个复根 ${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_p$ ，下面我们根据根的情况对解进行讨论

1.${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_p$ 均为互不相等实根

齐次差分方程的解为 $z_t=c_1{\lambda }_1^t+\cdots +c_p{\lambda }_p^t$

${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_p$ 均为实数，但存在相等的实数根

2.${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_d$ 是相等的实根，其他的实根互不相同

齐次差分方程的解为 $$z_t=(c_1+c_{2}t+\cdots +c_d{t}^{d-1}){\lambda }_1^t+c_{d+1}{\lambda }_{(}d+1{)}^t+\cdots +c_p{\lambda }_p^t$$

3.${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_p$存在复数

由于系数均为实数，所以复根必定共轭成对出现，设
$${\lambda }_1=a+ib=r{e}^{i\omega },{\lambda }_2=a-ib=r{e}^{-i\omega },r=\sqrt{a^2+b^2},\omega =arccos\frac{a}{r}$$
${\lambda }_3,{\lambda }_2,\dots ,{\lambda }_p$均为互不相等实根，此时

齐次差分方程的解为

$$z_t=r^t((c_1{e}^{it\omega }+c_2{e}^{-it\omega })+c_3{\lambda }_3^t+\cdots +c_p{\lambda }_p^t$$
非齐次差分方程通解问题

非齐次通解 $z_t$是齐次通解 $z_t^{\prime }$与非齐次特解 $z_t^{\prime \prime }$之和$z_t=z_t^{\prime }+z_t^{\prime \prime }$

AR( p )

AR是自回归模型的简称，AR模型的定义和性质在下面
$$\{\begin{array}{l}{x}_t={\varphi }_0+{\varphi }_1{x}_{t-1}+\dots +{\varphi }_p{x}_{t-p}+{\epsilon }_t\\ {\varphi }_p\ne 0\\ E\left({\epsilon }_t\right)=0,\mathrm{Var}\left({\epsilon }_t\right)={\sigma }_t^2,E\left({\epsilon }_t{\epsilon }_s\right)=0,s\ne t\\ E\left(X_t{\epsilon }_t\right)=0,\forall s<t\end{array}$$
用延迟算子表示为 $\Phi (L)x_t={\epsilon }_t$

中心化：
$$\mu =\frac{{\varphi }_0}{1-{\varphi }_1-\dots -{\varphi }_p},y_t=x_t=\mu$$
{ y t } \{y_t \} {yt​} 是 { x t } \{x_t \} {xt​} 的中心化序列

AR( p)的平稳性判别

中心化的序列平稳的条件是
$$t\to 0,x_t\to 0$$
观察齐次差分方程的解$$z_t=r^t(c_1{e}^{it\omega }+c_2{e}^{-it\omega })+c_3{\lambda }_3^t+\cdots +c_p{\lambda }_p^t$$，可以得出，序列平稳的等价条件是 $|{\lambda }_i|<1,i=1,2,\dots ,p$

这就是检验平稳性的特征根判别法

重要性质：

自相关系数（ACF）拖尾，偏自相关系数（PACF）p阶截尾

AR§的ACF拖尾的原因可以直观解释，$x_t={\varphi }_1{x}_{t-1}+{\varphi }_2{x}_{t-2}+\cdots +{\varphi }_p{x}_p+{\epsilon }_t$ 这个模型里面,$x_t$的值受到后面p期序列值的影响，但是$x_{t-1}$的序列值又受$x_{t-1-p}$的影响，这样下去,${\rho }_{0,1},{\rho }_{0,2},\dots {\rho }_{0,t}$ 的值以指数形式减少，所以是拖尾,

AR( p)的偏自相关系数p阶截尾，这个结论来自Yule-Walker方程

如果我们想计算滞后为k的偏自相关系数，我们需要知道$${\rho }_1,{\rho }_2,\dots ,{\rho }_{k-1},{\rho }_k$$ 的值，然后根据下面的方程组(Yule-Walker方程)求解
$${\rho }_l={\varphi }_{k1}{\rho }_{l-1}+{\varphi }_{k2}{\rho }_{l-2}+\dots +{\varphi }_{kk}{\rho }_{l-k},\forall l,1⩽l⩽k$$
用矩阵表示
$$\left(\begin{array}{cccc}1& {\rho }_1& \cdots & {\rho }_{k-1}\\ {\rho }_1& 1& \cdots & {\rho }_{k-2}\\ ⋮& ⋮& & ⋮\\ {\rho }_{k-1}& {\rho }_{k-2}& \cdots & 1\end{array}\right)\left(\begin{array}{c}{\varphi }_{k1}\\ {\varphi }_{k2}\\ ⋮\\ {\varphi }_{kk}\end{array}\right)=\left(\begin{array}{c}{\rho }_1\\ {\rho }_2\\ ⋮\\ {\rho }_k\end{array}\right)$$
这是一个线性方程组，根据Cramer法则就可以求解， 就是参数的解， 就是滞后为k的偏自相关系数
$${\varphi }_{kk}=\frac{D_k}{D}$$

$$D=\left|\begin{array}{cccc}1& {\rho }_1& \cdots & {\rho }_{k-1}\\ {\rho }_1& 1& \cdots & {\rho }_{k-2}\\ ⋮& ⋮& & ⋮\\ {\rho }_{k-1}& {\rho }_{k-2}& \cdots & 1\end{array}\right|,D_k=\left|\begin{array}{cccc}1& {\rho }_1& \cdots & {\rho }_1\\ {\rho }_1& 1& \cdots & {\rho }_2\\ ⋮& ⋮& & ⋮\\ {\rho }_{k-1}& {\rho }_{k-2}& \cdots & {\rho }_k\end{array}\right|$$

PACF截尾的原因是，当k>p时， $D_k=0$

下面简单证明一下这个结论

由于AR( p)模型的${\varphi }_p\ne 0$ ,向量 $\eta$可以表示成${\xi }_i(i=1,2,\dots ,p)$ 的非零线性组合,k>p时，$D_k$的前p个列向量正好是${\xi }_i(i=1,2,\dots ,p)$,最后一个列向量正好是$\eta$,所以$D_k=0,{\varphi }_{k}k=0$.

所以，AR( p)的PACF是p阶截尾

MA(q)

定义：具有下面的结构的模型是q阶移动平均模型，记为MA(q)
$$\{\begin{array}{l}{x}_t=\mu +{\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-{\theta }_2{\epsilon }_{t-2}-\dots -{\theta }_q{\epsilon }_{t-q}\\ {\theta }_p\ne 0\\ E\left({\epsilon }_t\right)=0,D\left({\epsilon }_t\right)={\sigma }_{\epsilon }^2,E\left({\epsilon }_t{\epsilon }_s\right)=0,s\ne t\end{array}$$
​

用延迟算子表示为 $x_t=\Theta (L){\epsilon }_t$

中心化： $y_t=x_t-\mu$

MA(q)任何时候都平稳

可逆性条件：

系数多项式 $\Theta (L)=0$的根都在单位圆外，这跟AR( p)的平稳性条件( $\Phi (B)=0$的根都在单位圆外)完全对偶

重要性质：

1.自协方差函数只和滞后阶数有关，是q阶截尾,自相关系数q阶截尾，${\rho }_k=\frac{{\gamma }_k}{{\gamma }_0}$
$${\gamma }_k=E(x_t{x}_{t-k})=E[({\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-\dots -{\theta }_q{\epsilon }_{t-q})({\epsilon }_{t-k}-{\theta }_1{\epsilon }_{t-k-1}-\dots -{\theta }_q{\epsilon }_{t-k-q})]$$

$$=\{\begin{array}{l}\left(1+{\theta }_1^2+\dots +{\theta }_q^2\right){\sigma }_{\epsilon }^2,k=0\\ \left(-{\theta }_k+{\sum }_{i=1}^{q-k}{\theta }_i{\theta }_{k+i}\right){\sigma }_{\epsilon }^2,1⩽k⩽q\\ 0,k>q\end{array}$$

2.偏自相关系数拖尾

这个结论来自MA(q)模型的可逆性，$MA(q)$模型可以写成$AR(\infty )$的形式,有AR§模型的PACF是p阶结尾，可以得出 $AR(\infty )$ 模型是无穷阶截尾，也就是拖尾，所以得出MA(q)模型的PACF拖尾的结论

ARMA(p,q)

定义：具有下面结构的模型叫自回归移动平均模型，记为ARMA(p,q)
$$\{\begin{array}{l}{x}_t={\varphi }_0+{\varphi }_1{x}_{t-1}+\dots +{\varphi }_p{x}_{t-p}+{\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-{\theta }_2{\epsilon }_{t-2}-\dots -{\theta }_q{\epsilon }_{t-q}\\ {\varphi }_p\ne 0,{\theta }_q\ne 0\\ E\left({\epsilon }_t\right)=0,D\left({\epsilon }_t\right)={\sigma }_{\epsilon }^2,E\left({\epsilon }_t{\epsilon }_s\right)=0,s\ne t\\ E\left(x_s{\epsilon }_t\right)=0,\forall s<t\end{array}$$
用延迟算子表示为$\Phi (L)x_t=\Theta (L){\epsilon }_t$

p=0,退化为MA(q),q=0,退化为AR( p)

平稳性条件：ARMA(p,q)的平稳性条件只由自回归部分决定，即$\Phi (B)=0$的根全在单位圆外

统计性质：自相关系数，偏自相关系数都不截尾

非平稳时间序列

差分运算

差分运算是未来使序列变得平稳

ARIMA(p,d,q)

对于一个非平稳序列，差分运算可以让序列变得平稳，所以引入ARIMA模型，通俗的说，是d阶差分运算后，变为ARMA(p,q)模型
$$\{\begin{array}{l}\Phi (L){\Delta }^d{x}_t=\Theta (L){\epsilon }_t\\ E\left({\epsilon }_t\right)=0,D\left({\epsilon }_t\right)={\sigma }_{\epsilon }^2,E\left({\epsilon }_t{\epsilon }_s\right)=0,s\ne t\\ E\left(x_s{\epsilon }_t\right)=0,\forall s<t\end{array}$$
${\Delta }^d=(1-L{)}^d$,是d阶差分

例题

例1

现有201个生产数据，拟合出合适的模型

81.9	89.4	79	81.4	84.8	85.9	88	80.3	82.6
83.5	80.2	85.2	87.2	83.5	84.3	82.9	84.7	82.9
81.5	83.4	87.7	81.8	79.6	85.8	77.9	89.7	85.4
86.3	80.7	83.8	90.5	84.5	82.4	86.7	83	81.8
89.3	79.3	82.7	88	79.6	87.8	83.6	79.5	83.3
88.4	86.6	84.6	79.7	86	84.2	83	84.8	83.6
81.8	85.9	88.2	83.5	87.2	83.7	87.3	83	90.5
80.7	83.1	86.5	90	77.5	84.7	84.6	87.2	80.5
86.1	82.6	85.4	84.7	82.8	81.9	83.6	86.8	84
84.2	82.8	83	82	84.7	84.4	88.9	82.4	83
85	82.2	81.6	86.2	85.4	82.1	81.4	85	85.8
84.2	83.5	86.5	85	80.4	85.7	86.7	86.7	82.3
86.4	82.5	82	79.5	86.7	80.5	91.7	81.6	83.9
85.6	84.8	78.4	89.9	85	86.2	83	85.4	84.4
84.5	86.2	85.6	83.2	85.7	83.5	80.1	82.2	88.6
82	85	85.2	85.3	84.3	82.3	89.7	84.8	83.1
80.6	87.4	86.8	83.5	86.2	84.1	82.3	84.8	86.6
83.5	78.1	88.8	81.9	83.3	80	87.2	83.3	86.6
79.5	84.1	82.2	90.8	86.5	79.7	81	87.2	81.6
84.4	84.4	82.2	88.9	80.9	85.1	87.1	84	76.5
82.7	85.1	83.3	90.4	81	80.3	79.8	89	83.7
80.9	87.3	81.1	85.6	86.6	80	86.6	83.3	83.1
82.3	86.7	80.2

z=[81.9000000000000	89.4000000000000	79	81.4000000000000	84.8000000000000	85.9000000000000	88	80.3000000000000	82.6000000000000	83.5000000000000	80.2000000000000	85.2000000000000	87.2000000000000	83.5000000000000	84.3000000000000	82.9000000000000	84.7000000000000	82.9000000000000	81.5000000000000	83.4000000000000	87.7000000000000	81.8000000000000	79.6000000000000	85.8000000000000	77.9000000000000	89.7000000000000	85.4000000000000	86.3000000000000	80.7000000000000	83.8000000000000	90.5000000000000	84.5000000000000	82.4000000000000	86.7000000000000	83	81.8000000000000	89.3000000000000	79.3000000000000	82.7000000000000	88	79.6000000000000	87.8000000000000	83.6000000000000	79.5000000000000	83.3000000000000	88.4000000000000	86.6000000000000	84.6000000000000	79.7000000000000	86	84.2000000000000	83	84.8000000000000	83.6000000000000	81.8000000000000	85.9000000000000	88.2000000000000	83.5000000000000	87.2000000000000	83.7000000000000	87.3000000000000	83	90.5000000000000	80.7000000000000	83.1000000000000	86.5000000000000	90	77.5000000000000	84.7000000000000	84.6000000000000	87.2000000000000	80.5000000000000	86.1000000000000	82.6000000000000	85.4000000000000	84.7000000000000	82.8000000000000	81.9000000000000	83.6000000000000	86.8000000000000	84	84.2000000000000	82.8000000000000	83	82	84.7000000000000	84.4000000000000	88.9000000000000	82.4000000000000	83	85	82.2000000000000	81.6000000000000	86.2000000000000	85.4000000000000	82.1000000000000	81.4000000000000	85	85.8000000000000	84.2000000000000	83.5000000000000	86.5000000000000	85	80.4000000000000	85.7000000000000	86.7000000000000	86.7000000000000	82.3000000000000	86.4000000000000	82.5000000000000	82	79.5000000000000	86.7000000000000	80.5000000000000	91.7000000000000	81.6000000000000	83.9000000000000	85.6000000000000	84.8000000000000	78.4000000000000	89.9000000000000	85	86.2000000000000	83	85.4000000000000	84.4000000000000	84.5000000000000	86.2000000000000	85.6000000000000	83.2000000000000	85.7000000000000	83.5000000000000	80.1000000000000	82.2000000000000	88.6000000000000	82	85	85.2000000000000	85.3000000000000	84.3000000000000	82.3000000000000	89.7000000000000	84.8000000000000	83.1000000000000	80.6000000000000	87.4000000000000	86.8000000000000	83.5000000000000	86.2000000000000	84.1000000000000	82.3000000000000	84.8000000000000	86.6000000000000	83.5000000000000	78.1000000000000	88.8000000000000	81.9000000000000	83.3000000000000	80	87.2000000000000	83.3000000000000	86.6000000000000	79.5000000000000	84.1000000000000	82.2000000000000	90.8000000000000	86.5000000000000	79.7000000000000	81	87.2000000000000	81.6000000000000	84.4000000000000	84.4000000000000	82.2000000000000	88.9000000000000	80.9000000000000	85.1000000000000	87.1000000000000	84	76.5000000000000	82.7000000000000	85.1000000000000	83.3000000000000	90.4000000000000	81	80.3000000000000	79.8000000000000	89	83.7000000000000	80.9000000000000	87.3000000000000	81.1000000000000	85.6000000000000	86.6000000000000	80	86.6000000000000	83.3000000000000	83.1000000000000	82.3000000000000	86.7000000000000	80.2000000000000];
figure
subplot(2,2,[1,3])%创建一个2行2列的窗格，把时序图放在第一个窗格和第三个窗格里面，也就是整个左侧
plot(z)
subplot(2,2,2)%创建一个2行2列的窗格，右上角是时间序列的自相关系数图（ACF）
autocorr(z)
subplot(2,2,4)%创建一个2行2列的窗格，右下角是时间序列的偏自相关系数（PACF）
parcorr(z)
%求解当前序列的ACF,PACF

运行结果：

AR3=arima(3,0,0);
resAR3=estimate(AR3,z')%由于estimate函数只能传入列向量，所以对z转置
%根据ACF,PACF拟合模型

运行结果：

拟合出的AR(3)模型是$$y_t=160.46-0.41438y_{t-1}-0.29869y_{t-2}-0.19421y_{t-3}+{\epsilon }_t$$

step=7;
Fore=forecast(resAR3,step,'Y0',z')
%预测后面7期数据

求解结果：

例2

1867~1938英国绵羊数量如下表所示，分析数据

2203	2360	2254	2165	2024	2078	2214	2292	2207	2119	2119	2137
2132	1955	1785	1747	1818	1909	1958	1892	1919	1853	1868	1991
2111	2119	1991	1859	1856	1924	1892	1916	1968	1928	1898	1850
1841	1824	1823	1843	1880	1968	2029	1996	1933	1805	1713	1726
1752	1795	1717	1648	1512	1338	1383	1344	1384	1484	1597	1686
1707	1640	1611	1632	1775	1850	1809	1653	1648	1665	1627	1791

z1=[2203	2360	2254	2165	2024	2078	2214	2292	2207	2119	2119	2137
2132	1955	1785	1747	1818	1909	1958	1892	1919	1853	1868	1991
2111	2119	1991	1859	1856	1924	1892	1916	1968	1928	1898	1850
1841	1824	1823	1843	1880	1968	2029	1996	1933	1805	1713	1726
1752	1795	1717	1648	1512	1338	1383	1344	1384	1484	1597	1686
1707	1640	1611	1632	1775	1850	1809	1653	1648	1665	1627	1791];
z=reshape(z1',[72,1]);
plot(z)
%观察可以发现，z不平稳，所以进行1阶差分
%观察序列图

h1=adftest(diff(z))
plot(diff(z,1))
%通过观察法和ADF检验，可以发现z平稳
%序列不平稳，1阶差分后利用ADF检验观察是否平稳

figure
subplot(2,2,[1,3])
plot(diff(z))
subplot(2,2,2)
autocorr(diff(z))
subplot(2,2,4)
parcorr(diff(z))
%1阶差分后得到平稳序列，求解ACF和PACF以判断模型

z=reshape(z1',[72,1]);%把题目给的表格，重构成列向量
u = iddata(diff(z));
test = [];
for p = 1:5 %自回归对应PACF,给定滞后长度上限p和q，一般取为T/10、ln(T)或T^(1/2),这里取T/10=12
    for q = 1:5 %移动平均对应ACF
        m = armax(u,[p q]);
        AIC = aic(m); %armax(p,q),计算AIC
        test = [test;p q AIC];
    end
end
for k = 1:size(test,1)
    if test(k,3) == min(test(:,3)) %选择AIC值最小的模型
        p_test = test(k,1);
        q_test = test(k,2);
        break;
    end
end
p_test
q_test
%ACF,PACF不截尾，模型应该为ARIMA(p,1,q),根据AIC准则确定阶数

ARIMA313=arima(3,1,3)
resARIMA313=estimate(ARIMA313,z)
%p=3,q=3,拟合ARIMA(3,1,3)

拟合出的ARIMA模型为$$\Delta y_t=-5.5017+0.35941y_{t-1}+0.47286y_{t-2}-0.45233y_{t-3}+0.074145{\epsilon }_{t-1}-0.89635{\epsilon }_{t-2}-0.17779{\epsilon }_{t-3}+{\epsilon }_t$$

step=7;
Fore=forecast(resARIMA313,step,'Y0',z)
%预测后面7期数据