时间序列

1. 确定性时间序列分析方法

1.1 移动平均法
设观测序列为$y_1,\cdots ,y_T$，取移动平均的项数N<T。一次移动平均值计算公式为$$M_t^{(1)}=\frac{1}{N}(y_t+y_{t-1}+\cdots +y_{t-N+1})=M_{t-1}^{(1)}+\frac{1}{N}(y_t-y_{t-N})$$二次移动平均值计算公式为$$M_t^{(2)}=\frac{1}{N}(M_t^{(1)}+\cdots +M_{t-N+1}^{(1)})=M_{t-1}^{(2)}+\frac{1}{N}(M_t^{(1)}-M_{t-N}^{(1)})$$当预测目标的基本趋势是在某一水平上下波动时，可用一次移动平均方法建立预测模型。其预测标准误差为$$S=\sqrt{\frac{{\sum }_{t=N+1}^T({\hat{y}}_t-y_t{)}^2}{T-N}}$$将最近N期序列值得平均值作为未来各期的预测结果，一般N取值范围$5⩽N⩽200$。选择最佳N值的一个有效方法是，比较若干模型的预测误差，预测标准差最小者为好。
当预测目标的基本趋势与某一线性模型相吻合时，常用二次移动平均法，但序列同时存在线性趋势和周期波动时，可用趋势移动平均法建立预测模型$${\hat{y}}_{T+m}=a_T+b_{T}m,m=1,2,\cdots$$式中$a_T=2M_T^{(1)}-M_T^{(2)},b_T=\frac{2}{N-1}(M_T^{(1)}-M_T^{(2)})$。
移动平均法只适合做近期预测，而且是预测目标的发展趋势变化不大的情况。如果目标的发展趋势存在其他的变化，采用简单的移动平均法就会产生较大的预测偏差和滞后。

1.2 指数平滑法
一次移动平均认为最近N期数据对未来值影响相同，都加权1/N，而N期以前的数据对未来值没有影响，加权为0。二次及更高次移动平均的权数却不是1/N，两端权数小，中间权数大，不符合一般系统的动态性。更切实际的方法是对各期观测值依时间顺序进行加权平均作为预测值，可使用指数平滑法。

1.2.1 一次指数平滑法

• 预测模型：设时间序列为$y_1,y_2,\cdots ,y_t,\cdots ,\alpha$为加权系数，$0<\alpha <1$，一次指数平滑公式为$$S_t^{(1)}=\alpha y_t+(1-\alpha )S_{t-1}^{(1)}$$以这种平滑值进行预测，就是一次指数平滑法，预测模型为${\hat{y}}_{t+1}=S_t^{(1)}$，即$${\hat{y}}_{t+1}=\alpha y_t+(1-\alpha ){\hat{y}}_t$$
• 加权系数的选择：由预测公式看出，$\alpha$值越大，新数据所占的比重就越大。将公式改写为$${\hat{y}}_{t+1}={\hat{y}}_t+\alpha (y_t-{\hat{y}}_t)$$可看出$\alpha$值越大，修正幅度越大。实际中应多取几个$\alpha$值进行试算，哪个预测误差小就采用哪个。
• 初始值的确定：用一次指数平滑法进行预测，除了选择合适的$\alpha$外，还要确定初始值$s_0^{(1)}$，初始值是由预测者估计或指定的。当时间序列的数据较多时（20个以上），初始值对以后的预测值影响很小，可选用第一期的数据为初始值。若时间序列的数据较少，一般以最初几期实际值的平均值作为初始值。

1.2.2 二次指数平滑法
一次指数平滑法虽然克服了移动平均法的缺点，但当时间序列的变动出现直线趋势时，用一次指数平滑法进行预测，仍存在明显的滞后偏差，因此也必须加以修正。二次指数平滑法利用滞后偏差的规律建立直线趋势模型，计算公式为$$S_t^{(2)}=\alpha S_t^{(1)}+(1-\alpha )S_{t-1}^{(2)}$$当时间序列 { y t } \left\{y_t\right\} {yt​}从某时期开始具有直线趋势时，可用直线趋势模型进行预测$${\hat{y}}_{t+m}=a_t+b_{t}m,m=1,2,\cdots$$式中$a_t=2S_t^{(1)}-S_t^{(2)},b_t=\frac{\alpha }{1-\alpha }(S_t^{(1)}-S_t^{(2)})$。

1.2.3 三次指数平滑法
当时间序列的变动表现为二次曲线趋势时，需要用三次指数平滑法，是在二次指数平滑的基础上，在进行一次平滑，计算公式为$$S_t^{(3)}=\alpha S_t^{(2)}+(1-\alpha )S_{t-1}^{(3)}$$三次指数平滑法的预测模型为$${\hat{y}}_{t+m}=a_t+b_{t}m+C_t{m}^2,m=1,2,\cdots$$式中$a_t=3S_t^{(1)}-3S_t^{(2)}+S_t^{(3)},b_t=\frac{\alpha }{2(1-\alpha {)}^2}[(6-5\alpha )S_t^{(1)}-2(5-4\alpha )S_t^{(2)}+(4-3\alpha )S_t^{(3)}],c_t=\frac{{\alpha }^2}{2(1-\alpha {)}^2}[S_t^{(1)}-2S_t^{(2)}+S_t^{(3)}]$。

指数平滑预测模型是以时刻t为起点，综合历史序列的信息，对未来进行预测。选择合适的加权参数$\alpha$是提高预测精度的关键环节，经验取值范围一般为0.1~0.3。若序列的基本趋势较稳，预测偏差由随机因素造成，则$\alpha$应取小一些，使预测模型能包含更多历史数据的信息。若预测目标的基本趋势已发生系统的变化，则$\alpha$值应取大一些，使预测模型适应预测目标的新变化。
由于指数平滑公式是递推公式，所以必须确定初始值$S_0^{(1)},S_0^{(2)},S_0^{(3)}$，可以取前3~5个数据的算术平均值作为初始值。

1.3 差分指数平滑法
当时间序列的变动具有直线趋势时，用一次指数平滑法会出现滞后偏差，其原因在于数据不满足模型要求。因此，也可以从数据变换的角度来考虑改进措施，差分方法是改变数据变动趋势的简易方法。

1.3.1 一阶差分指数平滑法
当时间序列呈直线增加时，可使用一阶差分指数平滑模型来预测，公式如下$$\begin{gathered} \nabla y_t=y_t-y_{t-1} \\ \nabla {\hat{y}}_{t+1}=\alpha \nabla y_t+(1-\alpha )\nabla {\hat{y}}_t \\ {\hat{y}}_{t+1}=\nabla {\hat{y}}_{t+1}+y_t \end{gathered}$$式中$\nabla$为差分记号。

dyt=diff(yt);		%求yt的一阶向前差分

1.3.2 二阶差分指数平滑模型
当时间序列呈二次曲线增加时，可使用二阶差分指数平滑模型来预测，公式如下$$\begin{gathered} \nabla y_t=y_t-y_{t-1} \\ {\nabla }^2{y}_t=\nabla y_t-\nabla y_{t-1} \\ {\nabla }^2{\hat{y}}_{t+1}=\alpha {\nabla }^2{y}_t+(1-\alpha ){\nabla }^2{\hat{y}}_t \\ {\hat{y}}_{t+1}={\nabla }^2{\hat{y}}_{t+1}+\nabla y_t+y_t \end{gathered}$$式中${\nabla }^2$为二阶差分。

差分方法和指数平滑法的联合应用，除了能克服一次指数平滑法的滞后偏差之外，对初始值的问题也有显著的改进。数据经过差分处理后，所产生的新序列基本上是平稳的，此时初始值取新序列的第一期数据对未来预测值不会有多大影响。但是，对于指数平滑法的加权系数$\alpha$的选择问题，以及只能逐期预测问题，差分指数平滑模型并未改进。

1.4 具有季节性特点的时间序列的预测
对于季节性时间序列的预测，要从数学上完全拟合其变化曲线是非常困难的。一种名为季节系数法的具体步骤如下：

1. 收集m年的每年各季度（每年n个季度）的时间序列样本数据$a_{ij}$。其中i表示年份的序号，j表示季度的序号。
2. 计算每年所有季度的算术平均值$\bar{a}$，即$$\bar{a}=\frac{1}{k}\sum _{i=1}^m\sum _{j=1}^n{a}_{ij},k=mn$$
3. 计算同季度数据的算术平均值$${\bar{a}}_{\cdot j}=\frac{1}{m}\sum _{i=1}^m{a}_{ij}$$
4. 计算季度系数$b_j={\bar{a}}_{\cdot j}/\bar{a}$
5. 预测计算：当时间序列是按季度列出时，先求出预测年份（下一年）的年加权平均$$y_{m+1}=\frac{{\sum }_{i=1}^m{w}_{i}y}{{\sum }_{i=1}^m{w}_i}$$式中$y_i={\sum }_{j=1}^n{a}_{ij}$为第i年的年合计数；$w_i$为第i年的权数，按自然数列取值，即$w_i=i$。
   再计算预测年份的季度平均值$${\bar{y}}_{m+1}=y_{m+1}/n$$最后预测年份第j季度的预测值为$$y_{m+1,j}=b_j{\bar{y}}_{m+1}$$

2. 平稳时间序列模型

平稳时间序列的统计特性不随时间的平移而变化，即均值和协方差不随时间的平移而变化。

2.1 时间序列的基本概念

• 定义1：设随机序列 { X t , t = 0 , ± 1 , ± 2 , ⋯   } \left\{X_t,t=0,\pm1,\pm2,\cdots\right\} {Xt​,t=0,±1,±2,⋯}满足$$\begin{gathered} (1)E(X_t)=\mu =常数 \\ (2){\gamma }_{t+k,t}={\gamma }_k(k=0,\pm 1,\pm 2,\cdots )与t无关 \end{gathered}$$则称$X_t$为平稳随机序列（平稳时间序列），简称平稳序列。其中${\gamma }_{t,s}=\mathrm{Cov}(X_t,X_s)$，称为自协方差函数。将${\gamma }_{t,s}$标准化后，为自相关函数$${\rho }_{t,s}=\frac{{\gamma }_{t,s}}{\sqrt{{\gamma }_{t,t}}\sqrt{{\gamma }_{s,s}}}=\frac{{\gamma }_{t,s}}{{\sigma }_t{\sigma }_s}$$

• 定义2：设平稳序列 { ε t , t = 0 , ± 1 , ± 2 , ⋯   } \left\{\varepsilon_t,t=0,\pm1,\pm2,\cdots\right\} {εt​,t=0,±1,±2,⋯}的自协方差函数为$${\gamma }_k={\sigma }^2{\delta }_{k,0}=\{\begin{array}{l}0,k\ne 0\\ {\sigma }^2,k=0\end{array}$$则称该序列为平稳白噪声序列。平稳白噪声序列的方差是常数${\sigma }^2$，且任意两个不同时点之间是不相关的，是一种最基本的平稳序列。
  平稳序列的样本自协方差函数有如下两种形式$${\hat{\gamma }}_k=\frac{1}{n}\sum _{t=1}^{n-k}(X_{t+k}-\bar{X})(X_t-\bar{X})或{\hat{\gamma }}_k^{\ast }=\frac{1}{n-k}\sum _{t=1}^{n-k}(X_{t+k}-\bar{X})(X_t-\bar{X}),0⩽k⩽n-1$$样本自相关函数为$${\hat{\rho }}_k=\frac{{\hat{\gamma }}_k}{{\hat{\gamma }}_0}或{\hat{\rho }}_k^{\ast }=\frac{{\hat{\gamma }}_k^{\ast }}{{\hat{\gamma }}_0},0⩽k⩽n-1$$

• 定义3：设 { ε t , t = 0 , ± 1 , ± 2 , ⋯   } \left\{\varepsilon_t,t=0,\pm1,\pm2,\cdots\right\} {εt​,t=0,±1,±2,⋯}是零均值平稳白噪声，$\mathrm{Var}({\epsilon }_t)={\sigma }_{\epsilon }^2$，若 { G k   , k = 0 , 1 , 2 , ⋯   } \left\{G_k\ ,k=0,1,2,\cdots\right\} {Gk​ ,k=0,1,2,⋯}是一数列，满足$$\sum _{k=0}^{\infty }|G_k|<+\infty ,G_0=1$$定义随机序列$$X_t=\sum _{k=1}^{\infty }{G}_k{\epsilon }_{t-k}$$则称$X_t$为随机线性序列，可证明$X_t$是平稳序列。

• 定义4：设 { X t , t = 0 , ± 1 , ± 2 , ⋯   } \left\{X_t,t=0,\pm1,\pm2,\cdots\right\} {Xt​,t=0,±1,±2,⋯}是零均值平稳序列，如果已知 { X t − 1 , X t − 2 , ⋯   , X t − k } \left\{X_{t-1},X_{t-2},\cdots,X_{t-k}\right\} {Xt−1​,Xt−2​,⋯,Xt−k​}的值，要对$X_t$作出预报。考虑线性最小均方估计，选择系数${\phi }_{k,1},{\phi }_{k,2},\cdots ,{\phi }_{k,k}$，使得$$\min \delta =E[(X_t-\sum _{j=1}^k{\phi }_{k,j}{X}_{t-j}{)}^2]$$由$\frac{\partial \delta }{\partial {\phi }_{k,j}}=0$，可解得${\phi }_{k,j}$应满足的线性方程组$$\left[\begin{array}{cccc}1& {\rho }_1& \cdots & {\rho }_{k-1}\\ {\rho }_1& 1& \cdots & {\rho }_{k-2}\\ ⋮& ⋮& \ddots & ⋮\\ {\rho }_{k-1}& {\rho }_{k-2}& \cdots & 1\end{array}\right]\left[\begin{array}{c}{\phi }_{k,1}\\ {\phi }_{k,2}\\ ⋮\\ {\phi }_{k,k}\end{array}\right]=\left[\begin{array}{c}{\rho }_1\\ {\rho }_2\\ ⋮\\ {\rho }_k\end{array}\right]$$该式称为Yule-Walker方程， { φ k , k   , k = 1 , 2 , ⋯   } \left\{\varphi_{k,k}\ ,k=1,2,\cdots\right\} {φk,k​ ,k=1,2,⋯}称为$X_t$的偏相关系数。

2.2 ARMA时间序列
ARMA时间序列是一种重要的平稳时间序列，分为三种类型：AR序列（自回归序列）、MA序列（移动平均序列）、ARMA序列（自回归移动平均序列）。

2.2.1 AR（p）序列
设 { X t , t = 0 , ± 1 , ± 2 , ⋯   } \left\{X_t,t=0,\pm1,\pm2,\cdots\right\} {Xt​,t=0,±1,±2,⋯}是零均值平稳序列，满足下列模型$$X_t={\phi }_1{X}_{t-1}+{\phi }_2{X}_{t-2}+\cdots +{\phi }_p{X}_{t-p}+{\epsilon }_t$$式中${\epsilon }_t$为零均值、方差是${\sigma }_{\epsilon }^2$的平稳白噪声；$X_t$为阶数为p的自回归序列，简记为AR（p）序列；$\mathbf{\phi }=[{\phi }_1,{\phi }_2,\cdots ,{\phi }_p{]}^T$为自回归参数向量，其分量${\phi }_i$称为自回归系数。
引入后移算子B，定义为$$B{X}_t\equiv X_{t-1},B^k{X}_t\equiv X_{t-k}$$记算子多项式$$\phi (B)=1-{\phi }_{1}B-{\phi }_2{B}^2-\cdots -{\phi }_p{B}^p$$则模型公式可改写为$$\phi (B)X_t={\epsilon }_t$$

2.2.2 MA（q）序列
设 { X t , t = 0 , ± 1 , ± 2 , ⋯   } \left\{X_t,t=0,\pm1,\pm2,\cdots\right\} {Xt​,t=0,±1,±2,⋯}是零均值平稳序列，满足下列模型$$X_t={\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-{\theta }_2{\epsilon }_{t-2}-\cdots -{\theta }_q{\epsilon }_{t-q}$$式中${\epsilon }_t$为零均值、方差是${\sigma }_{\epsilon }^2$的平稳白噪声；$X_t$为阶数为q的移动平均序列，简记为MA（q）序列；$\mathbf{\theta }=[{\theta }_1,{\theta }_2,\cdots ,{\theta }_q{]}^T$为移动平均参数向量，其分量${\theta }_i$称为移动平均系数。
对线性后移算子B，有$$B{\epsilon }_t\equiv {\epsilon }_{t-1},B^k{\epsilon }_t\equiv {\epsilon }_{t-k}$$记算子多项式$$\theta (B)=1-{\theta }_{1}B-{\theta }_2{B}^2-\cdots -{\theta }_q{B}^q$$则模型公式可改写为$$X_t=\theta (B){\epsilon }_t$$

2.2.3 ARMA（p,q）序列
设 { X t , t = 0 , ± 1 , ± 2 , ⋯   } \left\{X_t,t=0,\pm1,\pm2,\cdots\right\} {Xt​,t=0,±1,±2,⋯}是零均值平稳序列，满足下列模型$$X_t-{\phi }_1{X}_{t-1}-\cdots -{\phi }_p{X}_{t-p}={\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-\cdots -{\theta }_q{\epsilon }_{t-q}$$式中${\epsilon }_t$为零均值、方差是${\sigma }_{\epsilon }^2$的平稳白噪声；$X_t$为阶数为p,q的自回归移动平均序列，简记为ARMA（p,q）序列。当q=0时是AR（p）序列，当p=0时是MA（q）序列。
应用算子多项式$\phi (B),\theta (B)$，模型公式可改写为$$\phi (B)X_t=\theta (B){\epsilon }_t$$对一般的平稳序列 { X t , t = 0 , ± 1 , ± 2 , ⋯   } \left\{X_t,t=0,\pm1,\pm2,\cdots\right\} {Xt​,t=0,±1,±2,⋯}，设其均值$E(X_t)=\mu$，满足下列模型$$(X_t-\mu )-{\phi }_1(X_{t-1}-\mu )-\cdots -{\phi }_p(X_{t-p}-\mu )={\epsilon }_t-{\theta }_1{\epsilon }_{t-1}-\cdots -{\theta }_q{\epsilon }_{t-q}$$应用算子多项式$\phi (B),\theta (B)$，模型公式可改写为$$\phi (B)(X_t-\mu )=\theta (B){\epsilon }_t$$关于算子多项式$\phi (B),\theta (B)$，通常还要作下列假定：

• $\phi (B)$和$\theta (B)$无公共因子，又${\phi }_p\ne 0,{\theta }_q\ne 0$
• $\phi (B)=0$的根全在单位圆外，这一条件称为模型的平稳性条件
• $\theta (B)=0$的根全在单位圆外，这一条件称为模型的可逆性条件

2.3 ARMA模型的构建及预报
在实际问题建模中，首先要进行模型的识别与定阶，即要判断AR（p）、MA（q）、ARMA（p,q）模型的类别，并估计阶数p,q。当模型定阶后，就要对模型参数$\mathbf{\phi }=[{\phi }_1,{\phi }_2,\cdots ,{\phi }_p{]}^T$和$\mathbf{\theta }=[{\theta }_1,{\theta }_2,\cdots ,{\theta }_q{]}^T$进行估计。定阶与参数估计完成后，还要对模型进行检验，即要检验${\epsilon }_t$是否为平稳白噪声。若检验通过，则ARMA时间序列的建模完成。

2.3.1 ARMA模型的构建

• ARMA模型定阶的AIC准则：选p,q，使得$$\min \mathrm{AIC}=n\ln {\hat{\sigma }}_{\epsilon }^2+2(p+q+1)$$其中n为样本容量；${\hat{\sigma }}_{\epsilon }^2$为${\sigma }_{\epsilon }^2$的估计，与p和q有关。当ARMA（p,q）序列含有未知均值参数$\mu$时，此时未知参数个数为k=p+q+2，AIC准则为：选p,q，使得$$\min \mathrm{AIC}=n\ln {\hat{\sigma }}_{\epsilon }^2+2(p+q+2)$$实际上，这两个公式有相同的最小值点$\hat{p},\hat{q}$。
• ARMA模型的参数估计：有矩估计、逆函数估计、最小二乘估计、条件最小二乘估计、最大似然估计等方法。可直接使用Matlab工具箱给出相关的参数估计。
• ARMA模型的${\chi }^2$检验：Matlab中作Ljung-Box检验的函数为lbqtest

2.3.2 ARMA（p,q）序列的预报
时间序列的m步预报，是根据 { X k , X k − 1 , ⋯ } \left\{\bm{X_k,X_{k-1},\cdots}\right\} {Xk​,Xk−1​,⋯}的取值对未来k+m时刻的随机变量${\mathbf{X}}_{\mathbf{k}+\mathbf{m}}(m>0)$作出估计。估计量记作${\hat{X}}_k(m)$，它是${\mathbf{X}}_{\mathbf{k}},{\mathbf{X}}_{\mathbf{k}-1},\cdots$的线性组合。

• AR（p）序列的预报：递推公式为$$\{\begin{array}{l}{\hat{X}}_k(1)={\phi }_1{X}_k+{\phi }_2{X}_{k-1}+\cdots +{\phi }_p{X}_{k-p+1}\\ {\hat{X}}_k(2)={\phi }_1{\hat{X}}_k(1)+{\phi }_2{X}_k+\cdots +{\phi }_p{X}_{k-p+2}\\ ⋮\\ {\hat{X}}_k(p)={\phi }_1{\hat{X}}_k(p-1)+{\phi }_2{\hat{X}}_k(p-2)+\cdots +{\phi }_{p-1}{\hat{X}}_k(1)+{\phi }_p{X}_k\\ {\hat{X}}_k(m)={\phi }_1{\hat{X}}_k(m-1)+{\phi }_2{\hat{X}}_k(m-2)+\cdots +{\phi }_p{\hat{X}}_k(m-p),m>p\end{array}$$由此可见，${\hat{X}}_k(m)(m>1)$仅依赖于$X_t$的k时刻以前的p个时刻的值。这是AR（p）序列预报的特点。
• MA（q）与ARMA（p,q）序列的预报：
  对于MA（q）序列的预报，有$${\hat{X}}_k(m)=0,m>q$$因此，只需要讨论${\hat{X}}_k(m),m=1,2,\cdots ,q$。为此，定义预报向量$${\hat{\mathbf{X}}}_k^{(q)}=[{\hat{X}}_k(1),{\hat{X}}_k(2),\cdots ,{\hat{X}}_k(q){]}^T$$递推预报是求${\hat{X}}_k^{(q)}$和${\hat{X}}_{k+1}^{(q)}$的，递推关系，对MA（q）序列，有$${\hat{\mathbf{X}}}_{k+1}^{(q)}=\left[\begin{array}{ccccc}{\theta }_1& 1& 0& \cdots & 0\\ {\theta }_2& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\theta }_{q-1}& 0& 0& \cdots & 1\\ {\theta }_q& 0& 0& 0& 0\end{array}\right]{\hat{\mathbf{X}}}_k^{(q)}-\left[\begin{array}{c}{\theta }_1\\ {\theta }_2\\ ⋮\\ {\theta }_q\end{array}\right]{\mathbf{X}}_{k+1}$$递推初值可取${\hat{\mathbf{X}}}_{k_0}^{(q)}=0(k_0较小)$。因为模型的可逆性保证了递推式渐进稳定，即当n充分大后，初始误差的影响可以逐渐消失。
  对于ARMA（p,q）序列的预报，有$${\hat{\mathbf{X}}}_{k+1}^{(q)}=\left[\begin{array}{ccccc}-G_1& 1& 0& \cdots & 0\\ -G_2& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -G_{q-1}& 0& 0& \cdots & 1\\ -G_q+{\phi }_q^{\ast }& {\phi }_{q-1}^{\ast }& {\phi }_{q-2}^{\ast }& \cdots & {\phi }_1^{\ast }\end{array}\right]{\hat{\mathbf{X}}}_k^{(q)}+\left[\begin{array}{c}{G}_1\\ G_2\\ ⋮\\ G_{q-1}\\ G_q\end{array}\right]{\mathbf{X}}_{k+1}+\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\\ {\sum }_{j=q+1}^p{\phi }_j^{\ast }{X}_{k+q+1-j}\end{array}\right]$$其中${\phi }_j^{\ast }=\{\begin{array}{l}{\phi }_j,j=1,2,\cdots ,p\\ 0,j>p\end{array}$；$G_j$满足${\mathbf{X}}_t={\sum }_{j=0}^{\infty }{G}_j{\epsilon }_{t-j}$；式中第三项当$p⩽q$时为0。递推初值可取${\hat{\mathbf{X}}}_{k_0}^{(q)}=0(k_0较小)$。

3. 时间序列的Matlab工具箱及命令

Matlab时间序列相关命令对应的工具箱有：系统辨识工具箱、计量经济学工具箱、金融工具箱。

3.1 系统辨识工具箱相关命令

• 观察数据的获取命令：idinput
• 数据预处理命令：detrend、idfilt、idresamp
• 模型结构的选择命令：struc、arxstruc、ivstruc、selstruc
• 参数估计命令：ar、arx、armax、ivx
• 模型预报与仿真的命令：compare、pe、predict、sim

%%%%%%%%%%%%%%%%%%%拟合AR(p)模型%%%%%%%%%%%%%%%%%%%%%%
for i=1:3
    m{i}=ar(x,i);		%x为时间序列数据的行向量，分别以AR(1)、AR(2)、AR(3)拟合模型
    myaic(i)=aic(m{i});	%计算AIC的值，选取最小的值作为AR(p)的p值
end
xp=predict(m,x);	%计算已知数据的预测值 
res=x-xp;		%计算残差向量，也可以使用命令res=resid(m,x)计算残差向量
h=lbqtest(res)		%对残差向量进行Ljung-Box检验
%%%%%%%%%%%%%%%%%%%拟合ARMA(p,q)模型%%%%%%%%%%%%%%%%%%%%%%
for i=0:3
    for j=0:3
        if i==0 & j==0
            continue	%ARMA(p,q)模型中，p,q不能同时为0
        end
    m=armax(x,[i,j]);	%拟合模型ARMA(i,j)，已知数据x必须是列向量
    myaic=aic(m);		%计算AIC指标
    fprintf('p=%d,q=%d,AIC=%f\n',i,j,myaic);	%显示计算结果
    end
end

3.2 计量经济学工具箱相关命令

• 模型参数获取或设置命令：garchget、garchset
• 建模和仿真等命令：garchfit、garchpred、garchsim

GARCH模型：GARCH表示广义自回归条件异方差。该模型分为均值方差和方差方差两部分。

• 均值方程形式：$$y_t=C+\sum _{i=1}^R{\phi }_i{y}_{t-i}+{\epsilon }_t+\sum _{j=1}^M{\theta }_j{\epsilon }_{t-j}+\sum _{k=1}^{N_x}{\beta }_k\mathbf{X}(t,k)$$式中${\phi }_i$为自回归系数；${\theta }_j$为移动平均系数；$\mathbf{X}$为解释回归矩阵，它的每一列是一个时间序列。
• 方差方程形式：$${\sigma }_t^2=K+\sum _{i=1}^P{G}_i{\sigma }_{t-i}^2+\sum _{j=1}^Q{A}_j{\epsilon }_{t-j}^2$$其系数满足以下约束条件：${\sum }_{i=1}^P{G}_i+{\sum }_{j=1}^Q{A}_j<1,K>0,G_i>0,A_j>0$。

model = garch(P,Q)		%根据P,Q值构建参数大小未知的garch模型
model = garch(Name,Vaule)		%根据属性、属性值构建参数确定的garch模型
%%%%%%%%%%%%%%%%通过给定参数构建一个garch模型，然后再用garch(1,1)模型预测参数%%%%%%%%%%%%%%
VarMd=garch('Constant',0.01,'GARCH',0.2,'ARCH',0.3);%指定模型方差方程的参数，当参数为向量时用用元胞数组作为参数值输入
Md=arima('Constant',0,'AR',0.8,'MA',0.4, 'Variance',VarMd);	%指定模型均值方差的参数
[y,e,v]=simulate(Md,10000);		%产生指定结构模型的10000个模拟数据
ToEstVarMd=garch(1,1);
ToEstMd=arima('ARLags',1,'MALags',1,'Constant',0,'Variance',ToEstVarMd);%指定ARMA模型的p,q值	，高阶时输入数组
[EstMd,EstParamCov,logL,info]=estimate(ToEstMd,y)	%模型拟合
	%EstMd为拟合得到的模型，EstParamCov为估计参数的协方差矩阵，logL为似然目标函数，info输入参数信息
res=infer(EstMd,y);		%计算残差
h=lbqtest(res)		%进行模型检验
yhat=forecast(EstMd,3,'Y0',y)		%预测未来的3个值

3.3 金融工具箱相关命令

• 构造时间序列数组的命令：fints、ascii2fts
• 用户图形界面解法命令：ftsgui、ftstool

4. ARIMA序列与季节性序列

在实际中遇到的时间序列往往具有趋势性、季节性、非平稳性。可使用差分方法消除其趋势性、季节性，使变换后的序列是平稳序列，并假设为ARMA序列，再用上面介绍的方法进行建模和预报。

4.1 ARIMA序列与预报
4.1.1 ARIMA序列模型
对于某一时间序列$X_t$，若其$\phi (B)=0$的根并非全在单位圆外（有$B⩽1$的根），则该序列为非平稳序列。此时进行一阶向后差分运算$$\nabla X_t=X_t-X_{t-1}=W_t$$可将原来的非平稳序列$X_t$转化为非平稳序列$W_t$。如果一阶差分还不能使时间序列平稳化，还可以进行更高阶差分，最后将序列化为平稳序列。
对d阶差分，有$${\nabla }^d{X}_t=(1-B{)}^d{X}_t$$式中${\nabla }^d$称为d阶差分算子。设$X_t$是非平稳序列，存在正整数d，使得$${\nabla }^d{X}_t=W_t$$而$W_t$是ARMA（p,q）序列，则称$X_t$是ARIMA（p,d,q）序列。此时，$X_t$满足$$\phi (B){\nabla }^d{X}_t=\theta (B){\epsilon }_t或\phi (B)({\nabla }^d{X}_t-\mu )=\theta (B){\epsilon }_t,t>d$$后者称$X_t$是一般ARIMA（p,d,q）序列。
若$X_t$的观测样本是$X_1,X_2,\cdots ,X_n$，经过d阶差分后，数据为n-d个。由d阶差分${\nabla }^d{X}_t$复原数据，需要给定初值$X_1,X_2,\cdots ,X_d$。

4.1.2 在确定模型时，常用以下方法：

• 先对$X_t$的样本$X_1,X_2,\cdots ,X_n$计算样本自相关函数和样本偏相关函数。
• 如果是截尾或拖尾，说明已服从ARMA模型。
• 若自相关函数和偏相关函数至少有一个不是截尾或拖尾的，说明$X_t$不是平稳的。
• 对非平稳序列作一阶差分$\nabla X_t,t=2,3,\cdots ,n$，并求其样本自相关函数和样本偏相关函数，再用上述方法讨论。
• 如此进行，直至判断${\nabla }^d{X}_t$是平稳序列为止。

截尾：在大于某个常数k后快速趋于0
拖尾：在大于某个常数k后在0附近随机波动

若初值$X_1,X_2,\cdots ,X_d$已知，由$W_t={\nabla }^d{X}_t,t=d+1,\cdots ,n$可以复原$X_t$，d=1,2时最常见。

• d=1时$$X_t=X_1+\sum _{j=1}^{t-1}{W}_{j+1}=X_k+\sum _{j=1}^{t-k}{W}_{j+k},t>k⩾1$$
• d=2时$$\begin{gathered} X_t=X_2+(t-2)(X_2-X_1)+\sum _{j=1}^{t-2}(t-j-1)W_{j+2} \\ =X_k+(t-k)(X_k-X_{k-1})+\sum _{j=1}^{t-k}(t-j-k+1)W_{j+k},t>k⩾2 \end{gathered}$$

4.1.3 ARIMA序列的预报
设$X_t$是ARIMA（p,d,q）序列，d=1和d=2的情形在实用中使最常见的。

• 当d=1时，有$\nabla {\hat{X}}_k(m)={\hat{W}}_k(m)$，预报公式为$${\hat{X}}_k(m)=X_k+\sum _{j=1}^m{\hat{W}}_k(j)$$
• 当d=2时，有${\nabla }^2{\hat{X}}_k(m)={\hat{W}}_k(m)$，预报公式为$${\hat{X}}_k(m)=X_k+m(X_k-X_{k-1})+\sum _{j=1}^m(m+1-j){\hat{W}}_k(j)$$

r11=autocorr(a)		%计算自相关函数，a为列向量，默认计算20个滞后时间
r12=parcorr(a)		%计算偏相关函数
[EstMd,EstParamCov,logL,info]=estimate(ToEstMd,da);	%模型拟合
numParams = sum(any(EstParamCov));		%计算拟合参数的个数
[aic(k),bic(k)]=aicbic(logL,numParams,n);	%计算AIC和BIC信息

4.2 季节性序列及其预报
由于季节性因素或其他周期因素引起的周期性变化的时间序列，称为季节性时间序列，相应的模型为季节性模型。
一般地，对周期为s的序列，可先进行差分运算，在进行ARIMA建模，即$$\begin{gathered} {\nabla }_s{X}_t=(1-B^s)X_t \\ {\nabla }_s^d{X}_t=(1-B^s{)}^d{X}_t \end{gathered}$$