文章目录
ARCH模型以及编程实现
ARCH模型的引入
传统的ARIMA模型拟合非平稳时间序列,残差序列 { ε t } \{\varepsilon_t\} {εt}通常满足下面的三个假定条件
-
零均值
E ( ε t ) = 0 E(\varepsilon_t)=0 E(εt)=0 -
纯随机
C o v ( ε t , ε t − i ) = 0 , ∀ i ≥ 1 Cov(\varepsilon_t,\varepsilon_{t-i})=0,\forall i\ge1 Cov(εt,εt−i)=0,∀i≥1 -
方差齐性
V a r ( ε t ) = σ t 2 Var(\varepsilon_t)=\sigma_t^2 Var(εt)=σt2
传统的ARIMA模型假定时间序列的长期方差是相同的,但是,通常情况下,时间序列的存在异方差的,即
V
a
r
(
ε
t
)
=
h
(
t
)
Var(\varepsilon_t)=h(t)
Var(εt)=h(t)
根据以上条件,我们可以得出,残差序列的方差就是它平方的期望,即
V
a
r
(
ε
t
)
=
E
(
ε
t
2
)
Var(\varepsilon_t)=E(\varepsilon_t^2)
Var(εt)=E(εt2)
我们可以根据残差平方图 ε t \varepsilon_t εt关于 t t t变化的二维坐标图,对方差齐性直观判断。
我们对存在异方差的时间序列有两种处理方式
1.如果已知异方差函数具体形式,进行方差齐性变化
2.如果不知异方差函数具体形式,拟合条件异方差模型,即ARCH
异方差性检验
Q检验
原假设
H
0
:
H_0:
H0:残差平方序列纯随机,即
ρ
1
=
ρ
2
=
⋯
=
ρ
q
=
0
\rho_1=\rho_2=\dots=\rho_q=0
ρ1=ρ2=⋯=ρq=0
备择假设
H
1
:
H_1:
H1:残差平方序列具有自相关性,即
ρ
1
,
ρ
2
,
…
,
ρ
q
\rho_1,\rho_2,\dots,\rho_q
ρ1,ρ2,…,ρq不全为
0
0
0
检验统计量
Q
(
q
)
=
n
(
n
+
2
)
∑
i
=
1
q
ρ
i
2
n
−
i
Q(q)=n(n+2)\sum_{i=1}^{q}\frac{\rho_i^2}{n-i}
Q(q)=n(n+2)i=1∑qn−iρi2
其中, n n n为观察序列长度, ρ i \rho_i ρi为残差序列延迟i的自相关系数,设
σ ^ 2 = ∑ i = 1 n ε i 2 n \hat{\sigma}^2=\frac{\sum_{i=1}^n\varepsilon_i^2}{n} σ^2=n∑i=1nεi2
则
ρ
i
=
∑
t
−
i
=
1
n
(
ε
t
2
−
σ
^
2
)
(
ε
t
−
i
2
−
σ
^
2
)
∑
t
=
1
n
(
ε
t
2
−
σ
^
2
)
2
\rho_i=\sqrt{\frac{\sum_{t-i=1}^{n}(\varepsilon_t^2-\hat{\sigma}^2)(\varepsilon_{t-i}^2-\hat{\sigma}^2)}{\sum_{t=1}^n(\varepsilon_t^2-\hat{\sigma}^2)^2}}
ρi=∑t=1n(εt2−σ^2)2∑t−i=1n(εt2−σ^2)(εt−i2−σ^2)
Q统计量近似服从于自由度为 q − 1 q-1 q−1的 χ 2 \chi^2 χ2分布,即 Q ( q ) ∼ χ 2 ( q − 1 ) Q(q)\sim\chi^2(q-1) Q(q)∼χ2(q−1)
若
Q
(
q
)
≥
χ
1
−
α
2
(
q
−
1
)
Q(q)\ge\chi_{1-\alpha}^2(q-1)
Q(q)≥χ1−α2(q−1),则拒绝原假设,认为存在异方差
若
Q
(
q
)
<
χ
1
−
α
2
(
q
−
1
)
Q(q)\lt\chi_{1-\alpha}^2(q-1)
Q(q)<χ1−α2(q−1),则接受原假设,认为不存在异方差
LM检验
原假设
H
0
:
H_0:
H0:残差平方序列纯随机,即
ρ
1
=
ρ
2
=
⋯
=
ρ
q
=
0
\rho_1=\rho_2=\dots=\rho_q=0
ρ1=ρ2=⋯=ρq=0
备择假设
H
1
:
H_1:
H1:残差平方序列具有自相关性,即
ρ
1
,
ρ
2
,
…
,
ρ
q
\rho_1,\rho_2,\dots,\rho_q
ρ1,ρ2,…,ρq不全为0
若
H
0
H_0
H0为真,则
L
M
(
q
)
=
W
T
W
LM(q)=W^TW
LM(q)=WTW
W
=
(
ρ
1
2
σ
2
,
ρ
2
2
σ
2
,
…
,
ρ
q
2
σ
2
)
W=(\frac{\rho_1^2}{\sigma^2},\frac{\rho_2^2}{\sigma^2},\dots,\frac{\rho_q^2}{\sigma^2})
W=(σ2ρ12,σ2ρ22,…,σ2ρq2)
其中,
ρ
i
=
∑
t
=
i
+
1
n
(
ε
t
2
−
σ
^
2
)
(
ε
t
−
i
2
−
σ
^
2
)
∑
t
=
1
n
(
ε
t
2
−
σ
^
2
)
2
,
σ
^
2
=
∑
i
=
1
n
ε
i
2
n
\rho_i=\sqrt{\frac{\sum_{t=i+1}^{n}(\varepsilon_t^2-\hat{\sigma}^2)(\varepsilon_{t-i}^2-\hat{\sigma}^2)}{\sum_{t=1}^n(\varepsilon_t^2-\hat{\sigma}^2)^2}}, \hat{\sigma}^2=\frac{\sum_{i=1}^n\varepsilon_i^2}{n}
ρi=∑t=1n(εt2−σ^2)2∑t=i+1n(εt2−σ^2)(εt−i2−σ^2),σ^2=n∑i=1nεi2
L M ( q ) LM(q) LM(q)统计量近似服从于自由度为 q − 1 q-1 q−1的 χ 2 \chi^2 χ2分布,即 L M ( q ) ∼ χ 2 ( q − 1 ) LM(q)\sim\chi^2(q-1) LM(q)∼χ2(q−1)
若
L
M
(
q
)
≥
χ
1
−
α
2
(
q
−
1
)
LM(q)\ge\chi_{1-\alpha}^2(q-1)
LM(q)≥χ1−α2(q−1),则拒绝原假设,认为存在异方差
若
L
M
(
q
)
<
χ
1
−
α
2
(
q
−
1
)
LM(q)\lt\chi_{1-\alpha}^2(q-1)
LM(q)<χ1−α2(q−1),则接受原假设,认为不存在异方差
A R C H ( q ) ARCH(q) ARCH(q)模型
具有
{
x
t
=
f
(
t
,
x
t
−
1
,
x
t
−
2
,
…
)
+
ε
t
ε
t
=
h
t
e
t
h
t
=
E
(
ε
t
2
)
=
ω
+
∑
j
=
1
q
λ
j
ε
t
−
j
2
\begin{cases} x_t=f(t,x_{t-1},x_{t-2},\dots)+\varepsilon_t& \\\\ \varepsilon_t=\sqrt{h_t}e_t & \\\\ h_t=E(\varepsilon_t^2)=\omega+\sum_{j=1}^q \lambda_j\varepsilon_{t-j}^2 & \end{cases}
⎩
⎨
⎧xt=f(t,xt−1,xt−2,…)+εtεt=htetht=E(εt2)=ω+∑j=1qλjεt−j2
结构的模型称为q阶自回归条件异方差模型,简记为ARCH(q).
G A R C H ( p , q ) GARCH(p,q) GARCH(p,q)模型
具有
{
x
t
=
f
(
t
,
x
t
−
1
,
x
t
−
2
,
…
)
+
ε
t
ε
t
=
h
t
e
t
h
t
=
ω
+
∑
i
=
1
p
η
i
h
t
−
i
+
∑
j
=
1
q
λ
j
ε
t
−
j
2
\begin{cases} x_t=f(t,x_{t-1},x_{t-2},\dots)+\varepsilon_t \\\\ \varepsilon_t=\sqrt{h_t}e_t\\\\ h_t=\omega+\sum_{i=1}^p \eta_ih_{t-i}+\sum_{j=1}^q \lambda_j\varepsilon_{t-j}^2\end{cases}
⎩
⎨
⎧xt=f(t,xt−1,xt−2,…)+εtεt=htetht=ω+∑i=1pηiht−i+∑j=1qλjεt−j2
的模型称为广义自回归条件异方差模型,简记为
G
A
R
C
H
(
p
,
q
)
GARCH(p,q)
GARCH(p,q).
A R ( m ) − G A R C H ( p , q ) AR(m)-GARCH(p,q) AR(m)−GARCH(p,q)模型
{ x t = f ( t , x t − 1 , x t − 2 , … ) + ε t ε t = ∑ k = 1 m β k ε t − k + v t v t = h t e t h t = ω + ∑ i = 1 p η i h t − i + ∑ j = 1 q λ j v t − j 2 \begin{cases}x_t=f(t,x_{t-1},x_{t-2},\dots)+\varepsilon_t\\\\ \varepsilon_t=\sum_{k=1}^m\beta_k\varepsilon_{t-k}+v_t\\\\ v_t=\sqrt{h_t}e_t\\\\ h_t=\omega+\sum_{i=1}^p \eta_ih_{t-i}+\sum_{j=1}^q \lambda_jv_{t-j}^2 \end{cases} ⎩ ⎨ ⎧xt=f(t,xt−1,xt−2,…)+εtεt=∑k=1mβkεt−k+vtvt=htetht=ω+∑i=1pηiht−i+∑j=1qλjvt−j2
GARCH的衍生模型
指数GARCH模型(EGARCH)
{ x t = f ( t , x t − 1 , x t − 2 , … ) + ε t v t = h t e t ln h t = ω + ∑ i = 1 p η i ln h t − i + ∑ j = 1 q λ j g ( e t ) g ( e t ) = θ e t + γ [ ∣ e t ∣ − E ∣ e t ∣ ] = { ( θ + γ ) e t − γ E ∣ e t ∣ , e t > 0 ( θ − γ ) e t − γ E ∣ e t ∣ , e t < 0 \begin{cases} x_t=f(t,x_{t-1},x_{t-2},\dots)+\varepsilon_t\\\\ v_t=\sqrt{h_t}e_t\\\\ \ln{h_t}=\omega+\sum_{i=1}^p \eta_i\ln{h_{t-i}}+\sum_{j=1}^q\lambda_j g(e_t) \\\\ g(e_t)=\theta e_t+\gamma[|e_t|-E|e_t|]=\begin{cases}(\theta+\gamma)e_t-\gamma E|e_t|,\quad e_t>0\\(\theta-\gamma)e_t-\gamma E|e_t|,\quad e_t<0\end{cases} \end{cases} ⎩ ⎨ ⎧xt=f(t,xt−1,xt−2,…)+εtvt=htetlnht=ω+∑i=1pηilnht−i+∑j=1qλjg(et)g(et)=θet+γ[∣et∣−E∣et∣]={(θ+γ)et−γE∣et∣,et>0(θ−γ)et−γE∣et∣,et<0
模型的Matlab方法
残差异方差检验
Matlab中,可以使用archtest函数检验残差异方差,使用方法如下。
h=archtest(res)
h=archtest(res,Name,Value)
[h,pValue]=archtest(____)
[h,pValue,stst,eValue]=archtest(____)
res:用于检验的残差序列,丢失数据用NaN代替
h:测试的结果,长度等于测试的次数,h=1表示拒绝残差没有异方差的假设,h=0说明残差没有异方差性
pValue:检验统计量的概率p值
stat:检验统计量
cValue:检验统计量的关键值
建立ARCH以及GARAH模型
Matlab中,可以使用garch函数建立ARCH以及GARCH模型,使用方法如下。
Mdl=garch
Mdl=garch(P,Q)
Mdl=garch(Name,Value)
使用Matlab建立时间序列的ARCH模型步骤如下:
读取数据→数据可视化→平稳性检验→绘制自相关和偏自相关函数图→AIC准则给ARIMA模型定阶→用AIC准则定阶GARCH模型→建立模型
一个例子
给定一个时间序列,选择适当的模型拟合。
10.77 | 13.3 | 16.64 | 19.54 | 18.97 | 20.52 | 24.36 |
---|---|---|---|---|---|---|
23.51 | 27.16 | 30.8 | 31.84 | 31.63 | 32.68 | 34.9 |
33.85 | 33.09 | 35.46 | 35.32 | 39.94 | 37.47 | 35.24 |
33.03 | 32.67 | 35.2 | 32.36 | 32.34 | 38.45 | 38.17 |
32.14 | 39.7 | 49.42 | 47.86 | 48.34 | 62.5 | 63.56 |
67.61 | 64.59 | 66.17 | 67.5 | 76.12 | 79.31 | 78.85 |
81.34 | 87.06 | 86.41 | 93.2 | 82.95 | 72.96 | 61.1 |
61.27 | 71.58 | 88.34 | 98.7 | 97.31 | 97.17 | 91.17 |
80.2 | 85.12 | 81.4 | 70.87 | 57.75 | 52.35 | 67.5 |
87.95 | 85.46 | 84.55 | 98.16 | 102.42 | 113.02 | 119.95 |
122.37 | 126.96 | 122.79 | 127.96 | 139.2 | 141.05 | 140.87 |
137.08 | 145.53 | 145.59 | 134.36 | 122.54 | 106.92 | 97.23 |
110.39 | 132.4 | 152.3 | 154.91 | 152.69 | 162.67 | 160.31 |
142.57 | 146.54 | 153.83 | 141.81 | 157.83 | 161.79 | 142.07 |
139.43 | 140.92 | 154.61 | 172.33 | 191.78 | 199.27 | 197.57 |
189.29 | 181.49 | 166.84 | 154.28 | 150.12 | 165.17 | 170.32 |
先建立一个m文件,命名为data.m
data.m
function y=data()
data0=[10.77 13.3 16.64 19.54 18.97 20.52 24.36
23.51 27.16 30.8 31.84 31.63 32.68 34.9
33.85 33.09 35.46 35.32 39.94 37.47 35.24
33.03 32.67 35.2 32.36 32.34 38.45 38.17
32.14 39.7 49.42 47.86 48.34 62.5 63.56
67.61 64.59 66.17 67.5 76.12 79.31 78.85
81.34 87.06 86.41 93.2 82.95 72.96 61.1
61.27 71.58 88.34 98.7 97.31 97.17 91.17
80.2 85.12 81.4 70.87 57.75 52.35 67.5
87.95 85.46 84.55 98.16 102.42 113.02 119.95
122.37 126.96 122.79 127.96 139.2 141.05 140.87
137.08 145.53 145.59 134.36 122.54 106.92 97.23
110.39 132.4 152.3 154.91 152.69 162.67 160.31
142.57 146.54 153.83 141.81 157.83 161.79 142.07
139.43 140.92 154.61 172.33 191.78 199.27 197.57
189.29 181.49 166.84 154.28 150.12 165.17 170.32
];
data_new=reshape(data0',[1,112]);
y=data0;
首先,我们获取数据并且观察时间趋势
%% 获取数据并且观察时间趋势
data0=data();
data_new=reshape(data0',[112,1]);
figure(1)
plot(data_new)
我们利用pptest函数对时间序列进行平稳性检验,hp是逻辑值,如果hp=0,那么时间序列不平稳,需要进行差分使序列平稳化。
%% 平稳性检验
[hp,hpValue,stat,cValue,reg]=pptest(data_new,'model','TS')
if hp==0
data_diff=diff(data_new);
[hp,hpValue,stat,cValue,reg]=pptest(data_diff,'model','TS')
end
[hp,hpValue,stat,cValue,reg]=adftest(data_diff,'model','TS')
运行上述代码发现hp=0,那么时间序列不平稳,需要进行差分使序列平稳化。
接下来观察自相关和偏自相关函数图,观察是否存在拖尾和结尾性质
%% 绘制自相关和偏自相关函数图
figure(2) % 差分前的时间序列
subplot(2,1,1)
autocorr(data_new)
title('自相关图')
subplot(2,1,2)
autocorr(data_new)
title('偏自相关图')
parcorr(data_new)
figure(3) % 差分后的时间序列
subplot(2,1,1)
autocorr(data_diff)
title('差分后的自相关图')
subplot(2,1,2)
parcorr(data_diff)
title('差分后的偏自相关图')
由于未出现截尾,所以我们需要拟合ARIMA模型
这里有一点比较重要,对于一个时间序列,又可能存在某个i,j,使得这个模型拟合ARMA(i,j)会出现单位根情况,这会导致指数爆炸,进而报错,遇到这种情况,为使程序正常运行,我们需要跳过这样的i,j。
%% AIC定阶
maxLag=5;
AICSet=zeros(maxLag,maxLag);
AICSet(3,3)=800;
AICSet(4,4)=800;
AICSet(5,3)=800;
%本题中i=3 j=3,i=4 j=4,i=5 j=3时会出现单位根导致指数爆炸的情况,所以跳过这几种模型
parfor i=1:maxLag
for j=1:maxLag
if i==3&&j==3
continue
end
if i==4&&j==4
continue
end
if i==5&&j==3
continue
end
model_ARIMA=arima('ARLags',i,'MALags',j);
model_ARIMA=arima(i,1,j);
[EstMD1,EstParamCov,logL,info]=estimate(model_ARIMA,data_diff,'display','off')
AICSet(i,j)=aicbic(logL,length(info.X));
end
end
figure(4)
heatmap(AICSet/1000);
xlabel('MA Lags')
ylabel('AR Lags')
title('AIC准则图')
[ARLag,MALag]=find(AICSet==min(min(AICSet)));
title(['Optimal AR and MA Lags are(' num2str(ARLag) ',' num2str(MALag) ')'])
%Optimal AR and MA Lags are(4,2)
%这说明我们拟合的模型是ARIMA(4,1,2)
运行上述代码,可以得到下图
运行下面的代码,建立我们的ARIMA(4,1,2)模型,找出残差序列,检验残差序列是否存在异方差,即ARCH效应。
%% 建立模型
model_ARIMA=arima('ARLags',ARLag,'D',1,'MALags',MALag);
fit=estimate(model_ARIMA,data_new)
[res,~,~]=infer(fit,data_new);
[h,pVal]=archtest(res)
% pValue=0.0298,这说明存在ARCH效应。
我们根据选小原则,基于ARIMA(4,1,2)模型,建立GARCH模型
%% 用AIC准则定阶GARCH模型
% maxLag表示最大的延迟阶数,我们从GARCH(1,1)开始一直试探到GARCH(3,3),找出最小的AIC对应的模型
maxLag=3;
AICSet=nan(maxLag,maxLag);
for p=1:maxLag
for q=1:maxLag
Model_ARCH=arima('ARLags',ARLag,'D',1,'MALags',MALag,'Variance',garch(p,q));
% GARCH是ARIMA的一种扩展,ARLags表示自AR的阶数,MALags表示MA的阶数,D表示差分的阶数,Variance表示异方差,用garch(o,q)表示
[~,~,logL,info]=estimate(Model_ARCH,data_new,'display','off');
AICSet(p,q)=aicbic(logL,length(info.X));
end
end
figure;
heatmap(AICSet/1000);
xlabel('ARCH Lags')
ylabel('GARCH Lags')
title('Akaike information criteria')
[garchLags,archLags]=find(AICSet==min(min(AICSet)));
title(['Optimal GARCH and ARCH Lags are (' num2str(garchLags) ',' num2str(archLags) ')']);
% Optimal GARCH and ARCH Lags are (1,1)
% 这说明我们拟合的模型是ARCH(1,1)
根据AIC选小原则,我们得出模型为GARCH(1,1)
下面建立GARCH(1,1)模型
%% 建立GARCH模型
Model=arima('ARLags',ARLag,'D',1,'MALags',MALag,'Variance',garch(garchLags,archLags));
fit=estimate(Model,data_new)
以下部分是运行结果
%% MODEL
%ARIMA(4,1,2) Model (Gaussian Distribution):
% Value StandardError TStatistic PValue
% ________ _____________ __________ _________
% Constant 1.3678 0.44728 3.0581 0.0022274
% AR{4} -0.11418 0.10001 -1.1417 0.25357
% MA{2} -0.12001 0.1037 -1.1573 0.24715
%GARCH(1,1) Conditional Variance Model (Gaussian Distribution):
% Value StandardError TStatistic PValue
% _______ _____________ __________ __________
% Constant 1.3382 1.7759 0.75355 0.45112
% GARCH{1} 0.68573 0.12785 5.3634 8.1686e-08
% ARCH{1} 0.31427 0.15141 2.0756 0.03793
因此,我们的模型是
{
y
t
=
1.3678
−
0.11418
y
t
−
4
−
0.12001
ε
t
−
2
ε
t
=
v
t
h
t
v
t
∼
N
(
0
,
1
)
h
t
=
1.3382
+
0.68573
h
t
−
1
+
0.31427
ε
t
−
1
2
\begin{cases} y_t=1.3678-0.11418y_{t-4}-0.12001\varepsilon _{t-2}& \\ \varepsilon _t=v_t\sqrt{h_t}& v_t\sim N(0,1)\\ h_t=1.3382+0.68573h_{t-1}+0.31427\varepsilon _{t-1}^{2}& \\ \end{cases}
⎩
⎨
⎧yt=1.3678−0.11418yt−4−0.12001εt−2εt=vththt=1.3382+0.68573ht−1+0.31427εt−12vt∼N(0,1)