时间序列模型

1. 时间序列模型

众所周知，时间序列的预测是需要假定在时间序列平稳随机过程的基础上，若为非平稳时间序列，容易造成对随机过程的伪回归，得到的回归函数及检验没有可信度， 那么到底什么是平稳时间序列 ？ 什么是随机过程（白噪声）？ 
  
  时间序列平稳性分为弱平稳和强平稳，弱平稳是序列的均值、标准差、协方差不随时间的变化而发生改变，强平稳则是序列的联合分布函数与时间的位移无关。 
  
 我们预测时间序列一般都用序列本身存在的前后依附关系，这种关系是时间序列预测的基础，我们这里假设本期数值只与前一期相关（与前t期相关的类似），则 平稳性的推导 过程如下
                                          
 所以当r小于1的时候，序列才为平稳序列，也称为随机游走过程，但这个至少均值还为零，是带漂移时间序列的特殊形式，那什么是带漂移序列，带漂移序列就是均值也随着时间t发生改变，即
                                          
 但是也由上图第二个公式可知，一阶差分之后的（上面三角形加Yt）的均值和方差保持不变(因为Ut服从正态分布)
  
 既然我们知道了什么是非平稳时间序列，那我们应该怎么检验呢？
  
 答案是单位根检验
  
  单位根检验由来 
                                          
 ADF检验是单位根检验的一种，也是比较常用的一种检验方法，是DF检验的优化，消除了序列的自相关性，内部使用为t统计量检验过程，若t大于临界值，则不能拒绝原假设r=1，序列非平稳，否则不能拒绝备泽假设，序列为平稳序列。
  
  白噪声 
  
 纯随机过程，是由一个不相关的随机变量的序列构成的，即不存在自相关，协方差=0。对于一个随机过程来讲，如果期望和方差均为常数，则称该过程为白噪声过程。之所以成为白噪声过程是因为它和白光的过程有些相似，白光的光谱在各个频率上有相同的强度，在各个频率上面的值相同。所以白噪声序列一定是平稳的。
  
  问题解答： 
  
 &白噪声和平稳性有什么区别？
  
 这么说吧，白噪声一定是平稳序列，因为方差和均值都不随时间的变化而变化，且不存在自相关，平稳性呢是如果不平稳的话就进行差分，差分的时候是对xt=x（t-1）+u的序列方差和均值，但实际如果一介差分是平稳的话，我们实际用的数据是（xt-x（t-1））的差值这个序列，所以还要对这个序列进行白噪声检验，但也有人说白噪声的意义不大，还有我看一本书上面说，白噪声是检验数据有用的信息有没有被提取完毕，如果是白噪声，说明信息提取完毕，剩下的全是随机扰动，无法预测和使用；另一个版本的解释是，如果通不过白噪声检验就要对其进行自相关和偏自相关模型识别，即ARMA中的q.p，说明白噪声比平稳序列多一个不存在自相关，这个问题我也没找到官方的解答资料  （所以不知道理解的对不对，欢迎大家指正）  
  
 &时间序列本身是与时间相关的，为什么又要求均值和方差与t不相关呢？
  
 可以这么理解，时间序列中的数据点的位置依赖于时间，即数据的取值依赖于时间的变化，但不一定是时间t的严格函数，而且数据点和t相关不代表方差与t相关，因为方差代表了数据在其均值上的离散程度 
   
                                          
 首先，我们看下上面这个回归模型，如果上面的模型合适，则u应该是平稳序列，则y对长期均衡关系的偏离是暂时的，即回归函数为长期均衡关系，相反，如果u为非平稳，则偏离会长时期不会消失，所以上面函数是否有价值，主要看u是否平稳
  
 也许上面的几个解释变量都是非平稳序列，但是他们的线性组合也许是平稳的，这就是协整的理论思想
  
 但是利用协整必须满足几个条件：（直接上图吧）
                                                                                  
  协整检验 
  
 利用残差平稳性检验，第一种是使用上面用到的ADF检验，但需要注意的是，协整检验的临界值不但与漂移项、趋势项有关，还与非平稳变量个数有关。我们一般用第二种DW检验，用协整回归得到的残差构造
                                          
 实际上，多变量时间序列的协整就相当于是做回归，但是要看时序在回归上的可信度，对残差进行检验，判断变量间是否为长期的均衡关系，并且这个地方要注意，使用协整的条件，必须是几个变量为同阶差分。
  
 最后说一个误差修正，因为长期关系的稳定关系一般都建立在短期动态的不断调整下得到并维持的，但由于变量的长期变量相互抵消，所以我们建立的模型看起来依然可信度高。那这个地方我们就要用误差修正模型来调节短期行为
  
  误差修正过程 ：
  
 首先建立长期关系模型，并使其得到平稳的残差序列，即模型合理。然后建立短期动态关系，即误差修正方程，将长期关系中的各变量差分之后重新构造，并将长期模型中的残差序列作为解释变量引入，对短期动态关系检验，逐项剔除，直到找到适合的方法为止
  
 理论案例：（直接上图）
  
 
  
                                          
  时间序列模型的建立过程： 
  
 首先，画出散点图观察并进行检验，检验序列是否是平稳序列，不平稳进行差分或者log变换，平稳则进行白噪声检验，没有通过白噪声的情况下就要进行模型识别，AR、MA和ARMA，确定后对模型的随机扰动项u进行检验，是否为白噪声序列，如果不是，则返回到前面，对模型重新识别。

2. 如何用神经网络进行时间序列预测

神经网络是可以用来预测时间序列。例如神经网络人口预测。已知1990至2009年的某地区人口数[11 28 30 42 44 56 49 60 50 63 56 74 76 65 92 105 124 117 132 128]。预测2010-2016年的某地区人口数。
具体实施过程：
%已知数据
t=1990:2009;
x=[11 28 30 42 44 56 49 60 50 63 56 74 76 65 92 105 124 117 132 128];
 % 自回归阶数
lag=3; 
%预测步数为fn
fn=length(t);
%输出数据
[f_out,iinput]=BP(x,lag,fn);   %BP（）神经网络预测函数 
[x' iinput']
R2=corrcoef(x,iinput)

%预测年份或某一时间段
%t1=2015:2016;
t1=length(x)+1:length(x)+7;
%预测步数为fn
fn=length(t1);     
[f_out,iinput]=BP(x,lag,fn);
P=vpa(f_out,5);
%预测数据
t1=2010:2016;
[t1' P']
% 画出预测图
figure(6),plot(t,x,'b*-'),hold on
plot(t(end):t1(end),[iinput(end),f_out],'rp-'),grid on
title('BP神经网络预测某地区人口数')
xlabel('年份'),ylabel('人口数');
legend('2009-2014年人口变化数','2014-2016年人口预测数');

3. 时间序列信号模型

图1-7 离散平稳随机信号x（n）的有理传输函数模型

随机序列主要采用自相关函数和功率谱密度函数进行研究。对于平稳随机序列，这些年来从时间序列分析角度，又提出另外一种研究方法，即时间序列信号模型法。这种模型是一个线性模型，它具有连续功率谱的特性，在功率谱估计方面，表现出很大的优点，对于研究平稳随机序列是一种很有效的方法。许多平稳随机序列都可以看成是由典型噪声源激励一个线性系统产生的，这种噪声源一般是白噪声序列源。假设该线性稳定系统的系统函数用H（z）表示，实际应用中所遇到的随机过程大多数可以用有理传输函数（系统函数）模型很好地逼近，如图1-7所示，图中，输入激励ε（n）是均值为零、方差为  的白噪声序列，线性系统传输函数为

地球物理信息处理基础

式中：bk是前馈（或滑动平均）支路的系数，称为MA（Moving Average）系数；ak是反馈（或自回归）支路的系数，称为AR（Autoregressive）系数。系统的输出序列x（n）是被建模的离散随机信号。该模型的输出和输入之间满足差分方程

地球物理信息处理基础

由式（1-106）知，输出功率谱和输入功率谱之间存在下列关系：

地球物理信息处理基础

这是因为

地球物理信息处理基础

或

地球物理信息处理基础

若h（n）是实的，则H*（1/z*）=H（z-1），于是

地球物理信息处理基础

以下讨论的都是这种情况。
由于｜H（ejω）｜的增益系数可并入  进行考虑，不失一般性，可假设a0=1和b0=1。

4. 时间序列模型的种类

 ARMA模型的全称是自回归移动平均(auto regression moving average)模型，它是目前最常用的拟合平稳序列的模型，它又可细分为AR模型(auto regression model)、MA模型(moving average model)和ARMA模型(auto regression moving average model)三大类。具有如下结构的模型称为 阶自回归模型，简记为：如果一个系统在某时刻的响应与其以前的响应无关，而与其以前进入系统的扰动存在一定的相关关系，这一类系统则称之为移动平均MA系统。这是因为是由一系列的及其滞后项的加权和构造而成。这里的“移动”指的变化，而“平均”指加权和。一般移动平均模型由部分构成，形成如下：为了分析的方便将其表述为与系统因素的延迟项一致，即将模型中各加号改为减号有：用滞后因子表示为：把具有如下结构的模型称为自回归移动平均模型，简记为：引进延迟算子，模型简记为：式中：，为阶自回归系数多项式。，为阶移动平均系数多项式。限制条件条件一：这个限制条件保证了模型的最高阶数。条件二：这个限制条件实际上是要求随机干扰序列 为零均值白噪声序列。条件三：这个限制条件说明当期的随机干扰与过去的序列值无关。 ARIMA模型又称自回归求和移动平均模型，当时间序列本身不是平稳的时候，如果它的增量，即的一次差分，稳定在零点附近，可以将看成是平稳序列。在实际的问题中，所遇到的多数非平稳序列可以通过一次或多次差分后成为平稳时间序列，则可以建立模型：这说明任何非平稳序列只要通过适当阶数的差分运算实现差分后平稳，就可以对差分后序列进行ARIMA模型拟合了。模型是指阶差分后自相关最高阶数为，移动平均最高阶数为的模型，通常它包含个独立的未知系数：。它可以用最小均方误差原则实现预测：用历史观察值的线性函数表示为：式中，的值由下列等式确定：如果把记为广义自相关函数，有容易验证的值满足如下递推公式：那么，的真实值为：由于的不可获取性，所以的估计值只能为：真实值与预测值之间的均方误差为：要使均方误差最小，当且仅当，所以在均方误差最小原则下，期预报值为：预测误差为：真实值等于预测值加上预测误差：其中，预测误差的均值和方差分别为：

5. 三种时间序列模型

（1）如果除a0=1外所有其它的AR系数都等于零，则式（1-124）成为

地球物理信息处理基础

这种模型称为q阶滑动平均模型或简称为MA（q）模型（Moving Average Model），其系统函数（传输函数）为

地球物理信息处理基础

模型输出功率谱为

地球物理信息处理基础

或

地球物理信息处理基础

这是一个全零点模型，因为它只有零点，没有极点（除了原点以外）。如果模型的全部零点都在单位圆内，则是一个最小相位系统，且模型是可逆的。
（2）如果除b0=1外所有其它的MA系数都等于零，则式（1-124）成为

地球物理信息处理基础

这种模型称为p阶自回归模型或简称为AR（p）模型（Autoregressive Model），其传输函数为

地球物理信息处理基础

模型输出功率谱为

地球物理信息处理基础

或

地球物理信息处理基础

显然，该模型只有极点，没有零点（除了原点以外），因此这是一个全极点模型，而且只有当极点都在单位圆内时，模型才稳定。
（3）设a0=1和b0=1，其余所有的ak和bk不全为零。在这种情况下，模型的差分方程、系统函数和输出功率谱分别用式（1-124）、式（1-123）和式（1-125）或式（1-126）表示。分子部分称为MA部分，而分母部分称为AR部分，这两部分分别满足稳定性和可逆性的条件。这是一个“极点—零点”模型，称为自回归滑动平均模型ARMA（p，q）模型（Autore-gressive Moving Average Model）。
在上面已谈到，实际中所遇到的功率谱可分为三种：一种是“平谱”，即白噪声谱，第二种是“线谱”，即由一个或多个正弦信号所组成的信号的功率谱，第三种介于二者之间，即既有峰点又有谷点的谱，这种谱称为ARMA谱。可以看出，AR模型能突出反映谱的峰值，而MA模型能突出反映谱的谷值。
沃尔德（Wold）分解定理阐明了上述三类模型之间的联系，即：任何广义平稳随机过程都可分解成一个可预测（确定）的部分和一个不可预测（完全随机）的部分。确定性随机过程是一个可以根据其过去的无限个取样值完全加以预测的随机过程。例如，一个由纯正弦信号（具有随机相位以保证广义平稳）和白噪声组成的随机过程，可以分解成一个纯随机成分（白噪声）和一个确定性成分（正弦信号）。或者可以把这种分解看成为把功率谱分解成一个表示白噪声的连续成分和一个表示正弦信号的离散成分（具有冲激信号的形式）。
Wold分解定理的一个推论是：如果功率谱完全是连续的，那么任何ARMA过程（Au-toregressive Moving Average Process）或AR过程（Autoregressive Process）可以用一个无限阶的MA过程（Moving Average Process）表示。Колмогоров（Kolmogorov）提出的一个具有类似结论的定理：任何ARMA或MA过程可以用一个无限阶的AR过程表示。这些定理很重要，因为如果选择了一个不合适的模型，但只要模型的阶数足够高，它仍然能够比较好地逼近被建模的随机过程。
估计ARMA或MA模型参数一般需要解一组非线性方程，而估计AR模型参数通常只需解一组线性方程，因此，AR模型得到了深入的研究和广泛应用。如果被估计过程是p阶自回归过程，那么用AR（p）模型即能很精确地模拟它；如果被估计过程是ARMA或MA过程，或者是高于p阶的AR过程，那么用AR（p）模型作为它们的模型时，虽然不可能很精确，但却可以尽可能地逼近它，关键是要选择足够高的阶数。证明如下：
假设MA模型为

地球物理信息处理基础

对上式进行Z变换得到
X（z）=B（z）W（z）
式中B（z）是MA信号模型的系统函数，或者说是bi（i=1，2，3，…）序列的Z变换。
设MA信号模型满足可逆性条件，即B-1（z）的存在，令
B-1（z）=G（z）=1+g1z-1+g2z-2+g3z-3+…
这样
X（z）G（z）=（1+g1z-1+g2z-2+g3z-3+…）X（z）=W（z）
则

地球物理信息处理基础

对上式进行Z反变换，得到
x（n）+g1x（n-1）+g2x（n-2）+g3x（n-3）+…=w（n）
上式就是x（n）的AR信号模型，因此证明了一个时间序列可以用有限阶MA信号模型表示时，也可以用无限阶的AR模型表示，对于ARMA模型也同样可以证明。
［例1-2］已知x（n）的功率谱为

地球物理信息处理基础

求出该模型的系统函数H（z）。
解：利用欧拉公式可以将Pxx（ejω）变为

地球物理信息处理基础

取z=ejω，则上式变为

地球物理信息处理基础

令  ，那么，  ，显然有理多项式B（z）的分子、分母都是最小相位的。所以有

地球物理信息处理基础

与式（1-120）相比较，得  。又由式（1-125）得到所求的系统函数

地球物理信息处理基础

6. 时间序列的分析模型

时间数列的组合模型1 加法模型：Y=T+S+C+I (Y,T 计量单位相同的总量指标)(S,C,I 对长期趋势产生的或正或负的偏差)2 乘法模型：Y=T·S·C·I(常用模型) (Y,T 计量单位相同的总量指标)(S,C,I 对原数列指标增加或减少的百分比)

7. 时间序列模型的用途

 根据时间序列模型可调整输入变量使系统发展过程保持在目标值上，即预测到过程要偏离目标时便可进行必要的控制。

8. 时间序列模型简介

  目录 
   时间序列是一列观测值  的集合, 其中每个观测值是在时段  观测所得(  是自然数 ). 给定时间序列  , 如果对任意的  , 它满足下列条件:   i.      ii.      iii.      我们把它叫做(弱)平稳(weakly stationary)序列.(下文我们简称平稳序列.)
   通俗地讲,  平稳序列的期望, 方差, 协方差不随时间变化 . 例如,   服从同一个分布时, 它是平稳的.
    例1  下图中的时间序列由  生成. 从直观上看, 这个序列是"平稳的".   
                                           
    例2  下图的中的时间序列由  生成, 其中  ,   . 它起初有明显地增长, 然后趋于平稳. 利用ADF检验(详情见下文), 我们发现该序列是平稳的(p-value < 0.01).   
                                           
    Remark  弱平稳性的"弱"主要体现在时间序列在全局上是平稳的, 即，时间序列局部是波动的，但整体上看是平稳的, 或者随着时间的变化其样本的均值收敛.
   我们用统计学中假设检验的方法来判断样本的平稳性. 常用的是Augmented Dickey-Fuller(ADF)检验  [1]  .
   在显著水平  的条件下, 我们可以通过计算p-value来接受或者拒绝  :
   Python3中 statsmodels.tas.stattools 中的 adfuller 函数  [3]  实现了ADF检验. 使用方法如下所示.
   前面之所以介绍平稳序列的概念及检验方法, 是因为它是很多基础的时间序列模型的前提假设. 在本节我们介绍一些常见的时间序列模型(更多内容可以参考  [4]  ,   [5]  ).
   AR代表自回归(Autoregression). 假设时间序列  是平稳的, 它可以被表示成如下形式:     
   MA代表移动平均(Moving Average). 假设时间序列  是平稳的, 它可以被表示成如下形式:     
   ARMA模型是AR和MA的组合. 假设同上. 它可以被表示为如下形式:     
   ARIMA模型是ARMA模型的推广, 全称是Autoregressive Integrated Moving Average. 当时间序列  不满足平稳性时, 我们通常使用 差分 的技巧把序列变得平稳, 然后再应用ARMA模型.
   参数  代表差分的阶数. 下面是差分的计算公式(  为差分算子):
    例3  下图是原始的时间序列. 通过观察, 它的均值有明显的上升趋势且不收敛, 因此不是平稳序列(ADF检验的p-value为0.94).
                                           对该序列进行一阶差分后, 我们得到如下平稳的时间序列(p-value为0.00).
                                           该记号代表季节性(或周期性)ARIMA模型, 详细的表达式可以参考  [4]  ( 4.1 Seasonal ARIMA models ), 其中
   我们可以把它看成两阶段模型: 第一阶段在全局使用ARIMA(p,d,q); 第二阶段通过指定周期长度  , 再利用ARIMA(P,Q,D)模型考虑周期之间的关系.
    例4  考虑如下周期性的平稳时间序列(  ).
                                           对序列进行周期性差分:   得到新的时间序列  如下图所示(红色部分)
                                           通过使用周期性差分, 我们可以把原有时间序列的周期性移除. 同理, 通过采用周期性的自回归和移动平均系数, 我们可以把周期之间的依赖关系考虑进模型.
    例5  考虑周期s=18的数据(蓝色曲线). 用  和  分别进行预测的结果如下.
                                           不考虑周期性的ARIMA模型的预测结果(灰色曲线)逐渐收敛到时间序列的均值. 由于序列是平稳的, 这样的预测结果符合我们的期望. 考虑到该时间序列有比较强的周期性, 且通过观察发现周期  . 在本例中, 我们仅使用周期差分,  最终得到了如图所示(红色曲线)的周期性预测结果.
   ARCH的全称是Autoregressive Conditionally Heteroscedasticity, 它可以用来考虑样本的方差随着时间变化(或震荡)的时间序列. 设时间序列  是平稳的,   模型可以被表示成如下形式:
        其中     
   GARCH即Generalized ARCH, 是ARCH模型的推广  [6]  . 设时间序列  是平稳的,   模型可以被表示成如下形式:        其中     
    Remark  ARCH/GARCH随机过程产生的数据是什么样的? 前面提到它们允许 样本的方差 随时间变化, 但是由于  必须满足平稳性(前提假设), 因此样本的方差从局部看是变化(震荡)的, 但从整体看应该是"平稳的"序列. 例如下图是一个  过程生成的时间序列(  ).
                                           VAR即Vector Autoregression, 它是多变量的自回归模型. 类似地, 我们有  , 它是  的向量版本. 需要注意的是, VARMA模型处理的时间序列可以有趋势. 我们不做详细的展开, 感兴趣的读者可以参考  [4]   章节11.2: Vector Autoregressive models VAR(p) models .
   给定时间序列的观测样本, 选定预测模型之后如何确定模型的参数? 本节我们介绍两种常用的方法: 1. 画出ACF/PACF图, 然后观察出  的值; 2. 通过计算相关的统计指标, 自动化地选择参数.
   ACF的全称是Autocorrelation Function. 对变量  , ACF的值代表  与  之间的相关性.     
   PACF的全称是Partial Autocorrelation Function. 对变量  , PACF的值代表已知  的条件 下,   与  之间的相关性.     
    例6  设  . 考虑下面三个模型生成的时间序列, 并计算相应的ACF/PACF.
                                                                                                                           基本思想是通过计算一些指标, 并选择参数使得相关的指标值尽可能小. 下面我们介绍一些常用的指标.
   为方便描述, 我们先定义一些记号.
     
   (AIC的改良版, 解决小样本过拟合的问题)     
   (也称为Schwartz Criterion, SBC, SBIC)
     
     
    Remark  建议在实际中综合考虑这些指标.
    Python3 code on Github