数学建模提高班专题5——时间序列建模

1、柴中林柴中林 2015/5/92015/5/9课 件 提 纲v1 差分方程的引例与概念v2 特殊差分方程的解v3 平衡点及其稳定性v4 差分方程组v5时间序列与其中的趋势分析v6自回归模型v7自回归模型识别及参数确定v8自回归模型预测及相关说明v9 建模练习题1 差分方程的引例与概念差分方程的引例与概念 例1. 某人贷款80万元买了一套房子，期限20年.已知贷款月利率为5 ，请问他每月要还贷多少？ 在高数中，我们研究的函数中的变量的取值大都是连续的（在连续区间上取值，如（-,+ ）, (-1,12）。但在经济管理和很多实际问题中，变量只能取1,2,3,这样的值。这些变量称为离散型变量。描述离散

2、型变量间关系的模型称为离散型模型。差分方程就是常见的一种离散型模型。 微分方程：连续变量间存在函数关系。知道了这个关系，就能够研究变量间的联系与变化规律。然而，这个关系是未知的，但我们可以建立起含自变量，因变量及其导数或微分的等式，这就是微分方程。通过对方程的研究以求得这个函数关系，或通过方程直接揭示变量间的联系就构成了微分方程的主要研究内容。 差分方程与微分方程是类似的。只是这里的变量是离散的。差分方程：含自变量，未知函数（因变量），未知函数差分的等式。建立差分方程，求解它的目的就是研究离散变量间的关系。 一般的，对有函数关系的两个变量，常用x当自变量，y当因变量。但在差分方程中，因自变量只

3、取整数值（如0,1,2，），我们更喜欢用n(或t)表示自变量，这时因变量可用x或y表示。其函数关系是x=x(n)，但我们更常用xn表示。当然，这个关系是不知道的，但我们常能得到的是如下的式子 F(n, xn, xn-1, xn-k)=0 （1） 或 G(n, xn, xn+1, xn+k)=0 （2） 或 H(n, xn, xn, kxn)=0 （3） 这种式子就是差分方程。 有时，（1）也写成如下的形式 xn=f (n, xn-1, xn-k) (4)因此，差分方程也称为递推关系。 考虑例1，用n表示月份(n=0表示贷款月份)，xn表示第n月还贷后还欠的钱，r, a分别表示银行月利率和月还钱

4、数， xn 表示了账户中欠钱数与月份间的函数关系（ 未知），但我们容易得到一个式子 xn +rxn -a = xn+1即 xn+1-(1+r) xn =-a (5) 此外，还有初始条件 x0=80（万元）及x240=0。 这就是贷款问题的差分方程模型。 变化建模比较微元法 对离散关系xn，其函数值构成序（数）列xn(x1,x2,x3,)。记 xn=xn+1-xn （序列后项减前项构成的序列）, 称为xn的一阶差分，2 xn = ( xn)= xn +1-xn = xn +2-2 xn +1+xn称为xn的二阶差分，依次类推。 对式 (5)：xn+1-(1+r) xn =-a ，也可将它写为：

5、xn-(1+r) xn-1=-a 或xn-rxn=-a （差分方程因此而得名）. 即同一个关系用不同视角不同符号式子会不同，但可以互化（它们是同一个东西）。 差分方程中的最高阶差分的阶或因变量的最大下标与最小下标之差称为差分方程的阶。 差分方程的解是函数，通常有无穷多个。通解是全部解的集合（体现在独立任意常数上，其个数与方程阶数相同）。另外，在实际问题中常会给出一些附加条件（如x0的值），称为初始条件。满足初始条件的具体的解就是特解。 差分方程问题的研究内容： 1 差分方程的建立（离散变量关系的建立，也可将连续问题离散化）； 2 差分方程的求解和分析。 差分方程在实际问题中有广泛的应用。 差分

6、方程的求解并不比微分方程容易，大部分差分方程是无法求解的。这里介绍最简单同时用处很大的一类特殊差分方程的求解。 常系数线性齐次差分方程，其一般形式为 xn+a1xn-1+akxn-k=0 (6) 其中a1,ak是常数。 方程（6）有解，其求解步骤为:步骤1: 求解对应的特征方程 k +a1 k-1+ak=0 （7） 步骤2: 根据步骤1的解的情况写出（6）的通解； 2 特殊差分方程的解特殊差分方程的解 情况1：若是（7）的一个单实根，则n是（6）的一个特解。若1, 2, k是（7）的k个全部不同的单实根，则（6）的通解为xn=C1 1n +C2 2n+Ck kn（ C1 ，C2 ，Ck 是任意

7、常数）。 情况2：若是（7）的k重实根，则n, nn, , nk-1n都是（6）的特解。 情况3：若=i是（7）的单重复根，则 ncos n与n sin n都是（6）的特解，其中 ,是的模与幅角主值。情况4 ：若=i是（7）的k重复根，则 n cos n, nn cos n, nk-1n cos n与n sin n, nn sin n, nk-1n sin n都是（6）的特解，其中 ,是的模与幅角。 最后，将各个特解如情况1那样与任意常数相组合就得（6）的通解。 常系数线性非齐次差分方程，其一般形式为 xn+a1xn-1+akxn-k= b(n)(0) (8) (8)的求解方法是先求相应齐次方

8、程的通解，记为xn*，再求（8）的一个特解，记为 xn (0) （方法：根据b(n) 的特点将xn (0)的形式设出，再用待定系数法确定其中的系数），于是（8）的通解为 xn = xn* + xn (0) 此外，不同于微分方程，对差分方程，当初始条件给定后，可迭代求得任意xn的（精确）值，从而可以对xn的变化规律进行作图分析。如对方程xn=f(n, xn-1, xn-k)，若x1, x2, xk 给定，就可以根据方程依次算出xk+1, xk+2, xk+3 来。 下面求解例1： xn+1-(1+r) xn =-a 。它是一阶常系数非齐次线性差分方程。先解相应的齐次方程xn+1=(1+r) xn

9、，特征方程为= 1+r，其通解为xn *=C (1+r)n （ C为任意常数）, 再求其一个特解。从方程看设xn为常数（记为x），代入得xn (0) =a/r, 于是得方程通解：xn =C (1+r)n +a/r。 代入初始条件得方程组 rarCraC/)1 (0/80240) 1)1/()1 (80240240rrra解之得大约是5731元.3 平衡点及其稳定性平衡点及其稳定性 差分方程虽可用迭代法进行数值计算，但计算总归只能进行有限步，其深层次的性质必须用其它工具进行分析，平衡点就是其中一个。 平衡点相当于稳定点或不动点，对方程xn=f(n, xn-1, xn-k) 来说就是若xn-1,

10、xn-k都取某一常数，比如a，那么xn也一定是a，从而xn+1, xn+2, xn+3, 也都将取a 。 平衡点就是所有xn都取相同的值，且能使方程成立的点，于是将xn=f(n, xn-1, xn-k) 中所有xn都换成x，得方程x=f(n, x, x) ，将其求解，每一个解就是一个平衡点。 设a是方程的一个平衡点， xn是方程的任一解，若总有则称a是差分方程的一个稳定的平衡点（为什么？）。 稳定的平衡点在实际问题中有重要的价值。 limnnxa 现考虑差分方程 xn+a1xn-1+akxn-k=0 ，并且其解是如下形式 xn=C1 1n +C2 2n+Ck kn 。 显然0是方程的一个平衡点

11、，不难发现对任意s若有|s|1，则必有 这说明0是稳定的平衡点，这也是一般差分方程平衡点稳定性的判别方法：若齐次方程的特征方程的根的绝对值都小于1，平衡点稳定。而若某个的绝对值大于1，平衡点不稳。当等于1时，有多种情况且实际意义不大，不做讨论。若特征根是复根，就用其模来判断。 lim0nnx 例2 考虑数学模型书中供需关系的蛛网模型： xk：第第k时段商时段商品数量；品数量；yk：第第k时段商品价格时段商品价格，需求函数，需求函数 yk=f(xk)，供需平衡点，供需平衡点为为P0(x0,y0)。当商品生产者的生产只盯着前一期价格（供应函数。当商品生产者的生产只盯着前一期价格（供应函数为为xk+

12、1=g(yk)）时，在平衡点附近各时段商品数量的差分方程模）时，在平衡点附近各时段商品数量的差分方程模型为型为xK+1+xk=(1+)x0. 其齐次方程的特征方程的特征根为其齐次方程的特征方程的特征根为- 。所以。所以| - |=1就稳定，否则就不稳。就稳定，否则就不稳。 而当商品生产者的生产同时盯着前面两期的价格（供应函数而当商品生产者的生产同时盯着前面两期的价格（供应函数为为xk+1=g(yk+yk-1)/2）时，在平衡点附近各时段商品数量的差分）时，在平衡点附近各时段商品数量的差分方程模型为方程模型为2xK+2+xk+1+xk=2(1+)x0. 其齐次方程的特征其齐次方程的特征方程为方程

13、为22 2+ + + +=0。特征根。特征根为为4/ )8)(22, 1 当 8时，根为实根，必有一根绝对值大于时，根为实根，必有一根绝对值大于1.当当08时根为复根。用复数的模来判断，可以得到当时根为复根。用复数的模来判断，可以得到当02时稳定，否则不稳。时稳定，否则不稳。 差分方程组（自变量一个，因变量多个，仅讨论线性） 线性差分方程组的一般形式为 其中aij和bi (i, j=1,2,n)都是常数。1111122112211222221122(1)( )( ).,( )(1)( )( ).( )(1)( )( ).,( )nnnnnnnnnnnx ta x ta xta xtbxta x

14、 ta xtaxtbxta x taxta xtb4 差分方程组差分方程组 令 记x(t)=(x1(t), x2(t), , xn(t)T, b=(b1, b2, , bn)T，则上述方程可记为x (t+1)=Ax(t)+b。 该式类似于前面的一阶常系数线性差分方程，可编程数值计算分析，也可利用线性代数理论（主要是特征值和特征向量）进行分析讨论。若x*（向量）是该方程的一个平衡点（ 即x*=A x* +b ），则它稳定的条件是A的所有特征值的绝对值都小于1，若某一个的绝对值大于1,就不稳。 111212122212nnnnnnaaaaaaAaaa5 时间序列与其中的趋势分析时间序列与其中的趋势

15、分析 时间序列：按时间（有时是长度或温度）顺序排列的随机变量序列，但在应用中又指将某个统计指标在不同时间上的各个数值，按时间先后顺序排列而形成的序列（一般等间隔）。 时间序列分析：根据观测得到的时间序列数据(其机理未知)，通过曲线拟合和参数估计来建立数学模型和理论，希望从中寻找出变量的变化规律，对未来的某些阶段进行预测。 时间序列有广泛的应用。 设yt 是时间序列，虽然它暗含了时间变量t，但它仅指采样的时间点。因此，一般的不能认为y是纯t的函数，从而按回归等其他理论去做。因为许多变量都随着时间的变化而变化，所以时间序列中也常常包含因时间变量而产生的趋势变化。另外，在时间序列中，相近的各项间往往

16、有很强的依赖关系：当前的数值对下面的数值有很强的影响（如股市，期货）。此外，每个数据还受到无法刻画捕捉的随机因素的影响。通常yt 可表为yt =f(t) +xt，其中f(t)表示随时间变化的确定性趋势，xt则主要由随机因素或其积累而形成，是一个平稳序列。 在yt =f(t) +xt中，趋势成分f(t)起着主导的作用。当它存在时， xt可以认为是随机误差，并予以忽略，故可以用回归方法确定f(t)中的参数，得到f(t) 。影响f(t)的因素 有长期趋势，季节变动（季节性规律作用产生的周期变化），循环变动（周期长短不固定的一种变化）以及不规则的变动等。通常，趋势成分主要讨论长期趋势和季节变动趋势，这

17、里也是。 当f(t)是由长期趋势决定的，其表达式可能是 线性趋势 f(t)=a+bt 二次曲线趋势 f(t)=b0+b1t+b2t2 或更高阶多项式趋势 幂函数曲线趋势 f(t)=atb 对数曲线趋势 f(t)=a+blnt 双曲线趋势 f(t)=a+b/t,或 1/f(t)=a+b/t 指数曲线趋势 f(t)=aebt修正指数曲线趋势 f(t)=L+aebt,或 f(t)=L+abt (a0,0b1) 龚泊兹曲线趋势 (0a0,0bq时全为0的性质称为q步截尾性。若它不能在某步之后截尾，而是随着k的增大而迅速衰减到0，受一负指数函数（如 y=e-kx ）控制，或如正弦函数似的震荡，称为拖尾性

18、。此外，由于随机性， k全为0是不可能的。因此，截尾是指k突然变的很小，并很接近于0. 注：AR(p)模型平稳的充要条件是它的p 个特征根都在单位圆内。MA(q)模型总是平稳的。 ARMA(p,q)的平稳性与其AR(p)部分相同。 当序列Xt非平稳时，说明趋势f(t)存在，除了季节趋势和明显的指数增长或阻滞增长趋势外，在短期内一般可用多项式函数近似。当为多项式函数时（如yt=a+bt+t），通过不断的差分就可得到一个平稳序列。因此，对序列进行差分（季节规律用季节差分），是将非平稳序列变为一个平稳序列的常用手段(一般不超过两次)。续例5，yt的图像明显上升，Dt具有明显周期性，都非平稳。它们的自

19、相关系数图像如下。由图像可知，非平稳图中有两条对称的蓝线，是随机变量的2线。落在线内说明可以接受相关系数为0（95%的置信度），线外则不可。7.3 白噪声检验 对于时间序列yt，需要把其中的规律或项间关联全部提炼出来，使得残差t（余下的部分）仅为一个白噪声。因此，检验残差序列是否为白噪声是判断模型是否合理以及建模是否需要终止的一个条件。设模型的残差序列为t。记21/)(tkntkttkr计算 其中n为数据个数，m为最大时滞（m视数据多少取n/4，n/10或n0.5）。Qm近似服从2(m) 分布。对给定显著性水平 ，若Qm大于2（m），则拒绝假设，否则接受（认为是白噪声）。此外，t是否为白噪声也

20、可通过其相关函数来判断。若其k和kk都很小，可认为是。否则，不是。 另：数模书中346页的投资问题给出了一个残差序列自相关性诊断方法。可画出t- t-1(r(1:n-1),r(2:n)的图像观察，也可用DW检验：（仅检验一阶相关性，但一般够了）。mkkmrnQ12)()/1 (22221nttntttDW 续例3. f(t)用三次多项式拟合。 22.372 .f(t)=142.2+102.7t-7.67t2+0.22t3 。残差序列为xt=yt-f(t)。其图像如右上。可以看出，序列基本平稳。但从t- t-1图像（右下）看，残差相邻项间有很强的相关性。用DW检验，算得的DW = 0.9757，

21、应通不过检验. 此时Q8 = 22.372，也不满足白噪声检验。这说明残差序列间尚有信息留待提取。7.4 自回归模型识别和参数确定自回归有AR，MA和ARMA三种模型，可从k 与k的特性来判别： AR(p)模型 k拖尾， kk滞后p阶后截尾。 MA(q)模型 k：滞后q阶后截尾， kk ：拖尾。 ARMA(p,q)模型 k ：拖尾， kk ：拖尾。 续例5，将Dt进行季节差分（zt=Dt+4-Dt），画图如下 由图像看，虽然zt没有明显的上升或下降趋势，以及季节特征，但并不很好的表现出围绕均值的波动。将zt再差分，仍记为 zt，其图像如下。可以看出，序列平稳（检验自己做）了。 下图是最后得到的

22、zt的自相关函数和偏相关函数图像。由图知，自相关函数可以认为是拖尾的，偏相关函数则是截尾的（在2或4处， zt接近于白噪声，大家去检验）。故应选AR模型。 下图是例3的残差xt的自相关函数和偏相关函数图像。由图知，自相关函数可以认为是拖尾的，偏相关函数则是截尾的（在2处）。故应选AR模型（显然不能认为是白噪声）。 k与kk的截尾处的（严格）判断。 k ：若在某个q0 （含）之前，k显著不为0.当q=q0， q0+1, q0+2, q0+M中满足式 的个数少于M的68.3%。或上面不等式右端乘2，但比例变为95.5%,则可近似认为k在q0处截尾，其中N为数据个数， M同上。对kk ，判断方法类似

23、，只是不等式是|kk|1/N0.5,或|kk|2/N0.5 。 截尾值q0可用来判别序列自回归或移动的阶数。Nqiik/)21 (|5 . 012 若k （要特别关注）与kk既不拖尾也不截尾，说明序列非平稳，或有季节特征，需进行相关处理。 虽然ARMA(p,q)模型具有一般性，但它也最复杂。另外，用ARMA(p,q)模型时，各t-k通常是未知的（不可观测量）。因此，当用该模型时，必须求出前面的各t，这不容易。 对AR(p)模型，其偏相关函数截尾的值（基本）就是回归的阶数。此外，也可用不同阶模式进行回归，残差平方和最小的值就是回归的阶数（看后面）。对MA(q)，方法是类似的。对ARMA (p，q

24、)模型，不能直接从相关函数得到大致的阶数。但残差平方和规则仍适用。方法是从低阶开始，向高阶拟合，在拟合的模型中选残差小者。或者遇到第一个残差可认为是白噪声的模式即停止。 注意：当用高阶自回归移动平均模型去拟合序列时，拟合的效果总会提高的（不可能降低），即残差平方和会下降。但到了一定阶数后，阶数的再提高产生的效果会是微小的，非实质性的。这时的拟合属于过拟合（拟合过度）。在建模中对模型还有一个“简约性”要求，即在精度相近的模型中我们要选择简单模型。为此，又有一个一般的定阶准则：AIC准则 。记 a2=模型的剩余平方和/(实际观察数据个数-模型中参数个数) 则 AIC(p,q)=loga2+2 (p

25、+q) /n。p,q的确定应使AIC(p,q)达到最小。该式子应该也适合AR和MA模型。但这个方法必须对多种模型求参数，拟合，算残差，计算量大，故仅对ARMA用它，且尽量避免用它。7.3 模型的定阶 确定了模型的阶数，就要确定其中的各个系数。一个常用的准则是残差平方最小准则：模型中各个系数的确定应使得用模型计算各个时刻的x值时，残差（实际值与计算值的差）的平方和达到最小。对AR(p)模型，根据回归原理，可将xt-1, xt-2, , xt-p作为自变量，xt作为函数，用命令regress去做。还可用如下式子计算回归系数。 对MA (q)模型，一阶的仍可用回归，高阶的就不行了，只好根据定义来。比

26、如对MA (2)，有xt=t -1t-1-2t-2，即t = xt +1t-1+2t-2。设给定序列x1, x2, x15，我们令 x1=1, 则2= x2 +11= x2 +1x1, 3 = x3 +12+21 = x3 +1(x2 +1x1) +2x1 ， 如此可得到各个t，1，2的确定应使得下式最小pppppp21212121111117.4 回归系数的计算niiL1221),( 此外，对MA(q)模型，高阶的可用下式确定其系数 这是非线性方程，求解也不容易。此外，专门软件如SAS中也许有命令计算，也许新版的matlab也有，大家查查看。对ARMA (p，q)模型，可用类似于MA (q）

27、的方法确定系数。当然，更复杂。也可用如下近似方法，先用回归方法求出自回归的系数，再用xt-1 xt-1-2xt-2-pxt-p=t -1t-1-2t-2-qt-q这个移动平均模型来做。 当 p，q很大时模型会很复杂，计算也困难。但实际中这种模型是不多的（一般小于等于3，甚或2）。除非能显著的减少误差，我们都尽量用简单模型来做。qkqkqkqqkkkk, 01 ,12222122118 模型预测及相关说明 得到一个模型，当然要用模型。而其中一个重要的应用就是预测。对AR模型 xt=1 xt-1+2xt-2+ + pxt-p+t，舍去当前误差t（均值为0的正态分布，未知），利用最近的p个数值通过递

28、推不断的得到下面的值，并不断的更新最近的值。当然，预测的时段越长，效果越差。 对MA模型xt=t -1t-1-2t-2-qt-q ，因为各时段误差t未知，故模型不能直接用。用的方法是从最初开始，将之前的误差设定为0，用每次的实际值与预测值的差作为前面各时刻的误差，即t=xt +1t-1+2t-2+qt-q 。 之后，就利用模型预测未来的xt 。只是我们预测时每次都把当前误差t设定为0，因此，当q步之后，所有的xt都应是0. 对ARMA模型 xt=1 xt-1+2xt-2+pxt-p+t -1t-1-2t-2-qt-q，方法与MA类似，须先算出前面各时段的误差，才能对未来预测。但q步之后，预测就

29、与自回归一样了。 对一个时间序列，当它非平稳时，需要对原始序列进行处理，这相当于做变量替换得到新序列，而我们最终的目的是要得到一个原始序列的预测模型，所以当用回归方法得到处理后的序列的模型后，就需要进行变量回代以得到这个模型。 一些相关说明 自回归模型适合于序列的前后数据间具有依赖关系。应用自回归模型解决时间序列问题，首先必须觉得时间序列有自回归的特点，从模型辩识和回归结果看较满意才可用，否则没什么意义。 此外，在自回归模型中认为外界对序列除了白噪声外没有其他影响，但现实显然不全如此。当序列除了自身的相互影响外，还受到其它可处理的显著因素的影响，就要把它们结合进模型中。如数模书中346页的投资

30、问题模型。这已太难并超范围了，略去。 对最后的序列zt，用回归方法得模型zt=0.6022zt-1-0.4417zt-2+ 0.1477zt-3+ 0.2479zt-4。再检验模型残差是否是白噪声。因为数据仅有15个，较少，故取m=4，利用回归残差算得Q4= 3.4930.而查表值20.25(4)=9.49.前者小于后者，可认为残差属白噪声，序列中的信息完全提取。利用得到的模型回代得到原始序列的模型及模型的预测略去。续例59 建模练习题建模练习题 1 一对夫妇8年前贷款40万元买了一套房子，20年期，每月还贷2810.62元，请问月利率是多少（会得一个非线性方程，用近似方法求零点），现在还欠多少钱？如果这对夫妇现在每月可多还款1500元，因此打算申请缩短还贷时间，若利率不变，他们最多能提前几年还清? 2 下表是我国1990年1月-1996年12月的工业生产总值(单位：千亿；M:月份)，请据此建立合适的模型。 时间值时间值时间值时间值时间值时间值90M11.421491M31.893992M52.373193M72.862994M93.663795M114.091090M21.367491M41.969892M62.515893M82.864294M103.753495M124.650890M31.7197