ARCH模型以及编程实现

Logistic..

已于 2022-07-15 00:06:59 修改

阅读量4.1k

点赞数 11

文章标签：概率论机器学习 python

于 2022-07-15 00:05:12 首次发布

本文链接：https://blog.csdn.net/C1172440795/article/details/125794667

版权

文章目录

ARCH模型以及编程实现

ARCH模型以及编程实现

ARCH模型的引入

传统的ARIMA模型拟合非平稳时间序列，残差序列 $\{\varepsilon_t\}$ 通常满足下面的三个假定条件

零均值
$E(\varepsilon_t)=0$
纯随机
$Cov(\varepsilon_t,\varepsilon_{t-i})=0,\forall i\ge1$
方差齐性
$Var(\varepsilon_t)=\sigma_t^2$

传统的ARIMA模型假定时间序列的长期方差是相同的，但是，通常情况下，时间序列的存在异方差的，即
$Var(\varepsilon_t)=h(t)$
根据以上条件，我们可以得出，残差序列的方差就是它平方的期望，即
$Var(\varepsilon_t)=E(\varepsilon_t^2)$

我们可以根据残差平方图 $\varepsilon_t$ 关于 $t$ 变化的二维坐标图，对方差齐性直观判断。

我们对存在异方差的时间序列有两种处理方式

1.如果已知异方差函数具体形式，进行方差齐性变化

2.如果不知异方差函数具体形式，拟合条件异方差模型，即ARCH

异方差性检验

Q检验

原假设 $H_0:$ 残差平方序列纯随机，即 $\rho_1=\rho_2=\dots=\rho_q=0$
备择假设 $H_1:$ 残差平方序列具有自相关性，即 $\rho_1,\rho_2,\dots,\rho_q$ 不全为 $0$

检验统计量
$Q(q)=n(n+2)\sum_{i=1}^{q}\frac{\rho_i^2}{n-i}$

其中， $n$ 为观察序列长度， $\rho_i$ 为残差序列延迟i的自相关系数，设

$\hat{\sigma}^2=\frac{\sum_{i=1}^n\varepsilon_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}}$

Q统计量近似服从于自由度为 $q - 1$ 的 $\chi^2$ 分布，即 $Q(q)\sim\chi^2(q-1)$

若 $Q(q)\ge\chi_{1-\alpha}^2(q-1)$ ，则拒绝原假设，认为存在异方差
若 $Q(q)\lt\chi_{1-\alpha}^2(q-1)$ ，则接受原假设，认为不存在异方差

LM检验

原假设 $H_0:$ 残差平方序列纯随机，即 $\rho_1=\rho_2=\dots=\rho_q=0$
备择假设 $H_1:$ 残差平方序列具有自相关性，即 $\rho_1,\rho_2,\dots,\rho_q$ 不全为0

若 $H_0$ 为真，则 $LM(q)=W^TW$
$W=(\frac{\rho_1^2}{\sigma^2},\frac{\rho_2^2}{\sigma^2},\dots,\frac{\rho_q^2}{\sigma^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}}, \hat{\sigma}^2=\frac{\sum_{i=1}^n\varepsilon_i^2}{n}$

$L M (q)$ 统计量近似服从于自由度为 $q - 1$ 的 $\chi^2$ 分布，即 $LM(q)\sim\chi^2(q-1)$

若 $LM(q)\ge\chi_{1-\alpha}^2(q-1)$ ，则拒绝原假设，认为存在异方差
若 $LM(q)\lt\chi_{1-\alpha}^2(q-1)$ ，则接受原假设，认为不存在异方差

$A RC H (q)$ 模型

具有
$\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}$
结构的模型称为q阶自回归条件异方差模型，简记为ARCH(q).

$G A RC H (p, q)$ 模型

具有
$\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}$
的模型称为广义自回归条件异方差模型，简记为 $G A RC H (p, q)$ .

$A R (m) - G A RC H (p, q)$ 模型

$\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}$

GARCH的衍生模型

指数GARCH模型(EGARCH)

$\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}$

模型的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

因此，我们的模型是
$\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}$