ARMA模型和SAS求解.doc

第6讲 时间序列分析教材：应用时间序列分析课件（中国人民大学 王燕），SAS如何解及下载例程。时间序列分析(Time series analysis)是一种动态数据处理的统计方法。该方法基于随机过程理论和数理统计学方法，研究随机数据序列所遵从的统计规律，以用于解决实际问题。时间序列是把反映现象发展水平的统计指标数值，按照时间先后顺序排列起来所形成的一组统计数字序列。时间序列又称动态数列或时间数列。时间序列分析就是利用这组数列，应用数理统计方法加以处理，以预测未来事物的发展。时间序列分析是定量预测方法之一，它的基本原理：一是承认事物发展的延续性。应用过去数据，就能推测事物的发展趋势。二是考虑到事物发展的随机性。任何事物发展都可能受偶然因素影响，为此要利用统计分析中加权平均法对历史数据进行处理。该方法简单易行，便于掌握，但准确性差，一般只适用于短期预测。时间序列预测一般反映三种实际变化规律:趋势变化、周期性变化、随机性变化。 时间序列分析是根据系统观测得到的时间序列数据，通过曲线拟合和参数估计来建立数学模型的理论和方法。时间序列分析常用在国民经济宏观控制、区域综合发展规划、企业经营管理、市场潜量预测、气象预报、水文预报、地震前兆预报、农作物病虫灾害预报、环境污染控制、生态平衡、天文学和海洋学等方面。 时间序列分析主要用途:系统描述。根据对系统进行观测得到的时间序列数据，用曲线拟合方法对系统进行客观的描述。系统分析。当观测值取自两个以上变量时，可用一个时间序列中的变化去说明另一个时间序列中的变化，从而深入了解给定时间序列产生的机理。预测未来。一般用ARMA模型拟合时间序列，预测该时间序列未来值。决策和控制。根据时间序列模型可调整输入变量使系统发展过程保持在目标值上，即预测到过程要偏离目标时便可进行必要的控制。 基本步骤:用观测、调查、统计、抽样等方法取得被观测系统时间序列动态数据。根据动态数据作相关图，进行相关分析，求自相关函数。相关图能显示出变化的趋势和周期，并能发现跳点和拐点。跳点是指与其他数据不一致的观测值。如果跳点是正确的观测值,在建模时应考虑进去,如果是反常现象，则应把跳点调整到期望值。拐点则是指时间序列从上升趋势突然变为下降趋势的点。如果存在拐点，则在建模时必须用不同的模型去分段拟合该时间序列，例如采用门限回归模型。辨识合适的随机模型,进行曲线拟合,即用通用随机模型去拟合时间序列的观测数据。对于短的或简单的时间序列，可用趋势模型和季节模型加上误差来进行拟合。对于平稳时间序列，可用通用ARMA模型（自回归滑动平均模型）及其特殊情况的自回归模型、滑动平均模型或组合ARMA模型等来进行拟合。当观测值多于50个时一般都采用ARMA模型。对于非平稳时间序列则要先将观测到的时间序列进行差分运算，化为平稳时间序列，再用适当模型去拟合这个差分序列。本章重点：1）建立阶自回归模型：2）建立阶移动平均模型： 3）模型： 三个模型的拖尾、截尾性模型自相关系数偏自相关系数拖尾阶截尾阶截尾拖尾拖尾拖尾建模步骤：平稳非白噪声序列？计算ACF,PACFARMA模型识别估计模型中未知参数值模型优化模型预测模型检验即残差白噪声检验NY图6.0 自回归滑动平均(ARMA)模型建模步骤Y平滑处理N一. 几个概念随机过程: X(t);-t,其中X(t)是随机变量.随机序列: Xk;k=,-1,0,1,其中Xk是随机变量.特别当Xk=X(kh)时,序列Xk是过程X(t)的等间隔采样序列. 根据随机变量X和它的样本的定义, 我们有:样本序列:,x-1, x 0,x1,是序列xk的一个样本序列,又称为一个实现, 或一个观测序列等.请注意: 随机变量X的一个样本,就是一个数; 随机向量X的一个样本,就是一个向量数; 随机序列Xk的一个样本, 是一个无穷数列;在实际应用中, 我们无法记录无穷数列,从而在讨论随机序列Xk的样本时, 只能考虑一个样本的有限部分, 比如x1,x2,xn是序列Xk的一段观测值序列.在理论讨论时,为了方便又不得不涉及无穷数列. 这些都是学习和掌握时间序列分析时, 首先要认清的起点.序列的分布 :回忆随机变量X的定义便知,它的特征被它的概率分布所确定. 同样, 随机序列也被它的概率分布所确定.不过, 随机序列的分布是无穷个随机变量的概率分布,其复杂性可以想像. 这里为了避免涉及太深的概率论概念, 我们仅考虑最简单的,即XkN(mk,s2k), 它有密度fk(x)=(2ps2k)-1/2exp(x-mk)2/2s2k而且(Xk+1,Xk+2,Xk+m)有联合正态分布. 于是有:期望(均值): EXk=xfk(x)dx=mk,方差: Var(Xk)=E(Xk-mk)2=(x-mk)2fk(x)dx=s2k.自协方差: gkj=E(Xk-mk)(Xj-mj)=(x-mk)(y-mj)fkj(x,y)dxdy= E(Xj-mj)(Xk-mk)= gjk.（注E(X-EX)(Y-EY)展开=E(XY)-E(X)E(Y)=cov(X,Y)，协方差，不同事件之间的相关性度量）二. 平稳性检验1. 严平稳和宽平稳平稳时间序列有两种定义，根据限制条件的严格程度，分为:l 严平稳时间序列（strictly stationary）指序列所有的统计性质都不会随着时间的推移而发生变化，即(Xk+1,Xk+2,Xk+m)的联合分布（实际很难求）与k无关。l 宽平稳时间序列（week stationary）指序列的统计性质只要保证序列的二阶矩平稳就能保证序列的主要性质近似稳定。如果在任取时间、和时，时间序列满足如下三个条件：(6.1)(6.2)(6.3)则称为宽平稳时间序列。也称为弱平稳或二阶平稳。对于正态随机序列而言，由于联合概率分布仅由均值向量和协方差阵决定，即只要二阶矩平稳，就等于分布平稳了。2. 平稳时间序列的统计性质根据平稳时间序列的定义，可以推断出两个重要的统计性质：l 均值为常数。即式(6.2)的条件。l 自协方差只依赖于时间的平均长度（只与起点有关与跨度无关）。即式(6.3)的条件。如果定义自协方方差函数（autocovariance function）为：(6.4)那么它可由二维函数简化为一维函数，由此引出延迟自协方差函数：(6.5)容易推断出平稳时间序列一定具有常数方差：(6.6)如果定义时间序列自相关函数（autocorrelation function），简记为ACF：(6.7)由延迟自协方差函数的概念可以等价得到延迟自相关函数的概念：(6.8)容易验证自相关函数具有几个基本性质：l ；l ；l 自相关阵为对称非负定阵；l 非惟一性。注意区分：协方差函数和相关函数度量两个不同事件彼此之间的相互影响的程度。自协方差函数和自相关函数度量同一事件在两个不同时期之间的相互影响的程度。 3. 样本的估计值在平稳序列场合，序列的均值等于常数意味着原本含有可列多个随机变量的均值序列变成了只含有一个变量的常数序列，所以常数均值的估计值为(6.9)同样可以根据平稳序列二阶矩平稳的性质，得到基于样本计算出来的各种估计值。延迟自协方差函数的估计值：(6.10)总体方差的估计值：(6.11)延迟自相关函数的估计值：(6.12)4. 平稳性检验的方法对序列的平稳性检验有两种方法：一种是根据时序图和自相关图显示的特征做出判断的图检验方法；一是构造检验统计量进行假设检验的单位根检验（unit root test）方法。是建模的前提,可借助SAS(见例程)、SPSS、Eviews等软件实现。l 时序图和自相关图检验l 单位根检验（unit root test）所谓单位根检验就是通过检验时间序列自回归特征方程的特征根是在单位圆内还是在单位圆外（包括在单元圆上），来检验时间序列的平稳性。单位根检验统计量中最常用的是ADF检验统计量，又称增广DF检验（augmented Dickey-Fuller）。对任一p阶自回归AR(p)过程(6.13)它的特征方程为(6.14)如果该方程所有的特征根都在单位圆内，即则序列平稳。如果至少存在一个特征根不在单位圆内，不妨设，则序列非平稳，且自回归系数之和恰好等于1。即(6.15)因而，对于AR(p)过程可以通过检验自回归系数之和是否大于等于1来考察该序列的平稳性。设，那么原假设：（序列非平稳），ADF检验统计量：(6.16)式中，为参数的样本标准差。1979年，Dickey和Fuller使用蒙特卡洛模拟方法算出了检验统计量的临界值表。三. 纯随机性检验如果序列值彼此之间没有任何相关性，那就意味着该序列是一个没有记忆的数据序列，即过去的行为对未来的发展没有丝毫影响，这种序列我们称之为纯随机序列。从统计分析的角度而言，纯随机序列是没有任何分析价值的序列。因此，为了确保平稳序列还值不值得分析下去，需要对平稳序列进行纯随机性检验。5. 纯随机序列（白噪声序列）如果在任取时间和时，时间序列满足如下三个条件：(6.17)(6.18)(6.19)称此序列为纯随机序列，也称为白噪声（white noise）序列，简记为。之所以称之为白噪声序列是因为人们最初发现白光具有这种特性。比较平稳时间序列的定义，可看出白噪声序列一定是平稳序列，且是一种最简单的平稳序列。见图6.1所示是随机生成的1000个服从标准正态分布的白噪声序列观察值。图6.1 标准正态白噪声序列时序图根据白噪声序列的定义，白噪声序列具有三个重要的性质：l 常数均值（）;l 纯随机性（）;l 方差齐性（，即序列中每个度量的方差相等）。6. 纯随机性检验即白噪声检验(通常对残差序列做)。Barlett证明，如果一个时间序列是纯随机的，得到一个观察期数为的观察序列，那么该序列的延迟非零期的样本自相关系数将近似服从均值为零、方差为序列观察数倒数的正态分布，即(6.20)式中为延迟期数，为样本观察期数。根据Barlett定理，可以构造检验统计量和检验统计量来检验序列的纯随机性。原假设：延迟期数小于或等于期的序列值之间相互独立，即；备选假设：延迟期数小于或等于期的序列值之间有相关性，即至少存在某个。1) 检验统计量由Box和Pierce推导出的检验统计量为：(6.21)式中，为序列观察期数，为指定延迟期数。2) 检验统计量因为检验统计量在小样本场合时不太精确，所以Ljung和Box又推导出检验统计量为：(6.22)式中，为序列观察期数，为指定延迟期数。一般取值为6、12。为什么只需要检验前6期和前12期延迟的检验统计量就可以直接判断序列是否为白噪声序列呢？这是因为平稳序列通常具有短期相关性，只要序列时期足够长，自相关系数都会收敛于零。所以，如果序列值之间存在显著的相关关系，通常只存在在延迟时期比较短的序列值之间，而如果短期延迟的序列之间都不存在显著的相关关系，那么长期延迟之间就更不会存在显著的相关关系。 纯随机性检验小结：当Q统计量的P值大于给定的显著水平a（譬如=0.05）时，不拒绝原假设，即相关系数为0，认为序列是白噪声无任何信息可取，停止统计分析（建模）。四. 方法性工具7. 差分运算差分运算分为两种：步差分和阶差分。1) 步差分相距期的两个序列值之间的减法运算称为步差分运算，记为，表示与之间的减法运算，即：(6.23)2) 阶差分相距一期的两个序列值之间的减法运算称为1阶差分运算，记为，表示与之间的减法运算，即：(6.24)对1阶差分运算后序列再进行一次1阶差分运算称为2阶差分，记为，表示与之间的减法运算，即：(6.25)依此类推，对阶差分后序列再进行一次1阶差分运算称为阶差分，记为，表示与之间的减法运算，即：(6.26)8. 延迟算子延迟算子类似于一个时间指针，一个延迟算子乘以当前序列值，就相当于把当前序列值的时间向过去拨了一个时间刻度，记B为延迟算子，有(6.27)用延迟算子表示的步差分为：(6.28)用延迟算子表示的阶差分为：(6.29)五. ARMA模型ARMA模型的全称是自回归移动平均（auto regression moving average）模型，它是目前最常用的拟合平稳时间序列的模型。ARMA模型又可细分为AR模型、MA模型和ARMA模型三大类。9. 模型具有如下结构的模型称为阶自回归模型，简记为：(6.30)其中包含三个限制条件：模型的最高阶数为，即；随机干扰序列为零均值的白噪声序列，即；当期的随机干扰与过去的序列值无关，即。1) 中心化的模型当时，式(6.30)又称为中心化的模型。非中心化的序列都可以通过假设满足平稳性条件，在式(6.30)两边取期望，根据平稳时间序列均值为常数的性质，有，且因为为零均值的白噪声，有，所以：(6.31)如果把非中心化的序列减去上式(6.31)中的，则转化为中心化序列。特别地，对于中心化序列，有。引进延迟算子，设，又称为阶自回归系数多项式，则中心化模型可以简记为：(6.32)2) 模型的方差要得到平稳模型的方差，需要借助于Green函数的帮助。下面以求模型的方差为例来说明：将第二式代入第一式，有当我们继续将代入上式，一直到，可得到如果，设Green函数为，上式可改为(6.33)对求方差为(6.34)3) 模型的协方差对中心化的平稳模型在等号两边同乘，再求期望得到(6.35)由模型的限制条件，有，再根据平稳时间序列的统计性质，有自协方差函数只依赖于时间的平均长度而与时间的起止点无关，于是可由(6.35)式得到自协方差函数的递推公式：(6.36)例如，对于模型的自协方差函数的递推公式为：(6.37)4) 模型的自相关函数由于平稳时间序列有自相关函数，在自协方差函数的递推公式(6.36)等号两边同除以方差函数，就得到自相关函数的递推公式：(6.38)例如，对于模型的自相关函数的递推公式为：(6.39)根据式(6.38)可以推出，平稳模型的自相关函数有两个显著的性质：l 拖尾性指自相关函数始终有非零取值，不会在大于某个常数之后就恒等于零。l 负指数衰减随着时间的推移，自相关函数会迅速衰减，且以负指数（其中为自相关函数的差分方程的特征根）的速度在减小。见图6.2和图6.3所示是两个平稳模型的理论自相关图。图6.2 ACF按负指数单调收敛到零图6.3 ACF按正负相间地衰减到零5) 模型的偏自相关系数对于一个平稳模型，求出滞后自相关系数时，实际上得到的并不是与之间单纯的相关关系。因为这个还会受到中间个随机变量的影响，即这个随机变量既与又与具有相关关系。为了能单纯测度与之间的相关关系，引进了时间序列偏自相关函数（ partial autocorrelation function），简记为PACF。它是在剔除了中间个随机变量的干扰之后的滞后自相关系数，计算公式为：(6.40)式中，。如果我们用过去的期序列值对作阶自回归拟合，即(6.41)那么有。这说明滞后偏自相关系数实际上等于阶自回归模型第个回归系数的值。根据这个性质很容易计算PACF的值。在公式(6.41)中等号两边同乘，求期望并除以，得到(6.42)取前个方程构成的方程组：(6.43)该方程组被称为Yule-Walker方程。根据线性方程组求解的Gramer法则，有(6.44)式中：可以证明对于平稳模型，当时，有，这样。也就是说平稳模型的偏自相关系数具有步截尾性。见图6.4和图6.5所示是两个平稳模型的样本偏自相关图。图6.4 一个AR(1)模型n=101样本偏自相关函数PACF(k)图图6.5 一个AR(1)模型n=101样本偏自相关函数PACF(k)图由于样本的随机性，样本偏自相关系数不会和理论偏自相关系数一样严格截尾，但可以从图6.4和图6.5 中看出，两个平稳模型的样本偏自相关系数1阶显著不为零，1阶之后都近似为零。样本偏自相关图可以直观地验证平稳模型偏自相关系数具有步截尾性。10. 模型具有如下结构的模型称为阶移动平均，简记为： (6.45)其中包含两个限制条件：模型的最高阶数为，即；随机干扰序列为零均值的白噪声序列，即。1) 中心化的模型当时，式(6.45)又称为中心化的模型。非中心化的序列都可以通过假设满足平稳性条件，在式(6.45)两边取期望，根据平稳时间序列均值为常数的性质，有，且因为为零均值的白噪声，有，所以：(6.46)如果把非中心化的序列减去上式(6.46)中的，则转化为中心化序列。特别地，对于中心化序列，有。引进延迟算子，设，又称为阶自移动平均系数多项式，则中心化模型可以简记为：(6.47)2) 模型的方差平稳模型的方差为：(6.48)3) 模型的自协方差平稳模型的自协方差只与滞后阶数相关，且阶截尾。当时，；当时，；当时，有(6.49)4) 模型的自相关系数平稳模型的自相关系数为(6.50)5) 模型的偏自相关系数在中心化的平稳模型场合，滞后阶偏自相关系数为：(6.51)容易证明平稳模型的偏自相关系数拖尾性。见图6.6和图6.7所示是一个平稳模型的样本自相关图和样本偏自相关图。图6.6 一个MA(1)模型n=101样本自相关函数截尾图图6.7 一个MA(1)模型n=101样本偏自相关函数拖尾图6) 模型的可逆性容易验证当两个模型具有如下结构时：(6.52)根据公式(6.50)计算，它们的自相关系数正好相等。即不同的模型却拥有完全相同的自相关系数。这种自相关系数的不惟一性将会导致拟合模型和随机时间序列之间不会是一一对应关系。为了保证一个给定的自相关函数能够对应惟一的模型，我们需要给模型增加约束条件。这个约束条件称为的可逆性条件。把上式(6.52)中两个模型表示成两个自相关模型形式：(6.53)注意表示成自相关模型时运用公式，其中。显然，当时，模型1收敛，而模型2不收敛；当时，则模型2收敛，而模型1不收敛。若一个模型能够表示成收敛的模型形式，那么该模型称为可逆模型。一个自相关系数惟一对应一个可逆模型。11. 模型具有如下结构的模型称为自回归移动平均模型，简记为： (6.54)若，该模型称为中心化模型。模型的限制条件与模型、模型相同。引进延迟算子，中心化模型简记为：(6.55)式中：，称为阶自回归系数多项式，称为阶自移动平均系数多项式。显然，当时，模型就退化成模型；当时，模型就退化成模型。所以，模型和模型实际上是的特例，它们统称为模型。而模型的统计性质也正是模型和模型统计性质的有机组合。由于模型可以转化为无穷阶移动平均模型，所以模型的自相关系数不截尾。同理，由于模型也可以转化为无穷阶自回归模型，所以模型的偏自相关系数也不截尾。总结模型、模型和模型的自相关系数和偏自相关系数的规律，见表6.1所示。表6.1 拖尾性和截尾性模型自相关系数偏自相关系数拖尾阶截尾阶截尾拖尾拖尾拖尾假如某个时间序列观察值可以判定为平稳非白噪声序列，计算出样本自相关系数（ACF）和样本偏自相关系数（PACF）之后，就要根据它们表现出来的性质，选择阶数适当的模型拟合观察值序列。即根据样本的自相关系数和样本偏自相关系数性质估计自相关阶数和移动平均阶数。因此，这个过程也称为模型定阶过程或模型识别过程。由于样本的随机性，样本的自相关系数和偏自相关系数不会呈现出理论截尾的完美情况，本应截尾处仍会呈现出小值震荡的情况。同时，由于平稳时间序列通常都具有短期相关性，随着延迟阶数变大，自相关系数和偏自相关系数都会衰减至零值附近作小值波动。那么，如何判断自相关系数和偏自相关系数是截尾还是拖尾呢？以及如果为截尾那么相应的阶数为多少？ 通常分析人员是依据样本的自相关系数和偏自相关系数近似分布来作出尽可能合理的判断。Jankins和Watts已经证明样本自相关系数是总体自相关系数的有偏估计：(6.56)式中为延迟阶数，为样本容量。根据Bartlett公式计算样本自相关系数的方差近似等于：(6.57)当延迟阶数足够大时，；当样本容量充分大时，。所以样本自相关系数近似服从正态分布：(6.58)Quenouille证明，样本偏自相关系数也同样近似服从这个正态分布：(6.59)设显著水平取。如果样本自相关系数和样本偏自相关系数在最初的阶明显大于2倍标准差，而后几乎95%的系数都落在2倍标准差的范围内，且非零系数衰减为小值波动的过程非常突然，通常视为阶截尾；如果有超过5%的样本相关系数大于2倍标准差，或者非零系数衰减为小值波动的过程比较缓慢或连续，通常视为拖尾。六. 参数估计和检验对于一个非中心化，有(6.60)通过样本的自相关系数和偏自相关系数的性质，估计出自相关阶数和移动平均阶数。为模型定阶后，该模型共含有个未知参数：。参数用样本均值来估计总体均值（矩估计法）。对原序列中心化后，待估参数减少一个。对个未知参数的估计方法有三种：矩估计、极大似然估计和最小二乘估计。12. 参数的矩估计用时间序列样本数据计算出延迟1阶到阶的样本自相关系数，延迟阶的总体自相关系数为，公式中包含个未知参数变量。如果用计算出的样本自相关系数来估计总体自相系数，那么有个联立方程组：(6.61)从中解出个未知参数变量的值作为模型的参数估计值。这种方法称为参数的矩估计。白噪声序列的方差的矩估计，是用时间序列样本数据计算出样本方差来估计总体方差求得。模型的两边同时求方差，并把相应参数变量的估计值代入，可得白噪声序列的方差估计为：(6.62)13. 参数的极大似然估计当总体分布类型已知时，极大似然估计ML（maximum-likelihood）是常用的估计方法。极大似然估计的基本思想，是认为样本来自使该样本出现概率最大的总体。因此，未知参数的极大似然估计，就是使得似然函数（即联合密度函数）达到最大值的参数值。即：(6.63)在时间序列分析中，序列的总体分布通常是未知的。为了便于分析和计算，通常假设序列服从多元正态分布，它的联合密度函数是可导的。当似然函数关于参数可导时，常常可以通过求导方法来获得似然函数极大值对应的参数值。在求极大似然估计时，为了求导方便，常对似然函数取对数，然后对对数似然函数中的未知参数求偏导数，得到似然方程组。理论上，只要求解似然方程组即可得到未知参数的极大似然估计。但是在实际上是使用计算机经过复杂的迭代算法求出未知参数的极大似然估计。极大似然估计与矩估计的比较：矩估计的优点是不要求知道总体的分布，计算量小，估计思想简单直观。但缺点是只用到了样本自相关系数的信息，序列中的其他信息被忽略了，这导致矩估计方法是一种比较粗糙的估计方法，它的估计精度一般较差。因此，它常被作为极大似然估计和最小二乘估计的迭代计算的初始值。极大似然估计的优点是充分应用了每一个观察值所提供的信息，因而它的估计精度高，同时，还具有估计的一致性、渐近正态性和渐近有效性等优良统计性质，是一种非常优良的参数估计方法。14. 参数的最小二乘估计参数的最小二乘估计ULS（unconditional least squares）是使模型的残差平方和达到最小的那组参数值。即：(6.64)同极大似然估计一样，未知参数的最小二乘估计通常也是使用计算机借助迭代方法求出的。由于充分利用了序列的信息，因此最小二乘估计的精度最高。在实际运用中，最常用的是条件最小二乘估计CLS（conditional least squares）方法。它假定时间序列过去未观察到序列值等于序列均值，如果是中心化后的序列，则序列过去未观察到序列值等于零（）。根据这个假定可以的得到残差的有限项表达式：(6.65)于是残差平方和达到最小的那组参数值为：(6.66)在实际运用中，条件最小二乘估计CLS也是通过迭代法求出参数的估计值。15. 模型检验和参数检验在拟合好模型的参数之后，一般来说，都要对该拟合模型进行必要的显著性检验。包括：模型的显著性检验和参数的显著性检验。在ARMA模型场合，我们都使用统计量检验残差序列的自相关性。为了克服DW检验的有偏性，Durbin在1970年提出了DW统计量的两个修正统计量：Durbin t和Durbin h统计量，这两个统计量渐近等价。Durbin h统计量为：(6.67)式中，为观察值序列的长度；为延迟因变量系数的最小二乘估计的方差。修正后的有效地提高了检验精度，成为延迟因变量场合常用的自相关检验统计量。参数的显著性检验是要检验每一个模型参数是否显著非零。因为如果某个参数为零，模型中包含这个参数的乘积项就为零，可以简化模型。因此，这个检验的目的就是为了使模型最精简。原假设为：某个未知参数；备选假设为：。可以构造出检验未知参数显著性的检验统计量，其中为参数的个数。如果某个参数不显著，即表示所对应的那个自变量对因变量的影响不明显，该自变量就可以从拟合模型中剔除。剔除不显著参数对应的自变量后应重新拟合模型，最终模型将由一系列参数显著非零的自变量表示。七. 模型优化当一个拟合模型在指定的置信水平下通过了检验，说明了在这个置信水平下该拟合模型能有效地拟合时间序列观察值的波动。但是这种有效的拟合模型并不是惟一的。如果同一个时间序列可以构造两个拟合模型，且两个模型都显著有效，那么应该选择哪个拟合模型用于统计推断呢？通常采用AIC和SBC信息准则来进行模型优化。16. AIC准则AIC准则是由日本统计学家赤池弘次（Akaike）于1973年提出，AIC全称是最小信息量准则（an information criterion）。AIC准则是一种考评综合最优配置的指标，它是拟合精度和参数未知个数的加权函数：AIC=2ln(模型中极大似然函数值)+2(模型中未知参数个数)(6.68)使AIC函数达到最小值的模型被认为是最优模型。17. BIC准则AIC准则也有不足之处：如果时间序列很长，相关信息就越分散，需要多自变量复杂拟合模型才能使拟合精度比较高。在AIC准则中拟合误差等于，即拟合误差随样本容量放大。但是模型参数个数的惩罚因子却与无关，权重始终为常数2。因此在样本容量趋于无穷大时，由AIC准则选择的拟合模型不收敛于真实模型，它通常比真实模型所含的未知参数个数要多。为了弥补AIC准则的不足，Akaike于1976年提出BIC准则。而Schwartz在1978年根据Bays理论也得出同样的判别准则，称为SBC准则。SBC准则定义为：SBC=2ln(模型中极大似然函数值)+ln(n)(模型中未知参数个数)(6.69)它对AIC的改进就是将未知参数个数的惩罚权重由常数2变成了样本容量的对数。在所有通过检验的模型中使得AIC或SBC函数达到最小的模型为相对最优模型。之所以称为相对最优模型是因为不可能比较所有模型。八. 序列预测所谓预测就是利用时间序列已观察到的样本值对时间序列在未来某个时刻的取值进行估计。常用的预测方法是线性最小方差预测。线性是指预测值为观察值序列的线性函数，最小方差是指预测方差达到最小。根据模型的平稳性和可逆性，可以用Green函数的传递形式和逆转函数的逆转形式等价描述该序列：(6.70)(6.71)(6.70)式中，为Green函数：(6.72)式中：,(6.73)(6.71)式中，为逆转函数：(6.74)式中的和定义见(6.73)式。把式(6.71)代入式(6.70)中，可得到：(6.75)显然是历史数据的线性函数。不妨简记为：(6.76)那么，对于任意一个将来时刻而言，也可以表示成(6.76)式。但问题是已知，而未知。根据线性函数的可加性，所有未知信息都可以用已知信息的线性函数表示出来，并用该函数进行估计：(6.77)用衡量预测误差：(6.78)显然，预测的误差越小预测的精度就越高，目前最常用的预测原则是预测误差的方差最小法：(6.79)因为为的线性函数，所以也称为线性预测方差最小法。在线性预测方差最小法下得到的估计值是在序列已知的情况下得到的条件无偏最小方差估计值。且预测方差只与预测步长有关，而与预测起始点无关。但预测步长越大预测值的方差越大，因此只适合于短期预测。在正态假定下，估计值的的置信区间为：(6.80)九. proc arima过程 SAS/ETS软件中的ARIMA过程，是集模型识别、参数估计、相对最优模型选择、短期预测于一体的综合工具。一个完整的ARIMA过程是由identify(识别)、estimate(估计)、forecast(预测)三条命令组成。这三条命涵盖了平稳序列建模的每个步骤。它们既可以分开使用也可以联合使用。 ARIMA过程的语法格式如下：Proc arima data = ;Identify var = ststionarity = ( = ) nlag = k minic p = (0: n) q = (0: m);Estimate method = 参数估计方法 p = n q = m noint;Forecast lead = u id = 时间变量 out = results;Run 说明： 1) Identify命令输出5方面的信息：分析变量的描述性统计、样本自相关图、样本逆自相关图、样本偏自相关图和纯随机检验结果； 如果增加可选项minic短语，则可以得到一定范围内的最优模型定阶。如果增加可选项ststionarity短语，则可以得到单位根检验的结果。 2) Estimate命令输出5方面的信息：未知参数的估计值、拟合统计量的值、系数相关阵、残差自相关检验结果和拟合模型的具体形式； 如果增加可选项method短语，则可以选择参数估计的方法：ML（极大似然法）、ULS（最小二乘法）、CLS（条件最小二乘法系统默认）。 3) Forecast命令输出模型对序列的短期预测值，包括预测值的标准差、95%的置信下限、95%的置信上限等； 利用存储在临时数据集results中的数据，可以绘制拟合、预测图。十. ARIMAX模型的应用实例 下面以河南省历年国民生产总值数据（见表6.2）为例，介绍ARIMA模型的建模过程及SAS实现(见GDP序列arima模型程序.doc或gdpcx.sas)，并通过模型预测未来若干年内的GDP/PGDP变化情况。表6.2 河南省历年国民生产总值数据年份(Year)生产总值(亿元) (GDP)人均生产总值(元)(PGDP)年份(Year)生产总值(亿元) (GDP)人均生产总值(元)(PGDP)1978162.92232.319921279.751452.31979190.09266.719931662.761867.41980229.16316.719942224.432475.21981249.69340.119953002.743312.81982263.3035319963661.184007.41983327.95432.919974079.264430.11984370.04481.619984356.604695.11985451.74579.719994576.104893.71986502.91635.320005137.6654441987609.60755.820015640.115923.61988749.09909.920026168.736436.51989850.711012.320037048.597570.21990934.651090.620048815.099469.919911045.731201.21. 对数据进行平稳化处理与检验 在SAS中，使用Gplot过程作出GDP序列的线图如图6.8左。 图6.8 GDP及其对数差分序列线图 线图显示GDP序列有明显的增长趋势，说明序列不平稳。对序列取对数并做一次差分后，再做线图，线图显示序列的增长趋势基本消除（如图6.8右），初步判断为平稳序列。图6.9 对数差分后GDP序列的纯随机性与平稳性检验 进一步使用ARIMA过程的indentify命令对对数差分后的GDP序列做纯随机性与平稳性检验，显示如图6.9所示。 用QLB统计量作的c2检验结果表明：对数差分后的GDP序列的QLB统计量的P值为0.0045（Chisq该QLB统计量的P值另自相关函数图中：Lag延迟阶数 Covariance延迟阶数给定后的自协方差函数 Correlation延迟阶数给定后的自相关函数 Std Error自相关函数的标准差“.”2倍标注差范围单位根检验结果详解见王燕P234，主要检查哪个p值小于显著水平a：类型1 无常数均值、无趋势的p阶自回归过程： 类型2 有常数均值、无趋势的p阶自回归过程： 类型3 有常数均值、又有趋势的p阶自回归过程： 2. 模型定阶、参数估计与检验图6.10 相对最优模型定阶 采用相对最优模型识别，根据上述分析及序列的自相关和偏自相关图，适当选择m = 4，n = 2，使用indentify命令中的minic p = (0: n) q = (0: m)短语进行相对最优模型定阶。结果显示（图6.10），在p = 1，q = 4时，BIC函数值最小。执行ARIMA过程的Estimate p = 1 q = 4命令做参数检验，结果未能通过参数检验。让q在03之间取值，通过反复测试，只有ARMA(1, 3)模型与ARMA(1, 0)模型通过参数检验及模型检验，其检验结果及参数估计如图6.11所示。 图6.11 对数差分后的GDP序列的ARMA(1, 3)模型与ARMA(1, 0)模型参数检验与模型检验结果 从上面2个模型的检验结果可以看到，它们均为有效模型，但ARMA(1, 0)模型的AIC为-67，SBC为-65均比ARMA(1, 3)的AIC与SBC小，根据AIC准则和SBC准则，前者应该更有效，所以应选择前者作为预测模型。3. GDP对数序列拟合ARIMA（p, d, q）模型的具体形式图6.12 拟合ARIMA（p, d, q）模型 对对数差分后的GDP序列执行ARIMA过程的Estimate p = 1 q = 0命令，得到GDP对数序列的ARIMA(1, 1, 0)模型如图6.12所示。 GDP对数序列模型的口径为：其中，xt表示GDP序列，模型可写为：4. 预测与分析 使用ARIMA过程的forecast命令输出GDP对数序列的预测结果，经过指数变换，得到2005年GDP的预测值，预测值与实际统计数据3的对比结果如表6.3所示。表6.3 河南省GDP预测值与实际统计数据对比结果预测值ARIMA(1, 0)实际值相对误差10656.4110535.21.15% 可以看出，ARIMA(1, 1, 0)模型对GDP预测的相对误差为1.15%，这说明预测值与实际结果是比较接近的。 进一步对20062008年GDP预测见表6.4。表6.4 河南省未来3年GDP预测值年份GDP年增长率200612666.2220.23%200714928.6117.86%200817521.2717.37% 结果表明，未来3年我省GDP仍将保持高速增长的态势。河南省委七届十次全会和省委经济工作会议提出要实现“十一五”期间生产总值年均增长10%的目标，是比较客观，容易实现的。 可以逐年使用当年新数据添加到数据序列中，利用上述方法建模对未来若干年的GDP进行短期预测，以对我省的发展做出更加科学的规划。十一. proc autoreg过程自回归过程autoreg用于估计和预测误差项自相关或异方差的时间序列数据的线性回归模型。自回归误差模型被用来校正自相关系数和广义自回归条件异方差模型GARCH（generalized autoregressive conditional heteroskedastic），并且其变体如广义的ARCH（GARCH）、方差无穷的GARCH（IGARCH）、指数的GARCH（EGARCH）和依均值的GARCH（GARCH-M）被用于异方差的建模和校正。自回归过程autoreg可以拟合任意阶的自回归误差模型，并且可以拟合子集自回归模型。为了诊断自相关性，过程产生广义Durbin-Watson（DW）统计量和其边缘概率。普通回归分析假定误差方差对于所有观察是相同的，但当误差方差不相同时，数据被称为异方差，此时普通最小二乘法估计不是有效的，同时也影响预测值置信区间的精确性。Autoreg过程能检验异方差，并且提供GARCH模型族来估计和校正数据易变性。对于带有自相关扰动和随时间变化的条件异方差模型，过程输出条件均值和条件方差的预测值。proc autoreg过程由下列语句控制：proc autoreg data=数据集 ;model 因变量=独立回归变量列表 ;output out=数据集 ;by 变量列表;run ; 其中，至少要有一个model语句。每个model语句都可跟随一个output语句。18. proc autoreg语句。l outest数据集名把估计参数输出到指定数据集中。l covout 数据集名把估计参数的协方差阵输出到指定数据集中。19. model语句的。l center通过减去均值中心化因变量并且取消模型的均值参数。l noint取消模型的均值参数。l nlag数值/（数值列表）指定自回归误差的阶或者自回归误差的时间间隔的子集。例如，nlag3与nlag（1 2 3）作用相同，但与nlag（1 3）等不同。l garch（）指定广义条件异方差GARCH模型的类型。例如，定义GARCH（2，1）回归模型时，可用下面SAS语句：l model y=x1 x2 /garch=（q=2,p=1）;l 请特别注意SAS系统的自回归参数符号q和p与我们前面所述公式中的符号p和q正好相反。定义GARCH-M（1，1）回归模型时，可用下面SAS语句：l model y=x1 x2 /garch=（q=2,p=1,mean）;l type=选择值，指定GARCH模型的类型：选择值为noineq时指定无约束GARCH模型，缺省值；选择值为nonneg时指定非负约束GARCH模型；选择值为stn时指定约束GARCH模型系数的和小于1；选择值为integ时指定IGARCH模型；选择值为exp时指定EGARCH模型。选项noint取消条件异方差模型中的均值参数。选项tr对GARCH模型的估计使用信赖区域方法，缺省值为对偶拟牛顿法l all要求打印所有输出选项。l archtest要求用portmantea Q检验统计量和Engle的拉格朗日乘子LM（Lagrange multiplier test）检验是否存在条件异方差情况，即是否有ARCH效应。l coef打印前几条观察的变换系数。l corrb打印参数估计的估计相关系数。l covb打印参数估计的估计协方差。l dwn打印直到阶n的DW统计量，缺省值n为1。l dwprob打印DW统计量的p值。当误差自由度大于300时dwprob选项被忽略。l ginv打印Yule-Walker解的自协方差的Toeplitz矩阵的逆。l itprint打印每步迭代的目标函数和参数估计。l lagdetp打印DW t统计量，它用于检验存在