计量经济学模型及R语言应用

4 经典时间序列模型

4.1 时间序列的基本概念

4.1.1 时间序列的含义

　　从统计意义上讲，所谓时间序列就是将某一个指标在不同时间上的不同数值，按照时间的先后顺序排列而成的数列。这些时间序列由于受到各种偶然因素的影响，往往表现出随机性，彼此之间存在着统计上的相互依赖关系，而时间序列中的“时间”可以具有不同的物理意义,例如长度、温度、速度等等。
　　从数学意义上讲，如果我们对某一过程中的某一个变量或一组变量X(t)进行观察测量，在一系列时刻$t_1$,$t_2$,$t_3$，…，$t_4$（t为自变量，且$t_1$< $t_2$< $t_4$<…< $X_n$）得到的离散有序数集合$X_{t1}$，$X_{t2}$，$X_{t3}$，... ，$X_{tn}$,... ,$X_{tn}$称为离散数字时间序列，即随机过程的一次样本实现。设X($t:t{\in}T$)是一个随机过程，$X_{ti}$（i=1,2,…）是在时刻i对过程X(t)的观测值，则$X_{ti}$（i=1,2,…）称为一次样本实现，也就是一个时间序列。
　　从系统意义上看，时间序列就是某一系统在不同时间（地点、条件等）的响应。这个定义从系统运行的观点出发，不仅指出时间序列是按一定顺序排列而成的；这里的“一定顺序”既可以是时间顺序，也可以是具有各种不同意义的物理量，如代表温度，速度或其它单调递增地取值的物理量。可见，时间序列只强调顺序的重要性，而并非强调必须以时间顺序排列。
　　由于时间序列是所研究系统的历史行为的客观记录，因而它包含了系统结构特征及其运行规律。因此通过对时间序列的研究，可以认识所研究系统的结构特征（如周期波动的周期、振幅和趋势的种类等）；揭示其运行规律，进而用以预测、控制其未来行为；修正和重新设计系统（如改变其周期、参数），使之按照新的结构运行。
　　时间序列根据所研究的依据不同，有不同的分类。
　　1. 按研究的对象的多少分，有单变量时间序列和多变量时间序列。如果所研究的对象是一个变量，如某个国家的国内生产总值，即为单变量时间序列。如果所研究的对象是多个变量，如按年、月顺序排列的气温、气压、雨量数据，则为多变量时间序列。多变量时间序列不仅描述了各个变量的变化规律，而且还揭示了各变量间相互依存关系的动态规律性。
　　2. 按时间的连续性可将时间序列分为离散时间序列和连续时间序列。如果某一序列中的每一个序列值所对应的时间参数为间断点，则该序列就是一个离散时间序列。如果某一序列中的每个序列值所对应的时间参数为连续函数，则该序列就是一个连续时间序列。
　　3. 按序列的统计特性分，有平稳时间序列和非平稳时间序列两类。如果一个时间序列的概率分布与时间t无关，则称该序列为严格的（狭义的）平稳时间序列。如果序列的一、二阶矩存在，而且对任意时刻t满足：
　　（1）均值为常数
　　（2）协方差为时间间隔 的函数
　　则称该序列为宽平稳时间序列，也叫广义平稳时间序列。反之，不具有平稳性即序列均值或协方差与时间有关的序列称为非平稳序列。
　　一个平稳的随机时间序列过程，在直观上可以看作一条围绕其均值上下波动的曲线，在理论上，我们把具有下列性质的随机过程称为平稳随机时间序列：
　　(1) $E(y_t) =\mu$ (对所有t)
　　(2) V($y_t$) = ${\sigma}_y^{2}$= 常数 (对所有t)
　　(3) $COV(y_t, y_t+k)$ 仅与$y_t$和$y_t+k$相隔的时期数有关，而与时间点t无关(对所有t和k)。
　　随机过程是否具备平稳性对于时间序列预测来说十分重要，这一性质保证了随机过程的结构不会随时间变化，这是进行准确预测的必要条件。
　　4. 按序列的分布规律来分，有高斯型（Gaussian）时间序列和非高斯型（non-Gaussian）时间序列。服从高斯分布（正态分布）的时间序列叫做高斯型时间序列，否则叫做非高斯型时间序列。对于一些非高斯序列，往往可以通过适当的变换，可近似地看成是高斯型时间序列。
　　本章我们主要研究平稳时间序列的一些建模技术。
通俗来说，所谓时间序列的平稳性，是指时间序列的统计规律不会随着时间的推移而发生变化。直观上，一个平稳的时间序列可以看作一条围绕其均值上下波动的曲线。 最简单的随机时间序列yt是一个具有零均值同方差的独立同分布序列：

$$y_t = u_t　　{u_t}\sim N(0, {\sigma}^{2})$$
　　该序列常被称为是一个白噪声(whitenoise)。由于yt是一个具有相同的均值和方差序列，由平稳性的定义知白噪声序列是平稳的。
　　【例4.1】模拟数据---平稳时间序列。
　　设${u_t}\sim N(\mu, {\sigma}^{2})$，令$\mu= 0,{\sigma}^{2}=1^2$，下面将产生50个标准的白噪声序列。

set.seed(123)         
yt=rnorm(50,0,1)
plot(yt,type='b'); abline(h=0)

4.1.2 时间序列的相关性

　　这里所说的时间序列的相关性跟前面所讲的误差序列的相关性有所不同，这里是考证原始数据序列的自相关性，而前面讲的是残差序列的自相关性，但处理方法基本相同。
一、自相关函数
　　自相关函数是衡量序列yt中任意两个元素之间相关程度的量度。对随机过程${ y_t }$，元素$y_t$与$y_t+k$之间的自相关函数定义如下：

$${\rho}_{k}=\frac{cov(y_t,y_{t+k})}{\sqrt{var{y_t}var{y_t}}}（4.1-1）$$
　　自相关系数${\rho}_{k}$的序列{${\rho}_{k}$ }(k = 0，±1，±2，…)称为自相关函数（ACF）。
　　当$y_t$为平稳随机过程时
　　　　　　　　　　　 $${\rho}_{k}=\frac{\gamma _k}{\sigma _{y}^2}=\frac{\gamma _k}{\sigma _{0}}　　（4.1-2）$$ 　　　
　　其中$\sigma _0=V(y_t)=\sigma _{y}^2 　\gamma _k=cov(y_t,y_{t+k})$
　　由定义知，对任何随机过程${\rho}_{0}=1$，由公式知，${\rho}_{k}$是一个无量纲量。
　　在实际计算时，我们只能计算样本自相关函数，其样本自相关函数定义为
$$\hat{\rho _k}=\frac{\sum_{i=1}^{n-k}(y_t-\bar{y})(y_{t+k}-\bar{y})}{\sum_{i=1}^{n}(y_t-\bar{y})^2}　（4.1-3）$$
　　自相关度量了yt与其滞后yt+k之间的相关性。
　　在例4.1-1中，我们模拟了一组无自相关性的数据$y_t$，下面用我们定义的滞后函数L_函数来考证其相关性。

yt_1=L_(yt,na.is=T)             ＃一阶滞后
plot(yt,yt_1); abline(lm(yt~yt_1)) ＃相关图 

cor(yt,yt_1,"complete")         ＃相关系数

[1] 0.02567

　　从图上可以看出，$y_t$与其一阶滞后$y_{t-1}$的相关性很低，其相关系数（不是一阶自相关系数）为0.02567。
二、自相关系数的计算
　　要计算序列的自相关性，可用R语言自带的计算自相关函数（auto correlation function，简称ACF）acf，自相关函数的对象ac为自相关系数。

acf(yt)＄ac
>, , 1
          [,1]
 [1,]  1.00000
 [2,]  0.02558
 [3,] -0.17083
 [4,]  0.22624
 [5,] -0.10638
 [6,]  0.04813
 [7,] -0.03072
 [8,] -0.15668
 [9,]  0.06286
[10,] -0.13276
[11,] -0.08354
[12,]  0.06820
[13,] -0.26488
[14,] -0.07688
[15,]  0.16640
[16,] -0.19327
[17,] -0.01659

　　注意由简单相关$cor(y_t,y_{t_1},"complete")$算出的相关系数0.02567和由自相关函数计算的一阶自相关系数0.02558还是有些不同的（公式略有不同）。
　　同理由$cor(y_t,y_{t_1},"complete")$算出的二阶相关系数-0.1727和由自相关函数计算的二阶自相关系数-0.1708也是不同的。

4.1.3 序列自相关性判别

一、图示法
　　随机时间序列模型着重研究的是相关关系，因此自相关函数在时间序列模型中占有重要地位。但注意的是：在实际识别时，样本自相关函数$\hat{\rho}_k$ 只是总体自相关函数$\rho_k$的一个估计，由于样本的随机性，当k > q时，$\hat{\rho}_k$不会全为0，而是在0的上下波动。但可以证明，当k > q时，$\hat{\rho}_k$服从如下渐近正态分布: $\hat{\rho}_k ～N(0,1/n)$, 式中n表示样本容量。
　　因此，如果计算的$\hat{\rho}_k$ 满足：$\begin{vmatrix} \hat{\rho}_k \end{vmatrix} <\frac{1.96}{\sqrt{n}}$
　　我们就有95%的把握判断原时间序列不存在自相关性。
　　下面我们采用R语言自带的自相关函数来分析数据的自相关性。

acf(yt)

　
　　从自相关图上可以看出，所有自相关系数都落入可信区间之内，所以可以认为序列yt是一个不存在自相关的时间序列（符合我们的模拟）。

二、假设检验
　　序列自相关检验，即检验序列是否是白噪声。如果序列不是平稳过程，说明有些有用的动态结构没有进入时间序列模型，而这些结构对于描述时间序列往往是必不可少的。检验序列相关性的方法有很多，如3.3.2节的DW检验和LM检验。但DW检验有一个缺点，就是其只能检验一阶自相关，下面给出两种常用的检验序列各阶自相关的方法。
　　1. Box- Pierce检验
　　博克斯(G.E.P.Box)和皮尔斯(D.A.Pierce)提出的Q统计量可以检验时间序列的相关性。Q统计量定义为：

$$Q=n\sum_{k=1}^{m}\hat{\rho}_k^2　　（4.1-4）$$
　　其中n为样本容量，m为滞后长度。在大样本的情况下，它近似服从自由度为m的$\gamma^2$分布。若计算出来的Q大于在选定显著性水平下从的$\gamma^2$分布表中查出的临界Q值，则拒绝所有真实的($r_k$都为0)虚拟假设，这时序列不存在自相关；否则，序列存在自相关。
　　2. Ljung-Box检验
　　1978年扬-博克斯将博克斯(G.E.P.Box)和皮尔斯(D.A.Pierce)的Q统计量变形为LB统计量(Ljung-Box Statistic),其定义为：
$$LB=n(n+2)\sum_{k=1}^{m}(\frac{\hat{\rho}_k^2}{n-k})～\chi ^2(m)　　（4.1-5）$$
　　其检验过程与Q统计量检验过程一样，但LB统计量比Q统计量有更好的小样本性质(即在统计意义上更有效)，所以LB统计量常用来检验小样本序列的相关性。在利用R进行序列的相关性检验时，可以利用R作时间序列自相关图，直接分析序列的ACF，如果ACF都不显著，那么说明其序列不相关。R语言中stats程序包里的Box.test函数可以实现上述两个检验。

Box.test(yt)     ＃Box.test(yt,type='Box-Pierce')
    Box-Pierce test
data:  yt
X-squared = 0.0327, df = 1, p-value = 0.8564

Box.test(yt,type='Ljung-Box')

Box-Ljung test

data: yt
X-squared = 0.0347, df = 1, p-value = 0.8522

　　Box-Pierce和Ljung-Box方法都检验序列$y_t$没有自相关性（p>0.05)，跟模拟结果一致。 
　　关于时间序列平稳性的检验请看5.1.3节的单位根检验。
##4.2 时间序列自回归AR模型
　　建模的一个主要目的是为了预测，那么，随机过程的元素$y_t$对它的过去的依赖性就很重要。这使我们能够利用已经收集的样本观测值的过去信息预测变量的未来值。存在这种依赖性的简单例子是自回归（Auto Regression，简称AR）过程：
$$y_t=\phi_1y_{t-1}+\mu_t　　(4.2-1)$$                           
便是这样一种一阶自回归模型，其中ut为白噪声。
　　自回归模型描述序列{$y＿t$}在某一时刻t和前p时刻序列值之间的线性关系，表示为：
$$y_t=\phi_1y_{t-1}+\phi_2y_{t-2}+...+\phi_py_{t-p}+\mu_t 　　(4.2-2)$$
　　这里随机序列${u_t}$是白噪声，服从均值为0，方差为$\sigma_e^2$ 的正态分布, 且${\mu_t}$与序列${y_t}(k<t)$不相关，称模型(4.2-2)是p阶自回归模型，记为AR(p)。实参数$\phi_2,\phi_2,...,\phi_p$称为自回归系数，是模型的待估参数。
###4.2.1 AR模型的平稳性条件
　　只有产生时间序列的随机过程是平稳的，用自回归模型进行预测才有意义。因此，我们首先应研究自回归过程的平稳条件。
一、一阶自回归过程
　　对于一阶自回归过程(4.2-1)
　　
$$\begin{align*}
y&=\phi{y_{t-1}}+\mu_t\\ 
&=\mu_t+\phi(\phi{y_{t-2}}+\mu_{t-1})\\ 
&=\mu_t+\phi\mu_{t-1}+\phi^2(\phi{y_{t-3}}+\mu_{t-2})\\ 
&=\mu_t+\phi\mu_{t-1}+\phi^2\mu_{t-2}+{\phi^3}y_{t-3}\\ 
&...  　...\\ 
&=\mu_t+\phi\mu_{t-1}+\phi^2\mu_{t-2}+{\phi^3}\mu_{t-3}+...　　(4.2-3)
\end{align*}$$
　　可以看到，一阶自回归过程(4.2-1)表示成白噪声序列的线性组合。
　　由于$E(u_t)= 0$，所以$E(y_t)=0$平稳条件1显然满足。对(4.2-3)两端取方差：
$$Var(y_t)=\sigma_{\mu}^2(1+\phi^2+\phi^4+\phi^6+...)　　  (4.2-4)$$
　　仅当$|\phi| < 1$时，(4.2-4)才有
$$V(y_t)=\frac{\sigma_{\mu}^2}{1-\phi^2}　　(4.2-5)$$            
　　表明，只有当$|\phi| < 1$时，平稳条件2才成立。
　　由(4.2-3)有
$$y_{t+k}=\mu_{t+k}+\phi\mu_{t+k-1}+\phi^2\mu_{t+k-2}+\phi^k\mu_t+...　　(4.2-6)$$
$$\begin{align*}
COV(y_t,y_{t+k}&=E(y_ty_{t+k})=\phi^k\sigma_{\mu}^2+\phi^{k+2}\sigma_{\mu}^2+\phi^{k+4}\sigma_{\mu}^2+...\\
&=\sigma_{\mu}^2\phi^k(1+\phi^2+\phi^4+...)　　(4.2-7)
\end{align*}$$
　　当$|\phi| < 1$时，(4.2-7)便有
$$COV(y_t,y_{t+k}=\phi^k\frac{\sigma_{\mu}^2}{1-\phi^2}=\phi^k\sigma_y^2　　(4.2-8)$$
　　其中$\sigma_y^2=V(y_t)$ 。
　　(4.2-8)式表明仅与间隔时期数k有关，而与时间点t无关，平稳条件3成立。
　　综上所述，对于一阶自回归过程(4.2-1)，只要系数 的绝对值$|\phi| < 1$，便是平稳过程。
　　由ACF的性质可知，对于AR(p)模型，其ACF不能在某一步后为零(截尾)，而是按指数率衰减，称其具有拖尾性。
　　【例4.2-1】模拟研究。模拟AR(1)模型$y_t=0.8y_{t+1}+e_t,e_t～N(0,1)$

set.seed(1234)
n=50;
y1=rep(0,n); for(t in 2:n) { #也可用arima.sim(n=50,list(ar=0.8))，见下例
y1[t]=0.8*y1[t-1]+rnorm(1)
}
plot(y1,type='o') #散点图

![图片标题](https://leanote.com/api/file/getImage?fileId=5839a7ffab64413781014973)

acf(y1) ＃自相关图

![图片标题](https://leanote.com/api/file/getImage?fileId=5839a82fab6441366c014426)

Box.test(y1)

> Box-Pierce test
>data:  y1
X-squared = 28.78, df = 1, p-value = 8.088e-08

Box.test(y1,type='Ljung-Box')

> Box-Ljung test
>data:  y1
X-squared = 30.55, df = 1, p-value = 3.259e-08
　　各种检验都表明，序列存在自相关性。由于$\phi= 0.8 < 1$，所有说该序列应该是平稳的，这点也可以从其散点图上看出。
二、p阶自回归过程
　　将(4.2-2)改写成
    $$(1-\phi_1L-\phi_2L^2-...-\phi_pL^p)y_t=\mu_t　　(4.2-9)$$
　　引进算符多项式：
      $$\phi_p(L)=1-\phi_1L-\phi_2L^2-...-\phi_pL^p 　　(4.2-10)$$
　　则(4.2-9)可改写成：
  $$\phi_p(L)y_t=\mu_t或y_t=\phi_p^{-1}(L)y_t=\mu_t　　 (4.2-11)$$
　　这里L就是我们前面说的滞后算子，及$Ly_ty_{t-1}$等，在R语言中写为yt_1=L_(yt,1)，详见3.3.3节。
　　若(4.2-2)是平稳随机过程，则$\phi_p^{-1}(L)$必定收敛，即可表示为白噪声的无穷加权和。可以证明， $\phi_p^{-1}(L)$收敛的充要条件是算符多项式$\phi_p(L)$ 的特征方程
 $$\phi_p(z)=1-\phi_1z-\phi_2z^2-\phi_3z^3-...-\phi_pz^p=0　　(4.2-12)$$的根全部在复平面上单位圆周之外，或所有根的模$|z|>1$ 1。即p阶自回归过程的平稳条件为
$$\begin{vmatrix}
z
\end{vmatrix}=\sqrt{z_1^2+z_2^2}>1　　(4.2-13)$$
　　z1和z2分别在实部和虚部。
　　当p=1时，(4.2-12)写成
    $$1-{\phi}z=0$$           
　　解方程 $z=\frac{1}{\phi}$平稳条件：$\begin{vmatrix}
z
\end{vmatrix}=\frac{1}{\begin{vmatrix}
{\phi}
\end{vmatrix}}>1$ , 即$|\phi | < 1$
同前面的结论相同。
　　为了研究方便，如果不作特殊说明，总是假定：
　　（1）所有自回归过程都是平稳过程。
当发现时间序列是非平稳的，要清除非平稳性，一般采用差分法。只要对原始数据进行适当阶数的差分处理，便可消除非平稳性。
　　（2）自回归过程中每个元素的期望值都为0即$E(y_t)= 0$。如果实际的时间序列均值$\bar{y}\neq 0$ ，则可对它进行中心化$(y_t-\bar{y})$，中心化后的时间序列必然有零期望值。
###4.2.2 AR模型的自相关函数
　　一阶自相关过程AR(1)的自相关函数，利用(4.1-2)可直接写出
  $$\rho_k=\frac{\gamma_k}{\gamma_0}=\frac{\phi^k\sigma_y^2}{\sigma_y^2}=\phi^k　　(4.2-14)$$
　　AR(p)的自相关函数, 由于 ，将(4.2-2)代入得
$\gamma_k=cov(y_t,y_{t+k})=cov(y_t,y_{t-k})$
$$\begin{align*}
\gamma_k&=\phi_1COV(y_{t-1},y_{t-k})+\phi_2COV(y_{t-2},y_{t-k})+...+\phi_pCOV(y_{t-p},y_{t-k})\\
&=\phi_1\gamma_{k-1}+\phi_2\gamma_{k-2}+...+\phi_p\gamma_{k-p}(k>0)　　　　　　　(4.2-15)
\end{align*} $$
　　当k=0时，
$$\begin{align*}
\gamma_0&=V(y_t)=E(y_t^2)\\
&=E[y_t(\phi_1y_{t-1}+\phi_2y{t-2}+...+\phi_py_{k-p}+\mu_t)]\\
&=\phi_1\gamma_1+\phi_2\gamma_2+...+\phi_p\gamma_p+\sigma_{\mu}^2
\end{align*}$$            
　　对AR(1)便有$\gamma_0=\phi_1\gamma_1+\sigma_{\mu}^2$   ，再由(4.2-14)有$\gamma_1=\phi_1\gamma_0$ ，于是有
$$\gamma_0=\frac{\sigma_{\mu}^2}{1-\phi_1^2}=\sigma_y^2 　　(4.2-16)$$
　　此结果与(4.2-5)相同。
　　用 $\gamma_0=\sigma_y^2$除(4.2-16)式两端，得
 $$\rho=\phi_1\rho_{k-1}+\phi_2\rho_{k-2}+\phi_3\rho_{k-3}+...+\phi_p\rho_{k-p}　　(4.2-17)$$
　　(4.2-17)便是自回归过程AR(p)自相关函数的表达式(也称递推公式)。

acf(y1)$ac[1:9] ＃为方便显示仅列出前9个自相关系数

>[1] 1.0000 0.7587 0.5795 0.4509 0.2966 0.2312 0.1770 0.1792 0.1799
###4.2.3 AR模型的估计与识别
一、 Yule-Walker估计
　　在自相关函数表达式(4.2-17)中，令k=1,2,3,…,p，则得一组方程式，称之为Yule-Walker方程：
$$\left\{\begin{matrix}
\rho_1=\phi_1+\phi_2\rho_1+\phi_3\rho_2+...+\phi_p\rho_{p-1}\\ 
\rho_2=\phi_1\rho_1+\phi_2+\phi_3\rho_1+...+\phi_p\rho_{p-2}\\ 
\rho_3=\phi_1\rho_2+\phi_2\rho_1+\phi_3+...+\phi_p\rho_{p-3}\\ 
...\\ 
\rho_p=\phi_1\rho_{p-1}+\phi_2\rho_{p-2}+\phi_3\rho_{p-3}+...+\phi_p\\
\end{matrix}\right.　(4.2-18)$$
　　其矩阵表达式为：
$$\begin{bmatrix}
1 & \rho_1 &\ rho_2 &...  &\rho_{p-1} \\ 
\rho_1&1  & \rho_2  &...  &\rho_{p-2} \\ 
\rho_1 &\ rho_2  & 1 &...  &\rho_{p-3} \\ 
 ...& ... &...  & ... &... \\ 
 \rho_{p-1}&\ rho_{p-2} &\rho_{p-3}  & ... & 1
\end{bmatrix}\begin{bmatrix}
\phi _1\\ 
\phi _2\\ 
\phi _3\\ 
...\\ 
\phi _p
\end{bmatrix}=\begin{bmatrix}
\rho _1\\ 
\rho _2\\ 
\rho _3\\ 
...\\ 
\rho _p
\end{bmatrix} 　　(4.2-19)$$
　　简记为
$$P_p\Phi_p=\rho_p  　　(4.2-20)$$                                   
　　其中$P_p$为pXp矩阵，$\Phi_p$和$\rho_p$均为pX1阶向量。
　　利用Yule-Walker方程式，就可估计$\Phi_p$值。
　　显然，对AR(1)模型，有 $\phi_1=\rho_1$。
　　而，对AR(2)模型，由(4.2-19)便有
 $$\begin{pmatrix}
1 &\rho_1 \\ 
\rho_1 & 1
\end{pmatrix}\begin{pmatrix}
\phi_1\\ 
\phi_2
\end{pmatrix}=\begin{pmatrix}
\rho_1\\ 
\rho_2
\end{pmatrix}$$                                          
　　解之得
$$\begin{pmatrix}
\phi_1\\ 
\phi_2
\end{pmatrix}={\begin{pmatrix}
1 &\rho_1 \\ 
\rho_1 & 1
\end{pmatrix}}^{-1}\begin{pmatrix}
\rho_1\\ 
\rho_2
\end{pmatrix}$$  
　　于是
$$\left\{\begin{matrix}
\phi _1=\frac{\rho_1-\rho_1\rho_2}{1-\rho_1^2}\\ 
\phi _2=\frac{\rho_2-\rho_1^2}{1-\rho_1^2}
\end{matrix}\right.  　　(4.2-21)$$                                    　　如果将$\rho_1$ 和$\rho_2$ 的估计值$\hat{\rho_1}$和$\hat{\rho_2}$ 代入(4.2-21)，便可求得$\phi_1$ 和$\phi_2$  的估计值$\hat{\phi_1}$  和$\hat{\phi_2}$ 。
　　由(4.2-20)可得：
$$\Phi_p=P_p^{-1}\rho_p　　(4.2-22)$$  
　　$\Phi_p$中最后一个参数$\phi_p$ 称为偏自相关系数，序列${\phi_p }(p = 1，2，3，…)$称为偏自相关函数。(4.2-19)式表示，当自回归模型的阶为p时，则偏自相关函数$\phi_{p+1}$及其后的$\phi$值皆为零。例如，当自回归模型的阶数为2时，则 $\phi_3$及其后的$\phi$ 值皆为零。

ar(y1, method = "yule-walker")

>Call: ar(x = y1, method = "yule-walker")
>Coefficients:
    1  
0.759  
>Order selected 1  sigma^2 estimated as  0.935
　　模型的拟合结果$\hat{\phi_１}=0.759$,跟理论值$\phi_1=0.8$十分接近。
二、最小二乘估计
　　我们可以将(4.2-2)看成因变量为$y_t$，自变量为$y_{t-1}，y_{t-2}，…，y_{t-p}$的线性回归模型，并可用OLS法得出参数估计值。对于自回归模型(4.2-2)
$$y_t=\phi_0+\phi_1y_{t-1}+\phi_2y_{t-2}+..+\phi_py_{t-p}+\mu_t　　(4.2-23)$$
　　跟多元线性回归模型一样，可将其写成矩阵形式
$$Y_p=X_p\Phi_p+U 　　(4.2-24)$$                                 
　　对(4.2-24)应用最小二乘法，得参数估计
 $$\hat{\Phi}_p=(X_p^{'}X_p)^{-1}X_p^{'}Y_p　　(4.2-25)$$                
　　应该指出，估计量 $\hat{\Phi}_p$虽然不是无偏的，却是一致估计量，还是可以接受的。
　　应该指出，此方法所得到的参数估计值$\hat{\phi}_p$精度较低，只能作为首次估计，不宜实际使用。原因是一方面，当样本容量小于50时，样本自相关函数很可能低估，导致$\hat{\phi}_p$的误差扩大，另一方面，样本自相关函数提供的情况，不如时序原始资料本身。比较精确的估计，应该是由普通最小二乘法得到的参数估值$\hat{\phi}_p$ 。

ar(y1, method = "ols")

>Call: ar(x = y1, method = "ols")
>Coefficients:
    1  
0.793  
>Intercept: -0.0726 (0.126) 
>Order selected 1  sigma^2 estimated as  0.783
　　模型的拟合结果$\hat{\phi}_1 = 0.793$跟理论值$\phi_1=0.8$基本一致。
　　从AR模型(4.2-23)的形式看，其相当于$y_t$的滞后各阶回归模型，所以从形式上我们也可以用线性回归模型来求解AR模型，但理论上来说，这样做是有风险的，因为数据滞后会损失数据的。下面我们试用lm函数来求解AR(1)模型。

lm4.1=lm(y1~L_(y1)); lm4.1

>Call: lm(formula = y1 ~ L_(y1))
>Coefficients:
(Intercept)   　　    L_(y1)  
     -0.465      　　　　  0.793  

summary(lm4.1)

>Call: lm(formula = y1 ~ L_(y1))
>Residuals:
   Min     1Q Median     3Q    Max 
-1.877 -0.547 -0.089  0.437  2.865 
>Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -0.4651     0.2107   -2.21    0.032 ＊
>L_(y1)        0.7930     0.0898    8.83  1.5e-11 ＊＊＊
－－－
>Signif. codes:  0 ‘＊＊＊’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>Residual standard error: 0.904 on 47 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.624,     Adjusted R-squared:  0.616 
F-statistic:   78 on 1 and 47 DF,  p-value: 1.51e-11
　　模型的拟合结果跟理论值基本一致，但截距项显著，不符合拟合的模型，且该模型在建模过程中损失了一个数据(1 observation deleted due to missingness)。详细的分析见4.5节分布滞后与自回归模型。
　　由于R语言给出了一个模拟各种时间序列的函数(arima.sim)和拟合时间序列模型的函数(arima)，用它可以模拟AR模型、MA模型、ARMA模型以及ARIMA模型。所以为了统一和方便起见，我们下面的模拟和拟合都采用R自带的函数arima.sim和arima。在以后的估计中，最好统一用arima函数来估计这些参数，而且其估计效果较好。如果进一步还需了解模型的假设检验情况，可用我们改进的arima.test函数。
　　【例4.2-2】模拟研究。模拟AR(2)模型： $y_t = 0.7y_{t-1} - 0.5y_{t-2} + e_t$

set.seed(123)
y2=arima.sim(n=100,list(ar=c(0.7,-0.5))) ＃模拟AR(2)
plot(y2,type='o')

![图片标题](https://leanote.com/api/file/getImage?fileId=583a3767ab644137810150a3)

pacf(y2) ＃偏自相关图

![图片标题](https://leanote.com/api/file/getImage?fileId=583a37b7ab6441366c014bec)

pacf(y2)$ac[1:5] ＃前五个偏自相关系数

>[1]  0.471659 -0.486390  0.009216 -0.043939  0.138181
　　从偏自相关图中可以看出，只有前两个偏自相关系数超出可信区间界限意外，符合我们的理论模型（模拟序列）
```。
arima(y2,order=c(2,0,0))     ＃包含常数项
<div class="md-section-divider"></div>

Call: arima(x = y2, order = c(2, 0, 0))

Coefficients:
ar1 ar2 intercept
0.712 -0.494 -0.007
s.e. 0.087 0.087 0.112
sigma^2 estimated as 0.762: log likelihood = -128.7, aic = 265.4

arima(y2,order=c(2,0,0),include.mean=F)  ＃不包含常数项F即FALSE
<div class="md-section-divider"></div>

Call: arima(x = y2, order = c(2, 0, 0), include.mean = F)

Coefficients:
ar1 ar2
0.712 -0.494
s.e. 0.087 0.087

sigma^2 estimated as 0.762: log likelihood = -128.7, aic = 263.4

　　模型的拟合结果$\hat{\phi_1}=0.712，\hat{\phi_2}=-0.495$跟理论值$\phi_1 = 0.8$，$\phi_2 = -0.5$更接近。
　　注意：
　　（1）arima拟合函数可选择不同的方法进行参数估计，例如arima(prop, order = c(1,0,0), method="ML"))表示用最大似然估计（ML）进行AR(1)模型参数的估计，如参数method="CSS"，估计方法为条件最小二乘法，用条件最小二乘法时，不显示AIC。
　　（2）arima函数里的intercept估计的是均值，而不是截距！虽然intercept是截距的意思，这里如果用mean会更好。(the mean and the intercept are the same only when there is no AR term，均值和截距是相同的，只有在没有AR项的时候)，如果想得到截距，利用公式计算。$\phi_0 = (1-\phi_1).intercept$。
　　（3）arima函数也有一个缺点，就是没有给出系数的显著性检验，我们将其进行了扩展，给出了其检验统计量的函数arima.test（见包EmRa）。

AR2=arima(y2,order=c(2,0,0))
arima.test(AR2)             ＃ library(EmRa)
<div class="md-section-divider"></div>

arima(x = y2, order = c(2, 0, 0))

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ar1 0.7120 0.0872 8.1637 0.000
ar2 -0.4939 0.0875 -5.6470 0.000
intercept -0.0066 0.1121 -0.0589 0.953

sigma^2 estimated as 0.7619 :log likelihood = -128.7 aic = 265.4

AR2=arima(y2,order=c(2,0,0),include.mean=F)
arima.test(AR2)
<div class="md-section-divider"></div>

arima(x = y2, order = c(2, 0, 0), include.mean = F)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ar1 0.7120 0.0872 8.164 0
ar2 -0.4941 0.0874 -5.652 0

sigma^2 estimated as 0.7619 :log likelihood = -128.7 aic = 263.4

　　经检验，模型系数显著，而常数项不显著，跟模拟的理论模型参数基本一致。
三、自回归阶数p的识别
　　对AR(p)模型，关键是模型的识别即如何确定阶数p，一旦p值确定下来就转化为自回归阶数p已知的情况，问题就解决了。这里我们采用直观的偏自相关系数法。
　　自相关函数acf(k)给出了$X_t$与$X_{t-k}$的总体相关性，但总体相关性可能掩盖了变量间完全不同的隐含关系。
　　例如，在AR(1)随机过程中，$X_t$与$X_{t-2}$间有相关性可能主要是由于它们各自与$X_{t-1}$间的相关性带来的:

$$\rho_2=\phi^2=\rho_1^2=E(X_t{X_{t-1}})E(X_{t-1}{X_{t-2}})$$
即自相关函数中包含了这种所有的“间接”相关。
　　与之相反，$X_t$与$X_{t-k}$间的偏自相关函数(partial autocorrelation，简记为pacf)则是消除了中间变量$X_{t-1}，…，X_{t-k+1}$ 带来的间接相关后的直接相关性，它是在已知序列值$X_{t-1}，…，X_{t-k+1}$的条件下，$X_t$与$X_{t-k}$间关系的度量。
　　偏自相关的含义类似于偏相关，偏相关系数是指在排除了k-1个中间量滞后变量$y_{t-1},y_{t-2}……y_{t-k+1}$的影响后，$y_t$和$y_{t-k}$的自相关系数。在实际应用中偏自相关系数函数可由样本的自相关函数值通过递推的方法给出，其关系为：
$$\begin{matrix} \hat{\phi }_1=\hat{\rho }_1\\ \hat{\phi }_{k+1,j}=\hat{\phi }_{k,j}+\hat{\phi }_{k+1,k+1}\hat{\phi }_{k,k-j+1}　　j=1,2,...k\\ \hat{\phi }_{k+1,k+1}=(\hat{\rho }_{k+1}-\sum_{j=1}^{k}\hat{\rho }_{k+1-j}\hat{\phi }_{k,j})(1-\sum_{j=1}^{k}\hat{\rho }_j\hat{\phi }_{k,j})^{-1} 　　(4.2-26) \end{matrix}$$
　　可以证明，平稳的零均值序列${y_t}$为p阶自回归的充分必要条件是${y_t}$的偏自相关函数为p阶截尾(当$k\leq{p}$时，$\psi_k\neq0$，而当$k>p$时，$\psi_k=0$)，因此若序列${y_t}$的偏自相关函数在p步以后截尾，则可判断该序列为AR(p)。而由于任何可逆的MA(q)序列都可以转化成一个无限阶的AR序列，所以MA(q)序列的偏自相关函数也呈拖尾性。
　　在AR(1)中，从$X_t$中去掉$X_{t-1}$的影响，则只剩下随机扰动项$\varepsilon_t$，显然它与$X_{t-2}$无关，因此我们说$X_t$与$X_{t-2}$的偏自相关系数为零，记为
$$\phi_2 = Corr(\mu_t,X_{t-2}) = 0$$
　　同样地，在AR(p)过程中，对所有的k>p，$X_t$与$X_{t-k}$间的偏自相关系数为零。
　　AR(p)的一个主要特征是:k>p时，$\phi_k = Corr(X_t,X_{t-k}) = 0$ ，即$\phi_k$在p以后是截尾的。
　　一随机时间序列的识别原则：
　　若$X_t$的偏自相关函数在p以后截尾，即k>p时，$\phi_k=0$，而它的自相关函数$\phi_k$是拖尾的，则此序列是自回归AR(p)序列。
　　这种方法是在自回归阶数k逐步增加的过程中，通过对偏自相关系数$\phi_k$的显著性检验来确定适当阶数p的方法。偏自相关系数中的第 个系数$\phi_k$我们用$\phi_kk$表示。
　　为了对$\phi_kk$进行检验，必须知道OLS估计量的抽样分布。可以证明，对于大样本来讲，如果自回归过程(AR)的阶数为p，那么，在k > p时，偏自相关系数估计量近似服从期望值为0，方差为1/n的正态分布，这里的n为样本容量。
　　要判断在0.05显著性水平下$\phi_kk$是否为0，只要考察其估计值是否落在下面区间内：
$$[-\frac{1.96}{\sqrt{n}},\frac{1.96}{\sqrt{n}}]$$
　　如果估计值落在这个区间内，则$\phi_kk$不显著，即确认$\phi_kk= 0$，如果估计值落在此区间之外，则$\phi_kk$显著，即确认$\phi_kk\neq0$。

pacf(y1)$ac[1:5]         ＃模拟序列y1的偏自相关图及前5个偏自相关系数
<div class="md-section-divider"></div>

[1] 0.758746 0.008968 0.019635 -0.126983 0.096058

　　从中可以看出，y1序列只存在一阶自相关性，二阶及以上偏自相关系数很小，且都落入可信区间之内，符合AR(1)模型的拟合结果。

pacf(y2)$ac[1:5]  
＃模拟序列y2的偏自相关图及前5个偏自相关系数
<div class="md-section-divider"></div>

[1] 0.471659 -0.486390 0.009216 -0.043939 0.138181

　　从中可以看出，$y_2$序列的确只存在二阶自相关性，三阶及以上偏自相关系数很小，且都落入可信区间之内，符合拟合的AR(2)模型。
　　【例4.2-2】实证分析：股票收益率的研究。
下面对上海股票市场股票指数收益率的变动进行分析。本例选取2003年1月2日至2010年4月14日的上海证券交易所每日收盘价格指数作为样本数据，总共1776个数值，数据来源于雅虎财经，并指数的对数收益率$R_t$作为第t日的股票对数收益率，即这里的$R_t=ln(P_t/P_{t-1})*100$，$P_t$表示股票价格指数。

＃从SCIdata.csv文件中读取上证指数; Pt=read.csv("SCIdata.csv")$Pt
Pt=read.csv("http://emlab.jnu.edu.cn/EmRa/SCIdata.csv")$Pt    
Rt=log(Pt[-1]/L_(Pt))*100   ＃滞后函数L_会损失1组数据，需剔除Rt中第1组数据
plot(Pt,type='l')           ＃上证指数图 
plot(Rt,type='l');abline(h=0)  ＃上证指数对数收益率图
<div class="md-section-divider"></div>

　　从上证指数$P_t$图我们看到上证指数是不平稳序列，但从其对数收益率$R_t$的图上可以看出，上证指数的对数收益率$R_t$是一个平稳序列，下面我们来建立$R_t$的AR模型。

＃应用偏自相关函数确定模型的阶数
pacf(Rt)   ＃初步定仅有4阶偏自相关
<div class="md-section-divider"></div>

Box.test(Rt,lag=1)
<div class="md-section-divider"></div>

   Box-Pierce test

data: Rt
X-squared = 2.979, df = 1, p-value = 0.08436

Box.test(Rt,lag=2)
<div class="md-section-divider"></div>

   Box-Pierce test

data: Rt
X-squared = 4.544, df = 2, p-value = 0.1031

Box.test(Rt,lag=3)
<div class="md-section-divider"></div>

   Box-Pierce test

data: Rt
X-squared = 7.739, df = 3, p-value = 0.05172

Box.test(Rt,lag=4)
        Box-Pierce test
data:  Rt
X-squared = 14.87, df = 4, p-value = 0.004973
<div class="md-section-divider"></div>

　　Box-Pierce检验进一步确认该序列仅存在4阶偏自相关关系。

AR4=arima(Rt,order=c(4,0,0),include.mean=F)  ＃ 对数收益率均值为0
arima.test(AR4)
<div class="md-section-divider"></div>

arima(x = Rt, order = c(4, 0, 0), include.mean = F)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ar1 0.04272 0.02375 1.798 0.07211
ar2 -0.02982 0.02375 -1.256 0.20929
ar3 0.04416 0.02374 1.860 0.06288
ar4 0.06063 0.02375 2.553 0.01069

sigma^2 estimated as 3.93 :log likelihood = -3712 aic = 7434

　　于是有AR(4)模型：
$(1-0.04272L+0.02982L^2-0.04416 L^3-0.06063L^4) R_t = u_t$
即$R_t = 0.04272 R_{t-1} - 0.02982 R_{t-2}+0.04416 R_{t-3}+0.06063 R_{t-4} + u_t$
$\sigma_u^2 = 3.93$
　　从假设检验结果可以看到前三阶自回归系数都不显著，所以我们可建立仅包含4阶自回归的AR模型，其中fixed=c(0,0,0,NA)表示只拟合第4阶模型，transform.pars=FALSE使输出结果简单些。

AR4=arima(Rt,order=c(4,0,0),include.mean=F,fixed=c(0,0,0,NA) ,transform.pars=F)
arima.test(AR4)
<div class="md-section-divider"></div>

arima(x = Rt, order = c(4, 0, 0), include.mean = F, transform.pars = FALSE,
fixed = c(0, 0, 0, NA))

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ar1 0.00000 NA NA NA
ar2 0.00000 NA NA NA
ar3 0.00000 NA NA NA
ar4 0.06528 0.02375 2.749 0.005971

sigma^2 estimated as 3.947 :log likelihood = -3716 aic = 7436

　　于是得最终的AR(4)模型为：

$$(1-0.06528L^4) R_t = u_t 或 R_t = 0.06528 R_t-4 + u_t$$
　　　　$\sigma_u^2 = 3.947$
　　但从AIC的值看，两者差别不大，但后者模型更为简单。

4.3 时间序列移动平均MA模型

4.3.1 MA模型的基本形式

　　移动平均模型将序列${y_t}$表示为白噪声的线性加权。q阶移动平均模型表示为：

$$y_t = u_t+\theta_1{u_{t-1}}+\theta_2{u_{t-2}}+……\theta_p{u_{t-q}} 　　(4.3-1)$$
　　利用滞后算子则模型可简单表示为：
$$y_t = \Theta_q{(L)} u_t 　　(4.3-２)$$
　　或
$$u_t = \Theta_q^{-1}(L)y_t　　(4.3-３)$$
　　其中
$$\Theta_q(L) = 1+\theta_1L+\theta_2{L^2}+……\theta_q{L^q} 　　(4.3-4)$$
　　从上式可见，为使MA(q)过程可以转换成一个自回归过程，需要$\Theta_q(L)$收敛。而$\Theta_q(L)$收敛的充分必要条件是$\Theta_q{(L)}$的特征方程
$$\Theta_q(z) = 1+\theta_1z+\theta_2{z^2}+……\theta_q{z^q} 　　(4.3-4)$$
的所有的根的绝对值都大于1，即|z|>1，这个条件称为移动平均序列的可逆性(可转换性) 条件，可逆性条件是MA(q)序列必须满足的。
上式为 阶移动平均模型，记为MA(q)。c是常数项，实参数$\theta_1,\theta_2K, \theta_q$为移动平均系数，是模型的待估参数。

4.3.2 MA模型的阶数确定

　　对于给定的样本$y_1，y_2，…，y_n$，怎样为生成移动平均过程确定合适的阶数？为了回答这个问题，我们首先来研究反映移动平均过程特征的自相关函数。
一、自相关函数法
　　1. MA(1)模型
　　为了讨论方便，我们先研究MA(1)过程

$$y_t=u_t-\theta{u_{t-1}}　　 (4.3-5)$$
　　因为$u_t$和$u_{t-1}$ 是白噪声，所以，它的期望值为
$$E(y_t)=E(u_t)-\theta_1E(u_{t-1})=0$$
　　方差为
$$\begin{align*} var(y_t)&=E(u_t-\theta_1{u_{t-1}})^2\\ &=E(u_t^2)+\theta_1^2E(u_{t-1}^2)-2\theta_1E(u_tu_{t-1})\\ &=\sigma_u^2+\theta_1^2\sigma_u^2 \end{align*}$$
　　协方差
$$\begin{align*} cov(y_t,y_{t+k}&=E[(u_t-\theta_1{u_{t-1}})(u_{t+k}-\theta_1{u_{t+k-1}}\\ &=E(u_tu_{t+k})-\theta_1E(u_{t-1}u_{t+k})-\theta_1E(u_tu_{t+k-1})+\theta_1^2E(u_{t-1}u_{t+k-1})\\ &=\left\{\begin{matrix} (1+\theta_1^2)\sigma_u^2 &k=0 \\ -\theta_1\sigma_u^2&k=1 \\ 0&k=2,3,... \end{matrix}\right.　　(4.3-6) \end{align*}$$
　　以上讨论表明MA(1)过程是平稳的。由于$cov(y_t,y_{t+k}）$ 　　不依赖时间t，而只依赖ｋ，所以可以用$r_k$表示。于是，自相关函数为
$$\rho_k=\frac{r_k}{r_0}=\left\{\begin{matrix} \frac{-\theta_1}{1+\theta_1^2}&k=1\\ 0&k=2,3,... \end{matrix}\right.　　(4.3-7)$$
　　上式表明MA(1)只有1期记忆，即当 k＞1时$\rho_k=0$ 。
　　（2）MA(q)模型：
$$y_t=u_t-\theta_1u_{t-1} -\theta_2u_{t-2}-...--\theta_ku_{t-k}-...--\theta_qu_{t-q}(k\leq q) 　　(4.3-8)$$
　　与MA(1)的推导过程类似，可得结果：
$$r_0=\sigma_u^2(1+\theta_1^2+...+\theta_q^2)$$
$$r_k=\sigma_u^2(-\theta_k+\theta_1\theta_{k+1}+\theta_2\theta_{k+2}+...+\theta_{q-k}\theta_q)(k=1,2,...,q)$$
　　于是
$$\rho_k=\frac{r_k}{r_0}=\left\{\begin{matrix} 1&k=0\\ \frac{-\theta_k+\theta_1\theta_{k+1}+\theta_2\theta_{k+2}+...+\theta_{q-k}\theta_q}{1+\theta_1^2+...+\theta_q^2}&k=1,2,..,q\\ 0&k>q \end{matrix}\right.　　(4.3-9)$$

　　对平稳过程，方差$var(y_t=r_0$ 必须有限，因此要求$(1+\theta_1^2+...+\theta_q^2)<\infty$ ，对无穷阶移动平均过程要求$(1+\theta_1^2+\theta_2^2+...)<\infty$ 。这意味着$\theta_q$的绝对值必须随q的增大而减少。由(4.3-9)试算可以看出，MA(q)的$|\rho_k|$ 也将随 的增大而减少，与自回归过程不同的是当 k＞q时，$|\rho_k|=0$ 。这表明MA(q)只有q期记忆，即当 k＞q时，$|\rho_k|=0$。
二、阶数的确定
　　MA(q)只有q期记忆这一重要性质，可以帮助我们对模型进行识别。
　　1.先计算自相关函数$\rho_k$ 的估计值，

$$\hat{\rho_k}=\frac{\sum_{t=k+1}^{n}y_ty_{t-k}}{\sum_{t=1}^{n}y_t^2}　　(4.3-10)$$
　　这里根据我们的约定，$y_t$已经中心化了，即$\bar{y_t}=0$。利用上式计算出$\hat{\rho_1} ，\hat{\rho_2} ，……\hat{\rho_k}$。
　　2.以$\hat{\rho_k}(k=0,1,2,...)$为纵坐标，以k为横坐标作图，这样得出的图形称为样本自相关函数图。
　　3.将样本自相关函数图形与总体自相关函数图形比较，由于总体自相关函数图形与阶数p有关，所以通过比较便可确定p值进行显著检验。
可以证明，当样本容量很大时，$\hat{\rho_k}$近似服从期望值为0，方差为1/n的正态分布。
　　于是可以有与自回归过程类似的检验方法：
　　（1）构造95%的置信区间
$$[-\frac{1.96}{\sqrt{n}}, \frac{1.96}{\sqrt{n}}]\approx [-\frac{2}{\sqrt{n}},\frac{2}{\sqrt{n}}]$$
　　（2）计算样本的各阶自相关系数$\hat{\rho_k}$ 。
　　（3）考察$\hat{\rho_k}$ 是否落在这一区间之内。如果$\hat{\rho_k}$ 的数值落在此区间之外，表明 $\hat{\rho_k}$显著(即 $\hat{\rho_k}\neq0$)，否则不显著(即$\hat{\rho_k}=0$ )。若对 k≤q时， $\hat{\rho_k}$皆显著，对 k＞q时， $\hat{\rho_k}$皆不显著，则在0.05显著水平下，产生样本的移动平均的阶数确定为q。

4.3.3 MA模型的参数估计

　　设有移动平均过程MA(q)

$$y_t=u_t-\theta_1u_{t-1}-\theta_2u_{t-2}-...-\theta_qu_{t-q} 　　(4.3-11)$$
　　估计参数$\theta$ 的直观想法是利用$\rho_k$ 与$\theta$ 之间的关系
$$\rho_k=\frac{-\theta_k+\theta_1\theta_{k+1}+...+\theta_{q-k}\theta_q}{1+\theta_1^2+\theta_2^2+...+\theta_q^2}k=1,2,..,q 　　(4.3-12)$$
　　将此式中的 $\rho_k$用相应的估计值 $\hat{\rho_k}$代替，便得到关于 $\theta$ 的q个非线性代数方程，解这个非线性方程便可得到q个$\theta$ 的估计值。这样得出的参数估计量是一致的。
非线性方程组求解是比较困难，一般可用迭代法近似求解。由于实际问题中遇到的移动平均过程多数都是MA(1)或MA(2)等低阶过程，所以迭代过程也不算太繁。
　　这种估计方法的优点是简单易行，缺点是它只注意如何去拟合q个自相关函数，而没有考虑如何选择参数有利于提高预测精度。因此，往往把它作为首期估计值。为此可以考虑误差平方和：
$$S(\theta)=\sum_{t=1}^nu_t^2 　　(4.3-13)$$
　　选择适当的参数 $\theta$使$S(\theta)$ 达到极小，可得到参数$\theta$ 的一致估计量。
　　【例4.3-1】模拟研究。模型MR(1)估计的结果为$y_t = u_t - 0.8u_{t-1}$

set.seed(123456)
y3=arima.sim(n=100,list(ma=0.8))
plot(y3,type='o')
<div class="md-section-divider"></div>

acf(y3)$ac[1:5]
<div class="md-section-divider"></div>

[1] 1.00000 0.56589 0.05155 -0.05005 -0.04508

MA1=arima(y3,order=c(0,0,1),include.mean=F)
arima.test(MA1)
<div class="md-section-divider"></div>

arima(x = y3, order = c(0, 0, 1), include.mean = FALSE)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ma1 0.8984 0.0428 20.97 0

sigma^2 estimated as 0.9393 :log likelihood = -139.6 aic = 283.2

　　对$\rho_k$ 进行显著性检验，可以看出，只有$\hat{\rho_1}=0.56589$ 落在置信区间之外，其余$\hat{\rho_k}$ 皆落在区间之内，表明，该样本的移动平均过程阶数为1，即选定MA(1)过程：$y_t = u_t - \theta_1u_{t-1}$
　　　由模拟样本生成的MA(1)过程为$：y_t = u_t - 0.8984u_{t-1}$。可以看到，模拟结果与理论模型也比较接近：$：y_t = u_t - 0.8u_{t-1}$
　　【例4.3-2】股票收益率$R_t$的MA模型建立。
　　从acf函数图我们可以看到，股票收益率Rt的自相关图在4阶处截尾，可考虑拟合MA(4)模型进行分析。

acf(Rt)
<div class="md-section-divider"></div>

MA4=arima(Rt,order=c(0,0,4),include.mean=F)
arima.test(MA4)
<div class="md-section-divider"></div>

arima(x = Rt, order = c(0, 0, 4), include.mean = F)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ma1 0.04496 0.02371 1.8966 0.057882
ma2 -0.02133 0.02386 -0.8939 0.371369
ma3 0.04115 0.02490 1.6526 0.098414
ma4 0.06396 0.02401 2.6635 0.007732

sigma^2 estimated as 3.931 :log likelihood = -3712 aic = 7435

　　经检验，只有滞后4阶的模型是显著的，于是我们只拟合包含4阶的MA(4)模型：

MA4=arima(Rt,order=c(0,0,4),include.mean=F,fixed=c(0,0,0,NA))
arima.test(MA4)
<div class="md-section-divider"></div>

arima(x = Rt, order = c(0, 0, 4), include.mean = F, fixed = c(0,
0, 0, NA))

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ma1 0.00000 NA NA NA
ma2 0.00000 NA NA NA
ma3 0.00000 NA NA NA
ma4 0.06609 0.02391 2.764 0.005703

sigma^2 estimated as 3.947 :log likelihood = -3716 aic = 7436

　　于是得如下的MA(4)模型：
　　$y_t = (1-0.06609 L^4)u_t 或 y_t = u_t -0.06609u_{t-1}$
　　$\sigma_u^2 = 3.947$
　　但从AIC的值看，MA(4)模型跟AR(4)模型差别不大。

4.4 自回归移动平均ARMA模型

4.4.1 ARMA模型的概念

　　ARMA模型是一类常用的随机时间序列模型，由美国统计学家博克斯(G.E.P.Box)和詹金斯(G.M.Jenkins)创立，亦称B-J方法。它是一种精度较高的时间序列短期分析方法，其基本思想是：某些时间序列是依赖于时间t的一族随机变量，构成该时间序列的单个序列值虽然具有不确定性，但整个序列的变化却有一定的规律性，可以用相应的数学模型近似描述，通过对该数学模型的分析研究，能够从本质上认识时间序列的结构与特征，并得到最有效的预测结果。
　　如果平稳随机过程既具有自回归过程的特性又具有移动平均过程的特性，则不宜单独使用AR(p)或MA(q)模型，而需要两种模型混合使用。由于这种模型包含了自回归和移动平均两种成分，所以它的阶是二维的，由p和q两个数构成，其中p代表自回归成分的阶数，q代表移动平均成分的阶数，记作ARMA(p，q)，称作自回归移动平均混合模型或称为自回归移动平均模型。
　　自回归移动平均模型ARMA(p,q)的一般表达式为

$$y_t=\phi_1y_{t-1}+\phi_2y_{t-2} +...+\phi_py_{t-p}+u_t-\theta_1u_{t-1}--\theta_2u_{t-2}-...-\theta_qu_{t-q}　　(4.4-1)$$
　　该模型中p为自回归部分的阶数，q为移动平均部分的阶数，因此记为ARMA(p,q)。利用滞后算子多项式$\phi_p(L)$和$\Phi_q(L)$，ARMA(p,q)可改写为
$$\phi_p(L)y_t=\Phi_q(L)u_t 　　(4.4-2)$$
　　对ARMA(p,q)模型，其平稳条件与AR(p)的平稳条件相同，其可逆条件与MA(q)的可逆条件相同。
　　最简单的自回归移动平均模型是ARMA(1,1)，其具体形式为：
$$y_t=\phi_1y_{t-1}+u_t-\theta_1u_{t-1} 　　(4.4-3)$$
　　显然，$ARMA(0,q)=MA(q)，ARMA(p,0)=AR(p)$，因此，MA(q)和AR(p)分别看作ARMA(p，q)，当p=0和q=0时的特例。
　　ARMA(p,q)模型的优点是能以较少的参数描写单用AR(p)或MA(q)过程不能经济地描写的数据生成过程。在实际应用中，用ARMA(p,q)拟合实际数据时所需阶数较低，p和q的数值很少超过2。因此，ARMA模型在预测中具有很大的实用价值。

4.4.2 ARMA模型的相关分析

一、自相关分析
　　首先我们考察ARMA(p,q)中最简单的模型ARMA(1,1)：

$$y_t=\phi_1y_{t-1}+u_t-\theta_1u_{t-1} 　　(4.4-4)$$
$$\begin{align*} r_0&=E(y_t^2)=E(\phi_1y_{t-1}+u_t-\theta_1u_{t-1})^2\\ &=\phi_1^2r_0+\sigma_u^2+\theta_1^2\sigma_u^2-2\phi_1\theta_1\sigma_u^2 \end{align*}$$
$$r_0(1-\phi_1^2)=\sigma_u^2(1+\theta_1^2-2\phi_1\theta)$$
　　　　　　$r_0=\frac{1+\theta_1^2-2\phi_1\theta_1}{1-\phi_1^2}\sigma_u^2$
　 $$\begin{align*} 　r_1&=E(y_{t-1}y_t)=E[y_{t-1}(\phi_1y_{t-1}+u_t-\theta_1u_{t-1})]\\ &=\phi_1r_0-\theta_1\sigma_u^2\\ &=\frac{(\phi_1-\theta_1)(1-\phi_1\theta_1)}{1-\phi_1^2}\sigma_u^2\end{align*}$$
　　　$\begin{align*} 　　　r_2&=E[y_{t-2}y_t]\\ 　　　&=E[y_{t-2}(\phi_1y_{t-1}+u_t-\phi_1u_{t-1})]\\ 　　　&=\phi_1r_1 　　　\end{align*}$
　　同样有
　　　　　　　$r_3=\phi_1r_2$
　　　　　　　$......$
　　　　　　　$r_k=\phi_1r_{k-1}$
　　于是有自相关函数
$$\rho_1=\frac{r_1}{r_0}=\frac{(\phi_1-\theta_1)(1-\phi_1\theta_1)}{1+\theta_1^2-2\phi_1\theta_1}$$ $$\rho_k=\phi_1\rho_{k-1}(k\geq 2) 　　(4.4-5)$$
二、偏自相关分析
　　主要考虑一下ARMA(1,1)的偏自相关系数。
　　将AR(1)式推迟一期乘以$\theta_1$，推迟二期乘以$\theta_1^2 ……$得到
$$y_t=\phi_1y_{t-1}+u_t-\phi_1u_{t-1}$$
$$\theta_1y_{t-1}=\theta_1\phi_1y_{t-2}+\theta_1u_{t-1}-\theta_1^2u_{t-2}$$
$$\theta_1^2y_{t-2}=\theta_1^2\phi_1y_{t-3}+\theta_1^2u_{t-2}-\theta_1^3u_{t-3}$$
　　　　　　　　$...$

　　将以上各式两边对应相加：

$$y_t+theta_1y_{t-1}+\theta_1^2y_{t-2}+...=\phi_1y_{t-1}+\theta_1\phi_1y_{t-2}+\theta_1^2\phi_1y_{t-3}+...+u_t$$
　　整理得无限自回归模型：
$$y_t=(\phi_1-\theta_1)(y_{t-1}+\theta_1y_{t-2}+\theta_1^2y_{t-3}+...)+u_t$$
　　该模型与ARMA(1,1)等价。
　　记$\psi_1=(\phi_1-\theta_1),\psi_2=(\phi_1-\theta_1)\theta_1,\psi_3=(\phi_1-\theta_1)\theta_1^2,...$ 显然，$\psi$ 就是模型ARMA(1,1)的偏相关系数。由于$|\theta_1|<1,|\phi_1|<1$， ，知$\psi_k$的绝对值随k的增加而减少，其符号由乘积$(\phi_1-\theta_1)\theta_1^k(k=1,2,...)$ 来决定。关于零阶偏相关系数定义为$\psi_0=1$ 。
　　根据以上的讨论，我们可以把自相关函数和偏自相关函数的变化规律，画出图形。下列图给出了典型ARMA(1,1)过程的自相关函数acf和偏自相关函数pacf的图形。
　　【例4.4-1】模拟研究。模型ARMR(1,1)估计的结果为 　　$y_t = 0.8y_{t-1}+ u_t - 0.6u_{t-1}$

set.seed(12345)
y4=arima.sim(model=list(ar=0.8,ma=0.6),n=100)
plot(y4,type='o')
<div class="md-section-divider"></div>

acf(y4)
<div class="md-section-divider"></div>

pacf(y4)
<div class="md-section-divider"></div>

4.4.3 ARMA模型的统计推断

一、ARMA模型的参数估计
　　ARMA(p，q)的估计，需要用非线性估计法，要比AR(p)或MA(q)的估计困难得多。它的计算过程颇为发杂，但使用计算机则非常方便，在实际工作中，一般都是使用通用软件(如TSP)利用计算机进行估计。
　　我们仅对ARMA(p，q)的一种估计方法作简略介绍。
　　我们的问题是找出适当p,q使$\phi、\theta$ 误差

$$S(\phi、\theta)=\sum_{t=1}^{n}u_t^2=\sum_{t=1}^{n}(u_t|\phi,\theta,\omega )^2$$
达到最小。其中$\omega$是n个观测值的向量， 是q个系数$\theta_1,\theta_2,...,\theta_q$的向量， $\phi$是p个系数$\phi_1,\phi_2,...,\phi_p$ 的向量。 $u_t|\phi,\theta,\omega$表示$u_t$值的大小由条件$\phi,\theta,\omega$及 而定。
　　现在以向量$\beta$代表向量$\phi$和$\theta$，即$\beta=(\phi_1,\phi_2,...,\phi_p,\theta_1,\theta_2,...,\theta_q)$。则$S(\phi、\theta)$可写成：
$$S(\beta)=\sum_{t=1}^{n}(u_t|\omega,\beta) ^2=\sum_{t=1}^{n}u_t^2$$
　　我们将$u_t$ 在 $\beta$处泰劳(Taylor)展开：
$$u_t= (u_t|\omega,\beta_0)+\sum_{i=1}^{p+q}(\beta_i-\beta_{i,0})\frac{\partial{u_t}}{\partial{\beta_i}}|_{\beta=\beta_0}+ \frac{1}{2}\sum_{i=1}^{p+q}\sum_{j=1}^{p+q}(\beta_i-\beta_{i,0})(\beta_j-\beta_{j,0}\frac{\partial^2u_t}{\partial\beta_i\partial\beta_j}|_{\beta=\beta_0}+...$$
　　记$x_{i,t}=\frac{-\partial{u_t}}{\partial{\beta_i}}|_{\beta=\beta_0}$ ，上式右边取前二项作为近似值，于是有：
$$u_t=(u_t|\omega,\beta_0)-\sum_{i=1}^{p+q}(\beta_i-\beta_{i,0})x_{i,t}$$
　　或
$$(u_t|\omega,\beta_0)+\sum_{i=1}^{p+q}\beta_{i,0}x_{i,t}=\sum_{i=1}^{p+q}\beta_ix_{i,t}+u_t$$
　　将上式左端视为因变量，右端视为自变量和误差项。记
$$Y=\begin{bmatrix} [u_1|\omega,\beta_0]+\sum_{i=1}^{p+q}\beta_{i,0}x_{i,1}\\ [u_2|\omega,\beta_0]+\sum_{i=1}^{p+q}\beta_{i,0}x_{i,2}\\ ...\\ [u_n|\omega,\beta_0]+\sum_{i=1}^{p+q}\beta_{i,0}x_{i,n} \end{bmatrix}$$
$$X=\begin{bmatrix} x_{1,1}&x_{2,1}&...&x_{p+q,1}\\ x_{1,2}&x_{2,2}&...&x_{p+q,2}\\ \vdots&\vdots&&\vdots\\ x_{1,n}&x_{2,n}&...&x_{p+q,n} \end{bmatrix} ,\beta=\begin{bmatrix} \beta_1\\ \beta_2\\ \vdots\\ \beta_{p+q} \end{bmatrix},U=\begin{bmatrix} u_1\\ u_2\\ \vdots\\ u_n\\ \end{bmatrix}$$
　　其中Y 为nX1 向量，X为nX(p+q) 矩阵，$\beta$ 为 (p+q)X1向量，U为 mX1向量。则有
$$Y=X \beta+U 　　(4.4-8)$$
　　对模型应用OLS法，便可得到 的最小二乘估计量：
$$\hat{ \beta}=(X^{'}X)^{-1}X^{'}Y　　(4.4-9)$$
　　从原始向量估计值$\hat{ \beta_0}$ 开始，如此反复进行计算可依次得到$\hat{ \beta_1},...,\hat{ \beta_K}$直至差$|\beta_K-\hat{ \beta_K}|$ 小于预先指定的精度为止。

a11.y4=arima(y4,order=c(1,0,1),include.mean=F)
arima.test(a11.y4)
<div class="md-section-divider"></div>

arima(x = y4, order = c(1, 0, 1), include.mean = F)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ar1 0.8296 0.05768 14.382 6.675e-47
ma1 0.6249 0.11073 5.644 1.661e-08

sigma^2 estimated as 1.384 :log likelihood = -159.4 aic = 324.8

a21.y4=arima(y4,order=c(2,0,1),include.mean=F)
arima.test(a21.y4)
<div class="md-section-divider"></div>

arima(x = y4, order = c(2, 0, 1), include.mean = F)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
ar1 0.5479 0.1676 3.269 1.080e-03
ar2 0.2767 0.1643 1.684 9.226e-02
ma1 0.8762 0.1203 7.282 3.281e-13

sigma^2 estimated as 1.36 :log likelihood = -158.7 aic = 325.5

a22.y4=arima(y4,order=c(2,0,2),include.mean=F)
arima.test(a22.y4)
<div class="md-section-divider"></div>

arima(x = y4, order = c(2, 0, 2), include.mean = F)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ar1 0.04147 0.1397 0.2969 7.666e-01
ar2 0.67794 0.1049 6.4602 1.046e-10
ma1 1.44207 0.1739 8.2938 1.097e-16
ma2 0.44207 0.1674 2.6412 8.260e-03

sigma^2 estimated as 1.294 :log likelihood = -157.3 aic = 324.6

　　从上面的建模过程可以看到，ARMA(p,q)模型的阶数p和q的识别是拟合ARMA模型首先遇到的问题，即p和q到底取多大合适？

二、模型的识别
　　模型的识别是根据时间序列的样本自相关函数、偏自相关函数的特点来选择模型的类型，并初步判定时间序列所适合的阶数。
　　（1）AR(p)序列的识别：
如果序列${y_t}$的偏自相关函数在p步后截尾，则可判定该序列为AR(p)序列。
　　（2）MA(q)序列的识别：
如果序列${y_t}$的自相关函数q步后截尾，则可判定该序列为MA(q)序列。
　　（3）ARMA序列的识别：
如果序列${y_t}$的自相关函数和偏自相关函数皆无截断点，也即它们均是拖尾的，则可判定该序列为ARMA序列，如下表所示。
　　　　　　　　　　表4.4-1 ARMA(p,q)模型的acf和pacf理论模式

模型	acf	pacf
AR(p)	衰减趋于0或振荡	p阶截尾
MA(q)	q阶截尾	衰减趋于0或振荡
ARMA(p,q)	q阶后衰减趋于0或振荡	p阶后衰减趋于0或振荡

　　自相关函数和偏自相关函数只能初步断定序列${y_t}$为ARMA模型（因为自相关函数和偏自相关函数都是拖尾的）但是不能确定其阶数，这时需要采用一定的定阶准则，常用的定阶准则是AIC准则，也即赤池(Akaike)信息准则。
　　对于序列${y_t}$，取$\hat{\sigma}_e ^2$为模型拟合残差的方差：

$$AIC(p,q)＝ln\hat{\sigma}_e ^2 +2(p＋q)/T　　(4.4-6)$$ 　
　　对于由低阶到高阶不同的p、q取值分别建立模型并进行参数估计，比较各模型的AIC值，使其达到最小的$p_0,q_0$，这时的($p_0,q_0$)为最佳模型阶数，有
$$AIC(p_0,q_0)＝\min_{0\leq{n}\leq{N_o}} AIC(p,q)　　(4.4-7)$$
　　除了赤池信息准则(AIC)，还有Schwarz信息准则（SIC），即选择最佳的(p,q)使SIC达到最小。
　　ARMA模型的阶数确定是一个困难的事，目前还没有最好的方法，目前有两种R语言方法可供大家使用。
　　（1）扩展的自相关函数（Extended Autocorrelation Function，EACF）法。
　　eacf函数包含着TSA包中，需安装TSA包。
　　（2）自动arima模型确定法（Auto ARIMA）
　　 auto.arima函数包含着forecast包中，需安装forecast包。
从使用效果看，eacf函数不如auto.arima函数，因为eacf函数是一种直观方法，还需人为去判断，而auto.arima函数可自动给出ARMA(p.q)的p和q值。所以本书主张使用auto.arima函数来确定ARMA(p.q)的p和q。

library(forecast)           ＃需先安装包forecast
aa.y4=auto.arima(y4); 
aa.y4
<div class="md-section-divider"></div>

Series: y4
ARIMA(1,0,1) with non-zero mean

Coefficients:
ar1 ma1 intercept
0.741 0.653 1.881
s.e. 0.074 0.111 0.714

sigma^2 estimated as 1.33: log likelihood=-157.2
AIC=322.4 AICc=322.8 BIC=332.8

　　auto.arima函数自动给出了模型的参数p=1，q=1，符合我们模拟的模型。
　　【例4.4-2】股票收益率数据的ARMA模型建立。

library(forecast)            
aa.Rt=auto.arima(Rt); aa.Rt
<div class="md-section-divider"></div>

Series: Rt
ARIMA(4,0,4) with zero mean

Coefficients:
ar1 ar2 ar3 ar4 ma1 ma2 ma3 ma4
0.122 -0.256 0.423 0.681 -0.088 0.235 -0.381 -0.694
s.e. 0.190 0.156 0.154 0.190 0.180 0.147 0.146 0.177

sigma^2 estimated as 3.88: log likelihood=-3701
AIC=7421 AICc=7421 BIC=7470

　　auto.arima函数自动选出的模型为不包含常数项的ARMA(4,4)模型。

a44.Rt=arima(Rt,order=c(4,0,4),include.mean=F)
arima.test(a44.Rt)
<div class="md-section-divider"></div>

arima(x = Rt, order = c(4, 0, 4), include.mean = F)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ar1 0.1220 0.1899 0.6426 0.5205
ar2 -0.2561 0.1559 -1.6426 0.1005
ar3 0.4234 0.1544 2.7418 0.0061
ar4 0.6807 0.1899 3.5848 0.0003
ma1 -0.0877 0.1803 -0.4863 0.6268
ma2 0.2349 0.1468 1.5994 0.1097
ma3 -0.3810 0.1458 -2.6132 0.0090
ma4 -0.6939 0.1768 -3.9256 0.0001

sigma^2 estimated as 3.88 :log likelihood = -3701 aic = 7421

　　比较前面建立的AR(4)和MA(4)模型，可以看出，ARMA(4,4)模型是一个较为理想的模型，
$(1-0.1220L+0.2561L^2-0.4234L^3-0.6807L^4) R_t = (1-0.0877L+0.2349L^2-0.3810L^3-0.6939L^4)u_t$
或
$R_t = 0.1220R_{t-1} - 0.2561R_{t-2} +0.4234R_{t-3} + 0.6807R_{t-4} + u_t-0.0877u_t-1+0.2349u_{t-2}-0.3810u_t-3-0.6939u_{t-4 }$
$\sigma_u^2 = 3.88$
　　从AIC的值看，ARMA(4,4)模型的也较小。
三、ARMA模型的预测
　　若实际的数据生成过程是已知的，并且${y_t}$和${u_t}$的现在和过去各期的数值也已知，则可以以现在为原点，根据已掌握的信息，使用条件期望的方法对序列{yt}未来各期数值进行预测。预测式为：

$$E_ty_{t+j}=E(y_{t+j}|y_t,y_{t-1},y_{t-2},……,u_t,u_{t-1},u_{t-2},……)　　(4.4-8)$$
　　现在以ARMA(2,1)为例详述ARMA模型点预测。
　　对于参数都已知的ARMA(2,1)模型：
$$y_t=\phi_1y_{t-1}+\phi_2y_{t-2}+u_t+\theta{u_{t-1}}$$
　　前推j期，有：
$$y_{t+j} = \phi_1y_{t+j-1}+\phi_2y_{t+j-2}+u_{t+j}+\theta{u_{t+j-1}}$$
　　在t时已知信息条件下，求期望，得：
　　向前j步预测：
$$E_ty_{t+j}=\phi_1E_ty_{t+j-1}+\phi_2E_ty_{t+j-2} j≥2$$
　　j步预测误差：
$$e_{t(j)} = y_{t+j}-E_ty_{t+j}=\phi_1e_{t(j-1)}+\phi_2e_{t(j-2)}y_{t+j}+u_{t+j}+\theta{u_{t+j-1}}$$
　　在${ε_t}$服从正态分布的假设下，对于给定的置信水平$1-\alpha$，可建立$y_{t+j}$的预测区间为：
$$E_ty_{t+j}±Z{а/2}[var(e_{t(j)})]^{1/2}$$

predict(a44.Rt,n.ahead=5)  ＃向前预测五期
<div class="md-section-divider"></div>

$pred
Time Series:
Start = 1766
End = 1770
Frequency = 1
[1] 0.08948 -0.17490 -0.11923 0.25351 0.04832

$se
Time Series:
Start = 1766
End = 1770
Frequency = 1
[1] 1.970 1.971 1.971 1.972 1.972

　　也就是说，未来五天的预期对数收益率分别为0.08948，-0.17490，-0.11923，0.25351 0.04832。
四、ARMA模型的计算步骤：
　　1. 预处理
　　运用博克斯—詹金斯法的前提条件是：作为分析对象的时间序列是一组零均值的平稳随机序列。平稳随机序列的统计特性不随时间的推移而变化。直观地说，平稳随机序列的折线图无明显的上升或下降的趋势。但是，大量的社会经济现象随时间的推移，总表现出某种上升或下降的趋势，构成非零均值的非平稳的时间序列。对此的解决方法是在应用模型之前，对时间序列先进行零均值化和差分平稳化处理。
　　2. 模型的选择
　　首先，作出相应的时间序列图，判断时间序列有无趋势、异常点、缺失点和结构变化。如存在趋势，则需进行差分；异常点需要修正或者去除等。然后作样本自相关函数(ACF)图和偏自相关函数(PACF)图，并与理论ARMA模型的自相关函数(ACF)和偏自相关函数(PACF)图进行比较，选择可能合适的模型。选择可能合适的模型以后对模型进行估计和刷选。首先，用OLS、MLE等估计各模型参数，然后按照吝啬(principle of parsimony)原则在所估计出的模型中选出最终使用的模型，其中模型选择的准则可以采用AIC(Akaike information criterion)准则和BIC(Bayesian information criterion)准则。
　　3. 模型的定阶
　　博克斯—詹金斯法是以时间序列的自相关分析为基础的，以便识别时间序列的模式，实现建模的任务。自相关分析就是对时间序列求其本期与不同滞后期的一系列自相关系数和偏自相关系数，据以识别时间序列的特性。偏自相关系数被用来配合自相关系数，共同辨认适当的ARIMA模型。在自回归模型的识别中，我们可以使用偏自相关系数来初步判定模型的阶数；在移动平均模型中，我们可以用自相关系数来识别移动平均模型。
　　4. 模型参数的估计。
　　（1） p阶自回归模型AR(p)参数估计：可以使用Yule-Walker法、OLS或极大似然法(MLE)来估计，其中AR(p)的OLS和MLE估计是相同的。在R语言，已有基于(条件)OLS和MLE的函数arima。
　　（2） 阶移动平均模型MA(q)可用矩估计法，由方程组

$$r_k=\sigma^2(\theta_0.\theta_k+ \theta_1.\theta_{k+1}+...+\theta_q.\theta_q),0\leq{k}\leq{q}$$
求得参数 $$\theta_1,\theta_2,...,\theta_q$$ 的估计值。R语言的函数arima可以实现MA(q)模型的估计。
　　（3）ARMA(p，q)模型可用极大似然方法或者条件最小平方和法估计。R语言的函数arima同样可以解决ARMA(p,q)模型的估计问题。
　　5. 模型检验。
　　对所建立的ARMA模型优劣的检验，是通过对原始时间序列与所建立的ARMA模型之间的误差序列$e_t$ 进行检验来实现的。若误差序列$e_t$ 具有随机性，就意味着所建立的模型已包含了原始时间序列的所有趋势(包括周期性变动)，从而说明所建立的模型是对原始时间序列的合适的描述。反之，若误差序列 不具有随机性，则说明所建立的模型还有改善的余地，应对模型进行修正或重新建立新的模型。误差序列的这种随机性可利用博克斯—皮尔斯Q统计量来进行检验。Q统计量为：
$$Q=m\sum_{k=1}{n}r_k^2$$
　　式中：m为ARMA模型中所含的最大时滞；n为时间序列观测值的个数，且$Q～\chi^2(m)$ 。对于给定的置信度$1-\alpha$，将计算的Q值与$\chi^2(m)$ 比较：若$Q\leq\chi^2(m)$ ，则判定所选用的ARMA模型是合适的，若 $Q>\chi^2(m)$，则判定选用的ARMA模型不适合的。本文在R语言的基础上给出了相关的函数arima.test()。

4.5 分布滞后与自回归模型

　　在现实经济活动中，滞后现象是普遍存在的。比如货币供给的变化对经济影响很大，货币政策总是备受关注。货币政策的影响效应存在着时间上的滞后。在货币政策的传导过程中，货币扩张首先促使利率降低，或者一般价格水平的上升，这需要一段时间。这些因素对以GDP为代表的经济增长的影响，更是需要一段时间才能显示出来。只有经过一段时间以后，支出对利率的反应增强，投资、进出口和消费才会不断上升，货币政策才最终促使GDP增加。通常，货币扩张对GDP影响的最高点可能是在政策实施以后的一到两年间达到。这就要求我们在做经济分析时应该考虑时滞的影响。怎样才能把这类时间上滞后的经济关系纳入计量经济模型呢？

4.5.1 滞后效应与滞后变量模型

一、经济活动中的滞后现象
　　解释变量与被解释变量的因果联系不可能在短时间内完成，在这一过程中通常都存在时间滞后，也就是说解释变量需要通过一段时间才能完全作用于被解释变量。
　　此外，由于经济活动的惯性，一个经济指标以前的变化态势往往会延续到本期，从而形成被解释变量的当期变化同自身过去取值水平相关的情形。
　　这种被解释变量受自身或其它经济变量过去值影响的现象称为滞后效应。
　　【例4.5-1】 消费滞后。
　　消费者的消费水平，不仅依赖于当年的收入，还同以前的收入有关。一般来说，消费者不会把当年的收入全部花完。假定消费者将每一年40%的收入用于当年的花费，30%用于第二年的花费，20%用于第三年的花费，其余用作长期储蓄。于是改消费者的消费函数就可以表示成

$$Y_t = \alpha + 0.4 X_t + 0.3 X_{t-1}+ 0.2X_{t-2 }+ u_t　　(4.5-1)$$
　　其中，$Y_t、X_t$分别为第t年的消费和收入，$\alpha$ 为常数。
　　产生滞后效应的原因：
　　（1）心理因素：人们的心理定势，行为方式滞后于经济形势的变化，如中彩票的人不可能很快改变其生活方式。
　　（2）技术原因：如当年的产出在某种程度上依赖于过去若干期内投资形成的固定资产。
　　（3）制度原因：如定期存款到期才能提取，造成了它对社会购买力的影响具有滞后性。
二、滞后变量模型
　　滞后变量：是指过去时期的、对当前被解释变量产生影响的变量。滞后变量分为滞后解释变量与滞后被解释变量。
　　把滞后变量引入回归模型，这种回归模型称为滞后变量模型。
　　滞后变量模型的一般形式为
$$Y_t = \alpha +\beta_0X_t+\beta_1X_{t-1}+\beta_2X_{t-2}+...+ \beta_sX_{t-s}+\gamma_1Y_{t-1}+\gamma_2Y_{t-2}+...+\gamma_qY_{t-q}+u_t 　　 (4.5-2)$$
　　其中s,q分别为滞后解释变量和滞后被解释变量的滞后期长度。
　　滞后变量模型考虑了时间因素的作用，使静态分析的问题有可能成为动态分析。含有滞后解释变量的模型，又称动态模型（Dynamical Model）。
　　1. 分布滞后模型
　　模型中没有滞后被解释变量，仅有解释变量X的当期值及其若干期的滞后值。被解释变量受解释变量的影响分布在解释变量不同时期的滞后值上，即模型形如
$$Y_t = \alpha +\beta_0X_t+\beta_1X_{t-1}+\beta_2X_{t-2}+...+ \beta_sX_{t-s}+ u_t (4.5-3)$$
　　具有这种滞后分布结构的模型称为分布滞后模型，其中s为滞后长度。根据滞后长度s取为有限和无限，模型分别称为有限分布滞后模型和无限分布滞后模型。
　　在分布滞后模型中，各系数体现了解释变量的各个滞后值对被解释变量的不同影响程度，即通常所说的乘数效应：
　　$\beta_0$：称为短期乘数(short-run)或即期乘数(impact multiplier)，表示本期X变动一个单位对Y值的平均影响大小；
　　 $\beta_i（i=1,2…,s）$：称为延迟乘数或动态乘数，表示过去各时期X变动一个单位对Y值的平均影响大小；
　　$\sum_{i=0}^{s}\beta_i$：称为长期乘数（long-run）或总分布乘数（total distributed-lag multiplier），表示X变动一个单位时，由于滞后效应而形成的对Y总的影响大小。

　　2. 自回归分布滞后模型
　　如果滞后变量模型的解释变量仅包括自变量X的当期值和被解释变量的若干期滞后值，即模型形如

$$Y_t = \alpha +\beta_0X_t+\gamma_1Y_{t-1}+\gamma_2Y_{t-2}+...+\gamma_qY_{t-q}+u_t 　　(4.5-4)$$
　　则称这类模型为自回归分布滞后模型（autoregressive distributed lag model, ADL）：既含有Y对自身滞后变量的回归，还包括着X分布在不同时期的滞后变量，其中q称为自回归模型的阶数。
　　（1）有限自回归分布滞后模型：滞后期长度有限
　　（2）无限自回归分布滞后模型：滞后期无限，
模型
$$Y_t = \alpha_0+ \alpha_1X_t+ \alpha_2Y_{t-1}+u_t　　 (4.5-5)$$
称为一阶自回归分布滞后模型。

4.5.2 分布滞后模型的参数估计

一、分布滞后模型估计的困难
　　1. 无限期的分布滞后模型，由于样本观测值的有限性，使得无法直接对其进行估计。
　　2. 有限期的分布滞后模型，OLS会遇到如下问题：
　　（1）没有先验准则确定滞后期长度，即滞后长度难于确定的问题；
　　（2）如果滞后期较长，将缺乏足够的自由度进行估计和检验，即自由度问题；
　　（3）同名变量滞后值之间可能存在高度线性相关，即模型存在多重共线性。
　　处理的基本思想：人们提出了一系列的修正估计方法，但并不很完善。各种方法的基本思想大致相同：都是通过对各滞后变量加权，组成线性合成变量而有目的地减少滞后变量的数目，以缓解多重共线性，保证自由度。对于有限分布滞后模型，其基本思想是设法有目的地减少需要直接估计的模型参数个数，以缓解多重共线性，保证自由度。对于无限分布滞后模型，主要是通过适当的模型变换，使其转化为只需估计有限个参数的自回归模型。
二、经验加权估计法
　　所谓经验加权估计法，是根据实际经济问题的特点及经验判断，对滞后变量赋予一定的权数，利用这些权数构成各滞后变量的线性组合，以形成新的变量，再应用最小二乘法进行估计。
　　常见的滞后结构类型：递减滞后结构；$\Lambda$型滞后结构；不变滞后结构；
　　例如下面三组权数分别代表了递减滞后结构、$\Lambda$型滞后结构和不变滞后结构：
　　（1）1，1/2，1/4，1/8；
　　（2）1/4，1/2，2/3，1/4；
　　（3）1/4，1/4，1/4，1/4。
　　优点：简单易行、不损失自由度、避免多重共线性干扰及参数估计具有一致性。
　　缺点：设置权数的主观随意性较大，要求分析者对实际问题的特征有比较透彻的了解。通常的做法是，依据先验信息，多选几组权数分别估计多个模型，然后根据可决系数、F检验值、t检验值、估计标准误以及$R^2$值，从中选出最佳估计方程。
　　由于权重设置较难，所以实际中操作有一定的难度。
　　【例4.5-2】经验加权估计法：我们知道国内生产总值（GDP）与全社会固定资产投资总额（INV）有关，而且INV对GDP有明显的滞后效应。为表述方便，记$X_t=log(GDP)，Y_t=log(INV)$。
　　设定有限分布滞后模型为：

$$Y_t = \alpha +\beta_0X_t+\beta_1X_{t-1}+\beta_2X_{t-2}+ \beta_3X_{t-3}+u_t$$
　　首先我们用传统的最小二乘法进行分布滞后模型的估计。

Xt=log(GDP); Yt=log(INV)
summary(lm(Yt~Xt))
<div class="md-section-divider"></div>

Call: lm(formula = Yt ~ Xt)

Residuals:
Min 1Q Median 3Q Max
-0.18368 -0.09528 -0.00944 0.11965 0.19451

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.0167 0.0551 36.6 <2e-16 ＊＊＊
Xt 0.8143 0.0101 80.7 <2e-16 ＊＊＊
－－－
Signif.codes:0 ‘＊＊＊’ 0.001 ‘＊＊’ 0.01‘＊’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.117 on 34 degrees of freedom
Multiple R-squared: 0.995, Adjusted R-squared: 0.995
F-statistic: 6.51e+03 on 1 and 34 DF, p-value: <2e-16

summary(lm(Yt~Xt+L_(Xt)))
<div class="md-section-divider"></div>

Call: lm(formula = Yt ~ Xt + L_(Xt))

Residuals:
Min 1Q Median 3Q Max
-0.19229 -0.08511 0.00052 0.09734 0.14423

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.1027 0.0579 36.31 <2e-16 ＊＊＊
Xt 0.3316 0.1753 1.89 0.068 .
L_(Xt) 0.4836 0.1765 2.74 0.010 ＊＊
－－－
Signif. codes:0 ‘＊＊＊’ 0.001 ‘＊＊’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.106 on 32 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.996, Adjusted R-squared: 0.995
F-statistic: 3.63e+03 on 2 and 32 DF, p-value: <2e-16

　　从模型$Y_t ~ X_t +X_{t-1}$可知，模型$Y_t~X_t$有可能是一个伪回归（实际中也许的确如此）。
　　进一步从三阶滞后回归也可以看出，模型Yt~Xt有可能是一个伪回归。

summary(lm(Yt~Xt+L_(Xt,1)+L_(Xt,2)))
<div class="md-section-divider"></div>

Call: lm(formula = Yt ~ Xt + L_(Xt, 1) + L_(Xt, 2))

Residuals:
Min 1Q Median 3Q Max
-0.19268 -0.07633 0.00325 0.08862 0.14209

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.1427 0.0605 35.40 <2e-16 ＊＊＊
Xt 0.4649 0.1980 2.35 0.026＊
L_(Xt, 1) 0.0158 0.3403 0.05 0.963
L_(Xt, 2) 0.3344 0.1968 1.70 0.100 .
－－－
Signif.codes:0 ‘＊＊＊’ 0.001 ‘＊＊’ 0.01 ‘＊’ 0.05 ‘.’ 0.1 ‘’ 1

Residual standard error: 0.103 on 30 degrees of freedom
(2 observations deleted due to missingness)
Multiple R-squared: 0.996, Adjusted R-squared: 0.995
F-statistic: 2.38e+03 on 3 and 30 DF, p-value: <2e-16

summary(lm(Yt~Xt+L_(Xt,1)+L_(Xt,2)+L_(Xt,3)))
<div class="md-section-divider"></div>

Call: lm(formula = Yt ~ Xt + L_(Xt, 1, T) + L_(Xt, 2, T) + L_(Xt, 3, T))

Residuals:
Min 1Q Median 3Q Max
-0.1882 -0.0795 0.0120 0.0824 0.1588

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.1701 0.0684 31.75 <2e-16 ＊＊＊
Xt 0.4153 0.2096 1.98 0.057 .
L_(Xt, 1, T) 0.1246 0.3880 0.32 0.751
L_(Xt, 2, T) 0.1578 0.3861 0.41 0.686
L_(Xt, 3, T) 0.1162 0.2062 0.56 0.578
－－－
Signif.codes:0 ‘＊＊＊’ 0.001 ‘＊＊’ 0.01 ‘＊’ 0.05 ‘.’0.1 ‘ ’ 1

Residual standard error: 0.105 on 28 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.996, Adjusted R-squared: 0.995
F-statistic: 1.57e+03 on 4 and 28 DF, p-value: <2e-16

　　从结果可以看出，虽然模型的F值检验非常显著（p<0.001），可决系数$R^2$也很大（0.996），
　　但所有变量，包括滞后变量都不显著。
　　下面我们运用经验加权法，选择下列三组权数：
　　（1）1，1/2，1/4，1/8
　　（2）1/4，1/2，2/3，1/4 　　　　(4.5-6)
　　（3）1/4，1/4，1/4，1/4
　　分别估计上述模型，并从中选择最佳的方程。
　　记新的线性组合变量分别为：

$$\begin{matrix} Z_1 = 　　X_t + 1/2 X_{t-1} + 1/4 X_{t-2} + 1/8 X_{t-3}\\ Z_2 = 1/4 X_t + 1/2 X_{t-1} + 2/3 X_{t-2} + 1/4 X_{t-3} \\ Z_3 = 1/4 X_t + 1/4 X_{t-1} + 1/4 X_{t-2} + 1/4 X_{t-3} \end{matrix} 　　(4.5-7)$$
　　由上述公式生成线性组合变量$Z_1、Z_2、Z_3$的数据。然后分别估计如下经验加权模型。
$Y = \alpha + \beta{Z_k }+ u (k=1,2,3)$$
　　回归分析结果整理如下

Z1=Xt+1/2*L_(Xt,1)+1/4*L_(Xt,2)+1/8*L_(Xt,3)
Z2=1/4*Xt+1/2*L_(Xt,1)+2/3*L_(Xt,2)+1/4*L_(Xt,3)
Z3=1/4*Xt+1/4*L_(Xt,1)+1/4*L_(Xt,2)+1/4*L_(Xt,3)
lmZ1=lm(Yt~Z1); summary(lmZ1)
<div class="md-section-divider"></div>

Call: lm(formula = Yt ~ Z1)

Residuals:
Min 1Q Median 3Q Max
-0.18546 -0.09068 0.00824 0.08840 0.14311

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.14275 0.05473 39.1 <2e-16 ＊＊＊
Z1 0.43327 0.00527 82.2 <2e-16 ＊＊＊
－－－
Signif. codes:0‘＊＊＊’ 0.001 ‘＊＊’ 0.01 ‘＊’0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.101 on 31 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.995, Adjusted R-squared: 0.995
F-statistic: 6.75e+03 on 1 and 31 DF, p-value: <2e-16

lmZ2=lm(Yt~Z2); summary(lmZ2)
<div class="md-section-divider"></div>

Call: lm(formula = Yt ~ Z2)

Residuals:
Min 1Q Median 3Q Max
-0.1933 -0.0838 0.0210 0.0655 0.2229

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.25456 0.05742 39.3 <2e-16 ＊＊＊
Z2 0.48876 0.00639 76.5 <2e-16 ＊＊＊
－－－
Signif.codes:0 ‘＊＊＊’ 0.001 ‘＊＊’ 0.01 ‘＊’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.108 on 31 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.995, Adjusted R-squared: 0.995
F-statistic: 5.85e+03 on 1 and 31 DF, p-value: <2e-16

lmZ3=lm(Yt~Z3); summary(lmZ3)
<div class="md-section-divider"></div>

Call: lm(formula = Yt ~ Z3)

Residuals:
Min 1Q Median 3Q Max
-0.1922 -0.0839 0.0210 0.0623 0.2201

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.2424 0.0560 40.0 <2e-16 ＊＊＊
Z3 0.8154 0.0104 78.5 <2e-16 ＊＊＊
－－－

Signif. codes:0 ‘＊＊＊’0.001‘＊＊’ 0.01 ‘＊’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.106 on 31 degrees of freedom
(3 observations deleted due to missingness)
Multiple R-squared: 0.995, Adjusted R-squared: 0.995
F-statistic: 6.17e+03 on 1 and 31 DF, p-value: <2e-16

　　从上述回归分析结果可以看出，综合判断t检验值、剩余标准误、可决系数、F 检验值，可以认为：最佳的方程是模型一$Z_1$，即权数为（1,1/2,1/4,1/8）的分布滞后模型。

34.5.3 自回归模型及其估计

　　解释变量与被解释变量的因果联系不可能在短时间内完成，在这一过程中通常都存在时间滞后。由于经济活动的惯性，一个经济指标以前的变化态势往往会延续到本期，从而形成被解释变量的当期变化同自身过去取值水平相关的情形。
　　在实际问题中，有时需要对解释变量使用自回归模型进行分析。
　　一、库伊克模型
　　无限分布滞后模型中滞后项无限多，而样本观测总是有限的，因此不可能对其直接进行估计。要使模型估计能够顺利进行，必须施加一些约束或假定条件，将模型的结构作某种转化。库伊克（Koyck）变换就是其中较具代表性的方法。
　　对于如下无限分布滞后模型：

$$Y_t = \alpha+ \beta_0 X_t + \beta_1 X_{t-1 }+ \beta_2 X_t-2 + … + ut 　　 (4.5-8)$$
　　可以假定滞后解释变量$X_{t-i}$对被解释变量$Y_t$的影响随着滞后期i（i=0,1,2,…）的增加而按几何级数衰减。即滞后系数的衰减服从某种公比小于1的几何级数：
$$\beta_i = \beta_0\lambda^i (0<\lambda<1, i=0,1,2,…) 　　(4.5-9)$$
　　其中：$\beta_0$为常数，公比$\lambda$为待估参数，通常称为分布滞后衰减率，$\lambda$值越接近零，衰减速度越快。
将库伊克假定(4.5-9)式代入(4.5-8)式，得
$$Y_t=\alpha+\beta_0\sum_{i=0}^{\infty }\lambda^iX_{t-i}+u_t 　　(4.5-10)$$
　　 将(4.5-10)滞后一期，有
$$Y_{t-1}=\alpha+\beta_0\sum_{i=0}^{\infty }\lambda^{i-1}X_{t-i}+u_{t-1}　　(4.5-11)$$
　　对(4.5-11)式两边同乘$\lambda$并与(4.5-10)式相减得： $$Y_t-\lambda Y_{t-1}=(\alpha+\beta_0\sum_{i=0}^{\infty }\lambda^iX_{t-i}+u_t)-(\lambda\alpha+\beta_0\sum_{i=0}^{\infty }\lambda^iX_{t-i}+{\lambda}u_{t-1}$$

　　即

$$Y_t=\alpha(1-\lambda)+\beta_0X_t+{\lambda}Y_{t-1}+(u_t-{\lambda}u_{t-1})　　(4.5-12)$$
　　这就是库伊克模型。上述变换过程也叫库伊克变换。 令
$$\alpha^*= \alpha(1-\lambda),\beta_0^*=\beta_0$$
$$\beta_1^*=\lambda,u_t^*=u_t-{\lambda}u_{t-1}$$
　　则库伊克模型(4.3-11)式变为
$$Y_t=\alpha^*+ \beta_0^*X_t + \beta_1^*Y_{t-1} + u_t^* 　　(4.5-13)$$
　　这是一个一阶自回归模型。
　　库伊克变换的优点：
　　（1）以一个滞后被解释变量代替了大量的滞后解释变量，使模型结构得到极大简化，最大限度地保证了自由度，解决了滞后长度难以确定的问题；
　　 （2）滞后一期的被解释变量Yt-1与Xt 的线性相关程度将低于X的各滞后值之间的相关程度，从而在很大程度上缓解了多重共线性。
　　库伊克变换的缺陷：
　　（1）它假定无限滞后分布呈几何递减滞后结构。
这种假定对某些经济变量可能不适用，如固定资 产投资对总产出影响的滞后结构就不是这种类型。
　　（2）库伊克模型的随机扰动项形如
$$u_t^*=u_t-{\lambda}u_{t-1}$$
　　说明新模型的随机扰动项存在一阶自相关，且与解释变量相关。
　　（3）将随机变量作为解释变量引入了模型，不一定符合基本假定。
　　（4）库伊克变换是纯粹的数学运算结果，缺乏经济理论依据。
这些缺陷，特别是第二个缺陷，将给模型的参数估计带来定困难。
一个无限期分布滞后模型可以通过科伊克变换转化为自回归模型。事实上，许多滞后变量模型都可以转化为自回归模型，自回归模型是经济生活中更常见的模型。以自适应预期模型以及局部调整模型为例进行说明。
　　二、自适应预期模型
　　某些经济变量的变化会或多或少地受到另一些经济变量预期值的影响。为了处理这种经济现象，可以将解释变量预期值引入模型建立“期望模型”。
　　自适应预期假定：经济活动主体对某经济变量的预期，是通过一种简单的学习过程而形成的，其机理是，经济活动主体会根据自己过去在作预期时所犯错误的程度，来修正他们以后每一时期的预期，即按照过去预测偏差的某一比例对当前期望进行修正，使其适应新的经济环境。
在某些实际问题中，因变量$Y_t$并不取决于解释变量的当前实际值$X_t$，　　而取决于$X_t$的“预期水平”或“长期均衡水平”$X_t^e$。
例如，家庭本期消费水平，取决于本期收入的预期值；市场上某种商品供求量，决定于本期该商品价格的均衡值。因此，自适应预期模型最初表现形式是
$$Y_t = \beta_0 + \beta_1 X_t^e + u_t　　 (4.5-14)$$
　　 其中，$Y_t$为被解释变量，$X_t^e$为解释变量预期值，$u_t$为随机扰动项。
　　由于预期变量是不可实际观测的，往往作如下自适应预期假定:
$$X_t^e-X_{t-1}^e =r( X_t- X_{t-1}^e) 　　 (4.5-15)$$
　　其中：r为预期系数（coefficient ofexpectation）,$0{\leq}r{\leq}1$。
　　 该式的经济含义为：“经济行为者将根据过去的经验修改他们的预期”，即本期预期值的形成是一个逐步调整过程，本期预期值的增量是本期实际值与前一期预期值之差的一部分，其比例为r。这个假定还可写成：
$$X_t^e= rX_t+(1-r) X_{t-1}^e (4.5-15)$$
　　带入整理得
$$Y_t =\beta_0r+\beta_1rX_t+(1-r)Y_{t-1}+v_t　　(4.5-16)$$
　　其中，$v_t = u_t –(1-r)u_{t-1}$。
　　可见自适应预期模型转化为一阶自回归模型。
　　三、局部调整模型
　　在经济活动中，会遇到为了适应解释变量的变化，被解释变量有一个预期的最佳值与之对应的现象。
　　例如，企业为了确保生产或供应，必须保持一定的原材料储备，对应于一定的产量或销售量，存在着预期最佳库存量；为了确保一国经济健康发展，中央银行必须保持一定的货币供应，对应于一定的经济总量水平，应该有一个预期的最佳货币供应量。
　　局部调整模型的最初形式为
$$Y_t^e = \beta_0 +\beta_1X_t + u_t　　 (4.5-17)$$
　　这里，$Y_t^e$不可观测。
　　例如，企业为了保证生产和销售，必须保持一定的原材料储备。对应于一定的产量或销售量$X_t$，存在着预期的最佳库存$Y_t^e$。由于生产条件的波动，生产管理方面的原因，库存储备$Y_t$的实际变化量只是预期变化的一部分。
　　由于技术、制度、市场以及管理等各方面的限制，被解释变量的预期水平在单一周期内一般不会完全实现，而只能得到部分的调整。局部调整假设认为，被解释变量的实际变化仅仅是预期变化的一部分，即储备按预定水平逐步进行调整，故有如下局部调整假设：
$$Y_t-Y_{t-1}=\delta (Y_t^e-Y_{t-1} 或 Y_t = \delta Y_t^e+(1- \delta)Y_{t-1}$$
　　其中，$\delta$为调整系数，$0{\leq}\delta{\leq}1$，它代表调整速度。 越接近1，表明调整到预期最佳水平的速度越快。将上式代入$Y_t^e = \beta_0 +\beta_1X_t + u_t$得
$$Y_t= \delta \beta_0+ \delta \beta_1X_t +(1-\delta)Y_{t-1} +\delta{u_t } 　　(4.5-18)$$
　　满足局部调整假设的模型称为局部调整模型（Partial adjustment model）。在局部调整假设下，经过变形，局部调整模型可转化为一阶自回归模型。
　　四、三种模型的评价
　　1. 相同点
　　库伊克模型、自适应预期模型与局部调整模的最终形式都是一阶自回归模型，这样，对这三类模型的估计就转化为对相应一阶自回归模型的估计。
　　2. 区别
　　（1）导出模型的经济背景与思想不同，库伊克模型是在无限分布滞后模型的基础上根据库伊克几何分布滞后假定而导出的；自适应预期模型是由解释变量的自适应过程而得到的；局部调整模型则是对被解释变量的局部调整而得到的。
　　（2）由于模型的形成机理不同而导致随机误差项的结构有所不同,这一区别将对模型的估计带来一定影响。
　　五、自回归模型的估计
　　对于自回归模型
　 $$Y_t = \alpha+ \beta_0 X_t +\beta_1 X_{t-1} + \beta_2 X_{t-2} + … + \beta_q X_{t-q}+ u_t 　　 (4.5-19)$$
　　估计时的主要问题：滞后被解释变量的存在可能导致它与随机扰动项相关，以及随机扰动项出现序列相关性。因此，对自回归模型的估计主要需视滞后被解释变量与随机扰动项的不同关系进行估计。
　　1. 自回归模型估计的困难
　　 库伊克模型、自适应预期模型与局部调整模型，在模型结构上最终都可表示为一阶自回归。因此，对这三个模型的估计就转化为对一阶自回归模型的估计。
但是，上述一阶自回归模型的解释变量中含有滞后被解释变量$Y_{t-1}$，而$Y_{t-1}$是随机变量，它可能与随机扰动项相关；而且随机扰动项还可能自相关。模型可能违背古典假定，从而给模型的估计带来一定困难。
如果用最小二乘法直接估计自回归模型，则估计可能是有偏的，而且不是一致估计。 估计自回归模型需要解决两个问题：
　　（1）设法消除$Y_{t-1}$与$u_t$的相关性；
　　（2）检验$u_t$是否存在自相关。
　　2. 工具变量法
　　 所谓工具变量法，就是在进行参数估计的过程中选择适当的工具变量，代替回归模型中同随机扰动项存在相关性的解释变量。工具变量的选择应满足如下条件：
　　（1）与所代替的解释变量高度相关；
　　（2）与随机扰动项不相关；
　　（3）与其它解释变量不相关，以免出现多重共线性。
　　对于一阶自回归模型
$$Y_t = \alpha+ \beta_0 X_t + \beta_1Y_{t-1} + u_t$$
　　若$Y_{t-1}$与$u_t$同期相关，则OLS估计是有偏的，并且不是一致估计。
　　因此，对上述模型，通常采用工具变量法，即寻找一个新的经济变量$Z_t$，用来代替$Y_{t-1}$。 参数估计量具有一致性。
　　在实际估计中，一般用X的若干滞后的线性组合作为$Y_{t-1}$的工具变量:

　　由于原模型已假设随机扰动项ut与解释变量X及其滞后项不存在相关性，因此上述工具变量与$u_t$不再线性相关。一个更简单的情形是直接用$X_{t-1}$作为$Y_{t-1}$的工具变量。
　　若滞后被解释变量$Y_{t-1}$与随机扰动项ut同期无关（如局部调整模型），可直接使用OLS法进行估计，得到一致估计量。
上述工具变量法只解决了解释变量与$u_t$相关对参数估计所造成的影响，但没有解决$u_t$的自相关问题。
　　【例4.5-1】工具法求解自回归分布滞后模型。

Yt=log(TAX); Xt=log(GDP)
lm4.5=lm(Yt~Xt+L_(Xt))
summary(lm4.5)
<div class="md-section-divider"></div>

Call: lm(formula = Yt ~ Xt + L_(Xt))

Residuals:
Min 1Q Median 3Q Max
-0.4314 -0.0975 -0.0065 0.1368 0.5348

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.161 0.188 -11.51 6.5e-13 ＊＊８
Xt 0.882 0.653 1.35 0.19
L_(Xt) 0.155 0.651 0.24 0.81
－－－
Signif.codes:0 ‘＊＊＊’ 0.001‘＊＊＊’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.22 on 32 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.983, Adjusted R-squared: 0.982
F-statistic: 911 on 2 and 32 DF, p-value: <2e-16

　　滞后一期后X不显著，下面用工具变量法进来局部调整模型

summary(lm(Yt~Xt+L_(Yt)))
<div class="md-section-divider"></div>

Call: lm(formula = Yt ~ Xt + L_(Yt))

Residuals:
Min 1Q Median 3Q Max
-0.0944 -0.0508 -0.0281 0.0244 0.6085

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.2996 0.2339 -1.28 0.209
Xt 0.1905 0.0986 1.93 0.062 .
L_(Yt) 0.8219 0.0948 8.67 6.6e-10 ＊＊＊
－－－
Signif.codes:0 ‘＊＊＊’ 0.001 ‘＊＊’ 0.01 ‘＊’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.12 on 32 degrees of freedom
(1 observation deleted due to missingness)
Multiple R-squared: 0.995, Adjusted R-squared: 0.995
F-statistic: 3.08e+03 on 2 and 32 DF, p-value: <2e-16

　　事实上，对于自回归模型，ut项的自相关问题始终存在，对于此问题，至今没有完全有效的解决方法。唯一可做的，就是尽可能地建立“正确”的模型，以使序列相关性的程度减轻。
　　但本书通过研究发现，可通过ARMA模型来解决这类问题。

a.Yt=arima(Yt,order=c(1,0,0),xreg=Xt)
arima.test(a.Yt)

arima(x = Yt, order = c(1, 0, 0), xreg = Xt)

Coefficients:
Estimate Std. Error t value Pr(>|t|)
ar1 0.8392 0.08309 10.101 5.484e-24
intercept -2.1467 0.32861 -6.533 6.459e-11
Xt 1.0369 0.05059 20.497 2.278e-93

sigma^2 estimated as 0.0122 :log likelihood = 27.62 aic = -47.24

　　显然，该模型要比上面的模型好。

练习题4

一、手工解答题（纸质作业）
　　1. 举例说明经济现象中的自相关性的存在。什么是一阶自相关和高阶自相关？
　　2. 自相关与偏自相关有何不同？
　　3. 序列自相关图示法的判断依据是什么？
　　4. 产生分布滞后效应的原因有哪些？分布滞后模型的参数估计有哪些困难？
　　5. AR模型的平稳性条件是什么？
　　6. 伊库克模型、自适应预期模型和局部调整模型的相同点与不同点有哪些？
　　7. 如何根据自相关函数和偏相关系函数初步判断某个平稳随机过程为AR、MA、ARMA。
　　8. ARMA模型的分析步骤有哪些？
二、软件计算题（电子作业）
　　1. 下表数据是1970-1991年美国制造业固定厂房设备投资Y和销售量X，以10亿美元计价，且经过季节调整，根据该数据，判断厂房开支和销售量序列是否平稳？

　　２. 表中给出了某地区1980-2001年固定资产投资Y与销售额X的资料（单位：亿元）。

　　试就下列模型，按照一定的处理方法估计模型参数，并解释模型的经济意义，探测模型扰动项的一阶自相关性。
　　（1）设定模型 ，运用局部调整假定。
　　（2）设定模型 ，运用局部调整假定。
　　（3）设定模型 ，运用自适应预期假定。
　　（4）运用阿尔蒙多项式变换法，估计分布滞后模型：

　　３. 表中给出了某地区消费总额Y（亿元）和货币收入总额X（亿元）的年度资料，

　　分析该地区消费同收入的关系
　　（1）做Y关于X的回归，对回归结果进行分析判断；
　　（2）建立分布滞后模型，用库伊克变换转换为库伊克模型后进行估计，并对估计结果进行分析判断；
　　（3）建立局部调整——自适应期望综合模型进行分析。
　　４. 对全国居民消费价格指数进行分析。请读者选取1999年1月至2013年12月的全国居民消费价格指数CPI（月度数据，上年同月=100）作为样本数据，用R语言命令建立相应的AR模型、MA模型和ARM模型，并从中选一个合适的模型。
　　５. 股票收益率的研究。请读者选取2004年1月1日至2013年12月31日的沪深300指数作为样本数据，对我国证券去市场沪深300股票指数收益率的变动进行分析。并用R语言命令建立相应的AR模型、MA模型和ARM模型，并从中选一个合适的模型。