平稳时间序列以及MATLAB相关工具箱学习笔记

平稳时间序列以及MATLAB相关工具箱学习笔记

概念

（1）平稳序列

即序列的均值是个常数，与序列长度、起始位置无关。
直观看上去，该序列类似于围绕某一个值上下波动（该值为平稳序列均值）。

（2）平稳白噪声序列

所谓噪声，可以理解为一种非周期性的扰动。由于其自协方差函数值为0，于是其各序列值没有任何相关关系，所以说平稳白噪声序列是一段没有记忆的平稳序列。
在MATLAB中可用如下代码生成高斯分布下的平稳白噪声序列：

elps=randn(1,1000);
plot(elps)

基本的平稳时间序列

（1）$AR(p)$序列

需要说明的是：
1.$p$为模型阶数，需要进行定阶。
比说：随机产生时间序列$X_t+0.6X_{t-1}+0.2X_{t-2}={\xi }_t$中
$p=2$,则该序列为$AR(2)$序列。

2.$\psi =({\psi }_1,{\psi }_2,\cdots ,{\psi }_p{)}^T$为未知参数，需要对其进行估计。

（2）$MA(q)$序列

同样：
1.$q$为模型阶数，需要进行定阶。
2.$\theta =({\theta }_1,{\theta }_2,\cdots ,{\theta }_q{)}^T$为未知参数，需要对其进行估计。

（3）$ARMA(p,q)$序列

1.$p,q$为模型阶数，需要进行定阶。
2.$\theta =({\theta }_1,{\theta }_2,\cdots ,{\theta }_q{)}^T,\psi =({\psi }_1,{\psi }_2,\cdots ,{\psi }_p{)}^T$为未知参数，需要对其进行估计。

平稳时间序列建模流程

下面逐步分析流程并演示代码：

（1）平稳性检验

先导入一组数据，生成图像：

由图可知，其变动范围较大，均值可能为一常数（肉眼判断）在定义上有可能满足平稳序列的条件。
于是，引出了平稳性检验的第一种方法：$1.时序图检验$。
除了这种方法，更一般的方法为$2.自相关检验$
平稳序列通常具有短期相关性，该性质用自相关系数来描述就是：$随着延迟期数的增加，平稳序列的自相关系数会很快衰减为0$。

MATLAB检验自相关系数的函数为：

figure(1)
r11=autocorr(a)   %计算自相关函数
stem(r11)%绘制经线图
title('自相关系数')
figure(2)
r12=parcorr(a)   %计算偏相关函数
stem(r12)%绘制经线图
title('偏相关系数')

对源数据计算可得图像：

平稳系序列的自相关系数具有：截尾性、拖尾性、指数衰减性。根据以上两图可以进行直观理解。从以上图示，$该序列都满足平稳性的判别标准$。

但如上图所示，该序列的自相关系数分布图明显不展现截尾性、拖尾性、指数衰减性的性质，该序列为不平稳序列。
$3.Daniel检验法$。

MATLAB程序：

x0=[80.8 94.0 88.4 101.5 110.3 121.5 134.7 142.7]; 
n=length(x0);
[xsort,ind]=sort(x0) 
Rt(ind)=1:n %计算秩
t=1:n;
Qs=1-6/(n*(n^2-1))*sum((t-Rt).^2) %计算 Qs 的值
T=Qs*sqrt(n-2)/sqrt(1-Qs^2) %计算 T 统计量的值
t_0=tinv(0.975,n-2) %计算上 alpha/2 分位点

计算结果：
$T=11.0235>t_0=2.4469$该序列为$非平稳性序列$。
其次，若已知序列的类型($AR(p)$序列、$MA(q)$序列、$ARMA(p,q)$序列），可采用$4.单位根法$进行平稳性检验。
$\theta (B)及\psi （B）$的根均在单位圆内，则该序列平稳。
例子：已知随机随机时间序列：
$X_t+0.6X_{t-1}+0.2X_{t-2}={\xi }_t$
则$\theta (B)=1+0.6B+0.2B^2$

roots([0.2 0.6 1])

（2）构造差分序列

需要指出，并非是所有的非平稳序列都能进行差分使其平稳化，能进行差分的序列须为$带有一定趋势$的非平稳差分序列。

da=diff(a);      %计算1阶差分
r21=autocorr(da)  %计算自相关函数
r22=parcorr(da)   %计算偏相关函数
figure(1)
stem(r21)
title('1阶差分自相关')
figure(2)
stem(r22)
title('1阶差分自相关')

以上两幅图为原平稳序列差分后的结果，由图可知平稳序列的差分序列不会丧失其平稳特性，反倒越加平稳。若1阶差分满足不了，可采用2阶甚至多阶。

（3）模型定阶与拟合（参数估计）

一般地，若已知该序列的类型，可采用直接指定的方式定阶。还是这个例子：
$X_t+0.6X_{t-1}+0.2X_{t-2}={\xi }_t$

p=3;%指定大概率会的模型阶数
elps=randn(10000,1); x(1:2)=0;
for i=3:10000
    x(i)=-0.6*x(i-1)-0.2*x(i-2)+elps(i); %产生模拟数据
end
x=x';   %x必须为列向量
for i=1:p
    model{i}=ar(x,i);%采用工具箱里面的AR函数进行拟合
AIC(i)=aic(model{i});%根据AIC准则，确定最终
end
N=min(AIC);
P_great=find(AIC==N);

更一般地，对于一组未知类型的时间序列（已知某化工生产浓度数据、已知1997-2012年我国人均国内生产总值、1974-1981年布料销售量$\cdots$，可采用试探模型（即通过构造{$p=1:m,q=1:n$}的方法)进行定阶。这里，我们可以假设模型为计量经济学领域的CARCH模型：

这里我们需要用到MATLAB的两个常用模型指定函数：$$model=arima{(}^{\prime }MALag{s}^{\prime },[MA模型滞后项]{,}^{\prime }ARLag{s}^{\prime },[AR模型滞后项]{,}^{\prime }M{A}^{\prime },MA模型滞后项系数{,}^{\prime }A{R}^{\prime },MA模型滞后项系数{,}^{\prime }Constan{t}^{\prime },常数项）$$
$X_t+0.6X_{t-1}+0.2X_{t-2}={\xi }_t$
可用如下语句进行指定：

model=arima('ARLags',[1 2 3],'MAlags',1','Constant',0)

另外，指定模型之后需要进行拟合，具体函数如下：

输入：Md1为上面用arima指定的模型；y为需要拟合的时间序列；
输出：EstMd1为拟合得到的模型方程；EstParamCov为返回与估计参数；logL为目标函数值；info为汇总信息。
一段完整的定阶与拟合MATLAB程序如下：

for i=0:p
    for j=0:q
        if i==0 && j==0
            continue
        elseif i==0
            ToEstMd=arima('MALags',1:j);
        elseif j==0
            ToEstMd=arima('ARLags',1:i);
        else
         ToEstMd=arima('ARLags',1:i,'MALags',1:j);   
        end
        [EstMd,EstParamCov,logL,info]=estimate(ToEstMd,d);%模型拟合
        numParams=sum(any(EstParamCov));%计算拟合参数个数
        [aic(i,j),bic(i,j)]=aicbic(logL,numParams,n);
    end
end
N=min(min(bic));
[p_great,q_great]=find(N==bic);

（4）模型预报

这里也将有大篇幅的理论推导，这里将不再叙述。
AR预报：
ARMA预报：

 forecast_=forecast(EstMd,5,'Y0',d)   

这里，EstMd为拟合出来的模型，5为需要预报的值，’Y0‘为预测参量，d为序列。

（5）模型检验

模型检验常用的方法有：
时间序列检验
残差的 Portmanteau 检验
Lagrange 乘数检验
MATLAB残差检验的程序如下：

 res=infer(EstMd,d);%计算残差
 h=lbqtest(res);%模型检验