概述
第四章 平稳序列的拟合与预测
4.1 建模步骤
当我们拿到一组数据,想要通过这组数据来探究过去发生的规律,研究现在发生的强度并预测未来事情的发展走向,就需要从这些数据中挖掘出事物发生的特征,比如频率、强度、集中值等等,来提取“经验”,对过去的作出反应,弥补或者提取收益,而后建立一个合适的模型,来预测未来事物的发展,掌握未来动态,并及时调整,做好准备。
对于时间序列数据,当我们拿到数据时,我们首先要对数据进行预处理,事实上,所有的数据都是如此,经过预处理使之成为合适的,可以进行预测的,有用的数据。
时间序列数据的预处理是要进行平稳性检验,只有平稳的时间序列数据才是可以利用的,有效的数据,才能进行预测。还要进行纯随机检验,确保其是非白噪声序列。
4.2 模型识别
若数据是平稳的非白噪声序列,则可以考虑建模,但时间序列模型不止一种,具体要用哪种模型要根据数据的特征来判断,即先计算样本的ACF和PACF。根据样本的自相关系数和偏自相关系数的截尾和拖尾性来选择合适的模型拟合观察值序列。
定阶的基本原则如下:
ρ k ^ hat{rho_k} ρk^ | ϕ k k ^ hat{phi_{kk}} ϕkk^ | 模型定阶 |
---|---|---|
拖尾 | p阶截尾 | AR(p)模型 |
q阶截尾 | 拖尾 | MA(q)模型 |
拖尾 | 拖尾 | ARMA(p,q)模型 |
但是由于样本的随机性,样本的相关系数不会呈现出理论上的截尾的完美情况,应该截尾的样本自相关系数或偏自相关系数就会呈现出小值震荡。同时,又因为平稳时间序列通常具有的短期相关性,随着延迟阶数趋近于无穷, ρ k ^ hat{rho_k} ρk^与 ϕ k k ^ hat{phi_{kk}} ϕkk^都会衰减至零值附近作小值波动。
这就给我们留下了一个问题,判断拖尾的标准情况是什么?其实并没有,主要还是依靠分析人员的主观经验。根据诸多统计学家的证明,样本自相关系数是总体自相关系数的有偏估计,当k足够大时,根据平稳序列自相关系数呈负指数衰减,有 ρ k → 0 rho_krightarrow0 ρk→0,而当样本容量n充分大时,样本的自相关系数和偏自相关系数都是近似于正态分布,根据正态分布的性质,可以利用2倍标准差范围辅助判断。
如果样本自相关系数或偏自相关系数在最初的d阶明显超过2倍标准差范围,而后几乎95%的自相关系数都落在2倍标准差的范围以内,而且由非零自相关系数衰减为小值波动的过程非常突然,这时,通常视为自相关系数截尾。截尾阶数为d。
如果有超过5%的样本自相关系数落入2倍标准差范围之外,或者由显著非零的自相关系数衰减为小值波动的过程非常连续或比较缓慢,这时,通常视为自相关系数不截尾。
4.3 参数估计
对于非中心化的ARMA(p,q)模型,有 x t = μ + Θ q ( B ) Φ q ( B ) ϵ t x_t=mu+frac{Theta_q(B)}{Phi_q(B)}epsilon_t xt=μ+Φq(B)Θq(B)ϵt,式中 ϵ t ∼ W N ( 0 , σ ϵ 2 ) epsilon_tthicksim WN(0,sigma^2_{epsilon}) ϵt∼WN(0,σϵ2)【WN白噪声】。模型中共含有p+q+2个未知参数: θ 1 , θ 2 , ⋯ , θ q , ϕ 1 , ϕ 2 , ⋯ , ϕ p , μ , σ ϵ 2 theta_1,theta_2,cdots,theta_q,phi_1,phi_2,cdots,phi_p,mu,sigma_{epsilon}^2 θ1,θ2,⋯,θq,ϕ1,ϕ2,⋯,ϕp,μ,σϵ2。
参数 μ mu μ是序列均值,通常采用一阶中心距的估计方法,用样本均值估计总体均值即可得到它的估计值。
再将原序列中心化,则待估参数减少至p+q+1个。对这些参数有三种估计方法:距估计、极大似然估计、最小二乘估计。
4.3.1 距估计
用p+q个样本自相关系数估计总体自相关系数,构造p+q个方程组成的Yule-Walker方程组
{
ρ
1
(
ϕ
1
,
⋯
,
ϕ
p
,
θ
1
,
⋯
,
θ
q
)
=
ρ
1
^
⋮
ρ
p
+
q
(
ϕ
1
,
⋯
,
ϕ
p
,
θ
1
,
⋯
,
θ
q
)
=
ρ
p
+
q
^
begin{cases} rho_1(phi_1,cdots,phi_p,theta_1,cdots,theta_q)=hat{rho_1}\ vdots\ rho_{p+q}(phi_1,cdots,phi_p,theta_1,cdots,theta_q)=hat{rho_{p+q}}\ end{cases}
⎩⎪⎪⎨⎪⎪⎧ρ1(ϕ1,⋯,ϕp,θ1,⋯,θq)=ρ1^⋮ρp+q(ϕ1,⋯,ϕp,θ1,⋯,θq)=ρp+q^
从中解出的参数值
ϕ
1
^
,
⋯
,
ϕ
p
^
,
θ
1
^
,
⋯
,
θ
q
^
hat{phi_1},cdots,hat{phi_p},hat{theta_1},cdots,hat{theta_q}
ϕ1^,⋯,ϕp^,θ1^,⋯,θq^就是
ϕ
1
,
⋯
,
ϕ
p
,
θ
1
,
⋯
,
θ
q
phi_1,cdots,phi_p,theta_1,cdots,theta_q
ϕ1,⋯,ϕp,θ1,⋯,θq的距估计。
尤尔—沃克方程(Yule-Walker equation)是描述自回归序列参数与其协方差函数之间关系的方程。
所以,距估计其实就是建立多个自相关系数与参数的方程,而后解出参数(因为 ρ rho ρ可估或者说是已知)。
4.3.2 极大似然估计
在极大似然准则下(事件既然发生,则其概率就一定尽可能大),认为样本来自使该样本出现大概率最大的总体。因此未知参数的极大似然估计就是使得似然函数(即联合密度函数)达到最大的参数值。
L
(
β
1
^
,
⋯
,
β
k
^
;
x
1
,
⋯
,
x
n
)
=
m
a
x
{
p
(
x
1
,
⋯
,
x
n
)
;
β
1
,
⋯
,
β
k
}
L(hat{beta_1},cdots,hat{beta_k};x_1,cdots,x_n)=max{p(x_1,cdots,x_n);beta_1,cdots,beta_k}
L(β1^,⋯,βk^;x1,⋯,xn)=max{p(x1,⋯,xn);β1,⋯,βk}
但是,使用极大似然估计必须已知总体的分布函数,而在时间序列分析中,序列总体分布通常是未知的。为了便于分析计算,通常假设序列服从多元正态分布。
独立的多元正态分布:
以上在理论上是可行的,但实际上上式存在p+q个超越方程,通常需要利用迭代法才能求出未知参数的极大似然估计值。实际操作中可以借助计算机。
4.3.3 最小二乘估计
在ARMA(p,q)模型场合,记
β
~
=
(
ϕ
1
,
⋯
,
ϕ
p
,
θ
1
,
⋯
,
θ
q
)
′
F
t
(
β
~
)
=
ϕ
1
x
t
−
1
+
⋯
+
ϕ
p
x
t
−
p
−
θ
1
ϵ
t
−
1
−
⋯
−
θ
q
ϵ
t
−
q
tilde{beta}=(phi_1,cdots,phi_p,theta_1,cdots,theta_q)'\ F_t(tilde{beta})=phi_1x_{t-1}+cdots+phi_px_{t-p}-theta_1epsilon_{t-1}-cdots-theta_qepsilon_{t-q}
β~=(ϕ1,⋯,ϕp,θ1,⋯,θq)′Ft(β~)=ϕ1xt−1+⋯+ϕpxt−p−θ1ϵt−1−⋯−θqϵt−q
残差项为:
ϵ
t
=
x
t
−
F
t
(
β
~
)
epsilon_t=x_t-F_t(tilde{beta})
ϵt=xt−Ft(β~)
残差平方和为:
Q
(
β
~
)
=
∑
t
=
1
n
ϵ
t
2
=
∑
t
=
1
n
(
x
t
−
ϕ
1
x
t
−
1
+
⋯
+
ϕ
p
x
t
−
p
−
θ
1
ϵ
t
−
1
−
⋯
−
θ
q
ϵ
t
−
q
)
2
Q(tilde{beta})=sum_{t=1}^{n}epsilon_t^2\ =sum_{t=1}^{n}(x_t-phi_1x_{t-1}+cdots+phi_px_{t-p}-theta_1epsilon_{t-1}-cdots-theta_qepsilon_{t-q})^2
Q(β~)=t=1∑nϵt2=t=1∑n(xt−ϕ1xt−1+⋯+ϕpxt−p−θ1ϵt−1−⋯−θqϵt−q)2
使残差平方和达到最小的那组参数值即
β
~
tilde{beta}
β~的最小二乘估计。
在实际中,最常用的是条件最小二乘估计方法。它假定过去未观测到的序列值等于零,于是有
ϵ
t
=
Φ
(
B
)
Θ
(
B
)
x
t
=
x
t
−
∑
i
=
1
t
π
i
x
t
−
i
epsilon_t=frac{Phi(B)}{Theta(B)}x_t=x_t-sum_{i=1}^{t}pi_ix_{t-i}
ϵt=Θ(B)Φ(B)xt=xt−i=1∑tπixt−i
于是残差平方和为:
Q
(
β
~
)
=
∑
t
=
1
n
ϵ
t
2
=
∑
t
=
1
n
[
x
t
−
∑
i
=
1
t
π
i
x
t
−
i
]
2
Q(tilde{beta})=sum_{t=1}^{n}epsilon_t^2\ =sum_{t=1}^{n}[x_t-sum_{i=1}^{t}pi_ix_{t-i}]^2
Q(β~)=t=1∑nϵt2=t=1∑n[xt−i=1∑tπixt−i]2
4.4 模型检验
4.4.1 模型的显著性检验
模型的显著性检验主要是检验模型的有效性。且一个模型是否显著有效主要看它提取的信息是否充分。也即是残差项应该是白噪声序列。
4.4.2 参数的显著性检验
参数的显著性检验就是检验每一个未知参数是否显著非零。
4.5 模型优化
有时同一个序列可以有不同的拟合模型,且都显著有效,那么我们怎么选择呢?
这时就要用到AIC准则和SBC准则,即自小信息量准则,一般用软件拟合时会直接给出。
A
I
C
=
−
2
l
n
(
模
型
的
极
大
似
然
函
数
值
)
+
2
(
模
型
中
未
知
参
数
个
数
)
AIC=-2ln(模型的极大似然函数值)+2(模型中未知参数个数)
AIC=−2ln(模型的极大似然函数值)+2(模型中未知参数个数)
AIC准则的局限性在于受到样本容量的限制,当样本容量趋于无限大的时候,由AIC准则选择的模型不收敛于真实模型,它通常比真实模型所含的未知参数个数还要多。
BIC准则也称SBC准则弥补了这个不足。
S
B
C
=
−
2
l
n
(
模
型
的
极
大
似
然
函
数
值
)
+
(
l
n
n
)
(
模
型
中
未
知
参
数
个
数
)
SBC=-2ln(模型的极大似然函数值)+(ln n)(模型中未知参数个数)
SBC=−2ln(模型的极大似然函数值)+(lnn)(模型中未知参数个数)
4.6 序列预测
根据ARMA(p,q)模型的平稳性和可逆性,可以用传递形式和逆转形式等价描述该序列:
x
t
=
∑
i
=
0
∞
G
i
ϵ
t
−
i
ϵ
t
=
∑
j
=
0
∞
I
j
x
t
−
j
x_t=sum_{i=0}^{infty}{G_iepsilon_{t-i}}\ epsilon_t=sum_{j=0}^{infty}{I_jx_{t-j}}
xt=i=0∑∞Giϵt−iϵt=j=0∑∞Ijxt−j
其中
{
G
i
}
{G_i}
{Gi}是Green函数值;
{
I
j
}
{I_j}
{Ij}为逆转函数值。
将两式合并,有
x
t
=
∑
i
=
0
∞
∑
j
=
0
∞
G
i
I
j
x
t
−
i
−
j
x_t=sum_{i=0}^{infty}sum_{j=0}^{infty}G_iI_jx_{t-i-j}
xt=i=0∑∞j=0∑∞GiIjxt−i−j
简记为
x
t
=
∑
i
=
0
∞
C
i
x
t
−
1
−
i
x_t=sum_{i=0}^{infty}C_ix_{t-1-i}
xt=i=0∑∞Cixt−1−i
∀
l
⩾
1
forall lgeqslant1
∀l⩾1时,对任意未来t+l时刻,该时刻的序列值可以表示为它的历史数据
x
t
+
l
−
1
,
⋯
,
x
t
+
1
,
x
t
,
⋯
x_{t+l-1},cdots,x_{t+1},x_t,cdots
xt+l−1,⋯,xt+1,xt,⋯的线性函数
x
t
+
l
=
∑
i
=
0
∞
C
i
x
t
+
l
−
1
−
i
x_{t+l}=sum_{i=0}^{infty}{C_i}x_{t+l-1-i}
xt+l=i=0∑∞Cixt+l−1−i
但有趣的是
x
t
x_t
xt前的历史信息已知,
x
t
+
1
,
⋯
,
x
t
+
l
−
1
x_{t+1},cdots,x_{t+l-1}
xt+1,⋯,xt+l−1期的信息是未知的,但根据线性函数的可加性,所有的未知历史信息都可以用已知历史信息的线性函数表示出来。
如
x
t
+
2
=
∑
i
=
0
∞
C
i
x
t
+
1
−
i
=
C
0
x
t
+
1
+
∑
i
=
0
∞
C
i
+
1
x
t
−
i
x_{t+2}=sum_{i=0}^{infty}C_ix_{t+1-i}=C_0x_{t+1}+sum_{i=0}^{infty}C_{i+1}x_{t-i}
xt+2=∑i=0∞Cixt+1−i=C0xt+1+∑i=0∞Ci+1xt−i,但是因为
x
t
+
1
x_{t+1}
xt+1是未知的,而其又可以表示为
x
t
+
1
=
∑
i
=
0
∞
C
i
x
t
−
i
x_{t+1}=sum_{i=0}^{infty}C_ix_{t-i}
xt+1=∑i=0∞Cixt−i,故有
x
t
+
2
=
∑
i
=
0
∞
(
C
0
C
i
+
C
i
+
1
)
x
t
−
i
x_{t+2}=sum_{i=0}^{infty}(C_0C_i+C_{i+1})x_{t-i}
xt+2=i=0∑∞(C0Ci+Ci+1)xt−i
由此
x
t
+
2
x_{t+2}
xt+2最终表示为已知历史信息
x
t
,
x
t
−
1
,
⋯
x_t,x_{t-1},cdots
xt,xt−1,⋯的线性函数。同理,对于未来任意l时刻的序列值
x
t
+
l
x_{t+l}
xt+l最终都可以用已知历史信息
x
t
,
x
t
−
1
,
⋯
x_t,x_{t-1},cdots
xt,xt−1,⋯表示出来,并将该线性函数形式记为:
x
t
^
(
l
)
=
∑
i
=
0
∞
D
^
i
x
t
−
i
hat{x_t}(l)=sum_{i=0}^{infty}{hat{D}_ix_{t-i}}
xt^(l)=i=0∑∞D^ixt−i
x
^
t
(
l
)
hat{x}_t(l)
x^t(l)也称为序列
{
x
t
}
{x_t}
{xt}的第l步预测值。
代码:
A7<-read.csv(“A1_7.CSV”)
A8<-read.csv(“A1_8.csv”)
A9<-read.csv(“A1_9.csv”)
a7<-ts(A7
n
u
m
b
e
r
,
s
t
a
r
t
=
1900
,
f
r
e
q
u
e
n
c
y
=
1
)
a
8
<
−
t
s
(
A
8
number,start = 1900,frequency = 1) a8<-ts(A8
number,start=1900,frequency=1)a8<−ts(A8overshort,start = 1,frequency = 365)
a9<-ts(A9$change,start = 1880,frequency = 1)
a<-diff(a9,lag = 1)#对数据进行一阶差分
#install.packages(“aTSA”)
library(aTSA)
#class(A7)
adf.test(a7)#平稳性检验
acf(a7,lag=25)#绘制自相关图
acf(a7,lag=25)KaTeX parse error: Expected 'EOF', got '#' at position 4: acf#̲显示自相关系数 拖尾 pac…pacf#显示偏自相关系数 一阶截尾
plot(a8)
adf.test(a8)#ADF检验
Box.test(a8,lag = 6)
Box.test(a8,lag = 12)
for (k in 1:2) {print(Box.test(a8,lag = 6k,type = “Ljung-Box”))
}
acf(a8,lag=20)
pacf(a8)
plot(a9)
adf.test(a9)
plot(a)
acf(a)
pacf(a)
for (k in 1:2) {print(Box.test(a9,lag = 6k,type = “Ljung-Box”))
}
par(mfrow=c(1,2))
acf(a9)
pacf(a9)
###参数估计,使用函数arima()拟合
fit1<-arima(a7,order = c(1,0,0),include.mean=T,method = “CSS-ML”)#极大似然估计与条件最小二乘估计混合
fit1
fit1<-arima(a7,order = c(1,0,0),method = “CSS”)
fit2<-arima(a8,order = c(0,0,1),method = “CSS”)#条件最小二乘估计
fit2
fit3<-arima(a,order = c(1,0,1))
fit3
#####模型检验,使用aTSA包中的ts.diag函数
library(aTSA)
ts.diag(fit1)#模型的显著性检验
t=abs(fit1
c
o
e
f
)
/
s
q
r
t
(
d
i
a
g
(
f
i
t
1
coef)/sqrt(diag(fit1
coef)/sqrt(diag(fit1var.coef))
pt(t,length(a7)-length(fit1$coef),lower.tail = F)
###模型预测,使用forecast包中的forecast函数
install.packages(“forecast”)
#library(forecast)
forecast::forecast##因前面打开aTSA包
fore1<-forecast::forecast(fit1,h=10)
fore1
plot(fore1,lty=2)
lines(fore1$fitted,col=4)
未完…
最后
以上就是心灵美羽毛为你收集整理的第四章 平稳序列的拟合与预测第四章 平稳序列的拟合与预测的全部内容,希望文章能够帮你解决第四章 平稳序列的拟合与预测第四章 平稳序列的拟合与预测所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。
发表评论 取消回复