平稳时间序列以及MATLAB相关工具箱学习笔记
概念
(1)平稳序列
即序列的均值是个常数,与序列长度、起始位置无关。
直观看上去,该序列类似于围绕某一个值上下波动(该值为平稳序列均值)。
(2)平稳白噪声序列
所谓噪声,可以理解为一种非周期性的扰动。由于其自协方差函数值为0,于是其各序列值没有任何相关关系,所以说平稳白噪声序列是一段没有记忆的平稳序列。
在MATLAB中可用如下代码生成高斯分布下的平稳白噪声序列:
elps=randn(1,1000);
plot(elps)
基本的平稳时间序列
(1)AR(p)AR(p)AR(p)序列
需要说明的是:
1.ppp为模型阶数,需要进行定阶。
比说:随机产生时间序列Xt+0.6Xt−1+0.2Xt−2=ξt\color{#00F}{X_t+0.6X_{t-1}+0.2X_{t-2}=\xi_t}Xt+0.6Xt−1+0.2Xt−2=ξt中
p=2\color{#00F}{p=2}p=2,则该序列为AR(2)AR(2)AR(2)序列。
2.ψ=(ψ1,ψ2,⋯ ,ψp)T\psi=(\psi_1,\psi_2,\cdots,\psi_p)^Tψ=(ψ1,ψ2,⋯,ψp)T为未知参数,需要对其进行估计。
(2)MA(q)MA(q)MA(q)序列
同样:
1.qqq为模型阶数,需要进行定阶。
2.θ=(θ1,θ2,⋯ ,θq)T\theta=(\theta_1,\theta_2,\cdots,\theta_q)^Tθ=(θ1,θ2,⋯,θq)T为未知参数,需要对其进行估计。
(3)ARMA(p,q)ARMA(p,q)ARMA(p,q)序列
1.p,qp,qp,q为模型阶数,需要进行定阶。
2.θ=(θ1,θ2,⋯ ,θq)T,ψ=(ψ1,ψ2,⋯ ,ψp)T\theta=(\theta_1,\theta_2,\cdots,\theta_q)^T,\psi=(\psi_1,\psi_2,\cdots,\psi_p)^Tθ=(θ1,θ2,⋯,θq)T,ψ=(ψ1,ψ2,⋯,ψp)T为未知参数,需要对其进行估计。
平稳时间序列建模流程
下面逐步分析流程并演示代码:
(1)平稳性检验
先导入一组数据,生成图像:
由图可知,其变动范围较大,均值可能为一常数(肉眼判断)在定义上有可能满足平稳序列的条件。
于是,引出了平稳性检验的第一种方法:1.时序图检验\color{#00F}{1.时序图检验}1.时序图检验。
除了这种方法,更一般的方法为2.自相关检验\color{#00F}{2.自相关检验}2.自相关检验
平稳序列通常具有短期相关性,该性质用自相关系数来描述就是:随着延迟期数的增加,平稳序列的自相关系数会很快衰减为0\color{#F00}{随着延迟期数的增加,平稳序列的自相关系数会很快衰减为0}随着延迟期数的增加,平稳序列的自相关系数会很快衰减为0。
MATLAB检验自相关系数的函数为:
figure(1)
r11=autocorr(a) %计算自相关函数
stem(r11)%绘制经线图
title('自相关系数')
figure(2)
r12=parcorr(a) %计算偏相关函数
stem(r12)%绘制经线图
title('偏相关系数')
对源数据计算可得图像:
平稳系序列的自相关系数具有:截尾性、拖尾性、指数衰减性。根据以上两图可以进行直观理解。从以上图示,该序列都满足平稳性的判别标准\color{#00F}{该序列都满足平稳性的判别标准}该序列都满足平稳性的判别标准。
但如上图所示,该序列的自相关系数分布图明显不展现截尾性、拖尾性、指数衰减性的性质,该序列为不平稳序列。
3.Daniel检验法\color{#00F}{3.Daniel检验法}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>t0=2.4469T = 11.0235> t_0 =2.4469T=11.0235>t0=2.4469该序列为非平稳性序列\color{#D1F}{非平稳性序列}非平稳性序列。
其次,若已知序列的类型(AR(p)AR(p)AR(p)序列、MA(q)MA(q)MA(q)序列、ARMA(p,q)ARMA(p,q)ARMA(p,q)序列),可采用4.单位根法\color{#00F}{4.单位根法}4.单位根法进行平稳性检验。
θ(B)及ψ(B)\theta(B)及\psi(B)θ(B)及ψ(B)的根均在单位圆内,则该序列平稳。
例子:已知随机随机时间序列:
Xt+0.6Xt−1+0.2Xt−2=ξtX_t+0.6X_{t-1}+0.2X_{t-2}=\xi_tXt+0.6Xt−1+0.2Xt−2=ξt
则θ(B)=1+0.6B+0.2B2\theta(B)=1+0.6B+0.2B^2θ(B)=1+0.6B+0.2B2
roots([0.2 0.6 1])
(2)构造差分序列
需要指出,并非是所有的非平稳序列都能进行差分使其平稳化,能进行差分的序列须为带有一定趋势\color{#F00}{带有一定趋势}带有一定趋势的非平稳差分序列。
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)模型定阶与拟合(参数估计)
一般地,若已知该序列的类型,可采用直接指定的方式定阶。还是这个例子:
Xt+0.6Xt−1+0.2Xt−2=ξt\color{#00F}{X_t+0.6X_{t-1}+0.2X_{t-2}=\xi_t}Xt+0.6Xt−1+0.2Xt−2=ξ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:np=1:m,q=1:np=1:m,q=1:n}的方法)进行定阶。这里,我们可以假设模型为计量经济学领域的CARCH模型:
这里我们需要用到MATLAB的两个常用模型指定函数:model=arima(′MALags′,[MA模型滞后项],′ARLags′,[AR模型滞后项],′MA′,MA模型滞后项系数,′AR′,MA模型滞后项系数,′Constant′,常数项)model=arima('MALags',[MA模型滞后项],'ARLags',[AR模型滞后项],'MA',{MA模型滞后项系数},'AR',{MA模型滞后项系数},'Constant',常数项)model=arima(′MALags′,[MA模型滞后项],′ARLags′,[AR模型滞后项],′MA′,MA模型滞后项系数,′AR′,MA模型滞后项系数,′Constant′,常数项)
Xt+0.6Xt−1+0.2Xt−2=ξt\color{#00F}{X_t+0.6X_{t-1}+0.2X_{t-2}=\xi_t}Xt+0.6Xt−1+0.2Xt−2=ξ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);%模型检验