平稳时间序列以及MATLAB相关工具箱学习笔记
概念
(1)平稳序列
即序列的均值是个常数,与序列长度、起始位置无关。
直观看上去,该序列类似于围绕某一个值上下波动(该值为平稳序列均值)。
(2)平稳白噪声序列
所谓噪声,可以理解为一种非周期性的扰动。由于其自协方差函数值为0,于是其各序列值没有任何相关关系,所以说平稳白噪声序列是一段没有记忆的平稳序列。
在MATLAB中可用如下代码生成高斯分布下的平稳白噪声序列:
elps=randn(1,1000);
plot(elps)
基本的平稳时间序列
(1) A R ( p ) AR(p) AR(p)序列
需要说明的是:
1.
p
p
p为模型阶数,需要进行定阶。
比说:随机产生时间序列
X
t
+
0.6
X
t
−
1
+
0.2
X
t
−
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,则该序列为
A
R
(
2
)
AR(2)
AR(2)序列。
2. ψ = ( ψ 1 , ψ 2 , ⋯ , ψ p ) T \psi=(\psi_1,\psi_2,\cdots,\psi_p)^T ψ=(ψ1,ψ2,⋯,ψp)T为未知参数,需要对其进行估计。
(2)
M
A
(
q
)
MA(q)
MA(q)序列![在这里插入图片描述](https://i-blog.csdnimg.cn/blog_migrate/3dc38510eff130124d25a22716487c0e.png)
同样:
1.
q
q
q为模型阶数,需要进行定阶。
2.
θ
=
(
θ
1
,
θ
2
,
⋯
,
θ
q
)
T
\theta=(\theta_1,\theta_2,\cdots,\theta_q)^T
θ=(θ1,θ2,⋯,θq)T为未知参数,需要对其进行估计。
(3) A R M A ( p , q ) ARMA(p,q) ARMA(p,q)序列
1.
p
,
q
p,q
p,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.
D
a
n
i
e
l
检
验
法
\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
>
t
0
=
2.4469
T = 11.0235> t_0 =2.4469
T=11.0235>t0=2.4469该序列为
非
平
稳
性
序
列
\color{#D1F}{非平稳性序列}
非平稳性序列。
其次,若已知序列的类型(
A
R
(
p
)
AR(p)
AR(p)序列、
M
A
(
q
)
MA(q)
MA(q)序列、
A
R
M
A
(
p
,
q
)
ARMA(p,q)
ARMA(p,q)序列),可采用
4.
单
位
根
法
\color{#00F}{4.单位根法}
4.单位根法进行平稳性检验。
θ
(
B
)
及
ψ
(
B
)
\theta(B)及\psi(B)
θ(B)及ψ(B)的根均在单位圆内,则该序列平稳。
例子:已知随机随机时间序列:
X
t
+
0.6
X
t
−
1
+
0.2
X
t
−
2
=
ξ
t
X_t+0.6X_{t-1}+0.2X_{t-2}=\xi_t
Xt+0.6Xt−1+0.2Xt−2=ξt
则
θ
(
B
)
=
1
+
0.6
B
+
0.2
B
2
\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)模型定阶与拟合(参数估计)
一般地,若已知该序列的类型,可采用直接指定的方式定阶。还是这个例子:
X
t
+
0.6
X
t
−
1
+
0.2
X
t
−
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
:
n
p=1:m,q=1:n
p=1:m,q=1:n}的方法)进行定阶。这里,我们可以假设模型为计量经济学领域的CARCH模型:
这里我们需要用到MATLAB的两个常用模型指定函数:
m
o
d
e
l
=
a
r
i
m
a
(
′
M
A
L
a
g
s
′
,
[
M
A
模
型
滞
后
项
]
,
′
A
R
L
a
g
s
′
,
[
A
R
模
型
滞
后
项
]
,
′
M
A
′
,
M
A
模
型
滞
后
项
系
数
,
′
A
R
′
,
M
A
模
型
滞
后
项
系
数
,
′
C
o
n
s
t
a
n
t
′
,
常
数
项
)
model=arima('MALags',[MA模型滞后项],'ARLags',[AR模型滞后项],'MA',{MA模型滞后项系数},'AR',{MA模型滞后项系数},'Constant',常数项)
model=arima(′MALags′,[MA模型滞后项],′ARLags′,[AR模型滞后项],′MA′,MA模型滞后项系数,′AR′,MA模型滞后项系数,′Constant′,常数项)
X
t
+
0.6
X
t
−
1
+
0.2
X
t
−
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);%模型检验