【转自半美人】matlab 实现 garch 模型波动率估计

代码高亮问题

数据获取

数据处理

描述性统计

时间序列平稳性检验

相关和偏自相关

arch效应检验

建立 garch 模型

波动率估计

代码及文档地址

代码高亮问题

这边首先说一个问题，你们看到的 matlab 代码没有高亮，看上去很不舒服。但这个锅是 csdn 的，真垃圾，连个语法高亮都搞不好，matlab 虽然近些年慢慢被 python 取代，但是用户数量还是挺多，不算小众语言。因为没有代码高亮可能，看得非常不爽，所以我特点制作了电子版，可以点击下载

数据获取

使用 tushare 获取历史股票数据。首先，需要 注册tushare账号 ，然后获取自己的 token。没有 tushare 账号的，可以点去这里注册个https://tushare.pro/register?reg=126259账号，好像直接点击链接会提示风险，需要复制打地址栏打开。然后还需要从 tushare 网站获取 matlab sdk 放到代码运行目录下，以供调用。

close all

clear

clc

%% 读取数据

stockcode = '600551.SH';

diary([stockcode,'.txt'])

addpath(genpath(pwd));

% 如果没有 tushare 的 token 按我说的去注册个账号，然后替换成自己的token

token = '***************2ad39c08b3d******************';

api = pro_api(token);

start_time = '20170101';

end_time = '20180101';

ktype = 'D';

data = pro_bar(stockcode, api, start_time, end_time,ktype,'E','qfq');

数据处理

在获取了数据之后，我们需要对数据进行一些处理。

% 翻转数据

data = flipud(data);

temp = data.trade_date;

temp = char(temp);

temp = str2num(temp);

tradedate = datetime(temp,'ConvertFrom','yyyymmdd','format','yyyy-MM-dd');

%% 计算对数收益率

ret = price2ret(data.close);

% 绘制对数收益率图

x = tradedate(1:end-1);

figure;

plot(x,ret);

title([stockcode,'日收益率'])

saveas(gcf,[stockcode,'日收益率图.jpg']);

描述性统计

然后，对收益率数据进行描述性统计：

% 可以看出波动聚集

disp([stockcode,' 波动率预测']);

% 描述性统计

mu = mean(ret);

ret_median = median(ret);

ret_max = max(ret);

ret_min = min(ret);

sigma = std(ret);

ret_skew = skewness(ret);

ret_kurt = kurtosis(ret);

% h为测试结果,若h=0,则可以认为X是服从正态分布的;若h=1,则可以否定X服从正态分布;

% p为接受假设的概率值,P越接近于0,则可以拒绝是正态分布的原假设;

% jbstat为测试统计量的值;

% cv为是否拒绝原假设的临界值.

[h,p,jbstat,cv]=jbtest(ret);

disp(['收益序列均值为',num2str(mu)]);

disp(['收益序列中位数为',num2str(ret_median)]);

disp(['收益序列最大值为',num2str(ret_max)]);

disp(['收益序列最小值为',num2str(ret_min)]);

disp(['收益序列标准差为',num2str(ret_max)]);

disp(['收益序列偏度为',num2str(ret_skew)]);

disp(['收益序列峰度为',num2str(ret_kurt)]);

disp(['收益序列jb统计量为',num2str(jbstat)]);

相应的命令行结果输出：

600551.SH 波动率预测

警告: P is less than the smallest tabulated value, returning 0.001.

> In jbtest (line 136)

In garch_var (line 42)

收益序列均值为-0.0030367

收益序列中位数为-0.00062952

收益序列最大值为0.032405

收益序列最小值为-0.10542

收益序列标准差为0.032405

收益序列偏度为-2.1069

收益序列峰度为11.7364

收益序列jb统计量为572.3164

时间序列平稳性检验

我们利用 adf 检验对时间序列平稳性进行检验：

%% adf单位根检验

alpha=0.1;

disp('%%%%%%%%%%%%%%%%%%%%%%');

[h,pValue,stat,cValue] = adftest(ret,'alpha',alpha);

fprintf(['在%d%%的置信水平下，单位根检验统计量为%f,\n',...

'临界值为%f，概率为%f\n'],alpha*100,stat,cValue,pValue);

alpha=0.05;

disp('%%%%%%%%%%%%%%%%%%%%%%');

[h,pValue,stat,cValue] = adftest(ret,'alpha',alpha);

fprintf(['在%d%%的置信水平下，单位根检验统计量为%f,\n',...

'临界值为%f，概率为%f\n'],alpha*100,stat,cValue,pValue);

alpha=0.01;

disp('%%%%%%%%%%%%%%%%%%%%%%');

[h,pValue,stat,cValue] = adftest(ret,'alpha',alpha);

fprintf(['在%d%%的置信水平下，单位根检验统计量为%f,\n',...

'临界值为%f，概率为%f\n'],alpha*100,stat,cValue,pValue);

命令行结果输出：

%%%%%%%%%%%%%%%%%%%%%%

在10%的置信水平下，单位根检验统计量为-14.355258,

临界值为-1.615200，概率为0.001000

%%%%%%%%%%%%%%%%%%%%%%

在5%的置信水平下，单位根检验统计量为-14.355258,

临界值为-1.942690，概率为0.001000

%%%%%%%%%%%%%%%%%%%%%%

在1%的置信水平下，单位根检验统计量为-14.355258,

临界值为-2.581820，概率为0.001000

%%%%%%%%%%%%%%%%%%%%%%

相关和偏自相关

我们判断相关和偏自相关关系，主要是想初略判断后面 arch 效应检验时所要用到的滞后阶数：

% % 通过残差平方的偏自相关函数判断arch模型阶数

figure;

plot(x,at.^2);

name = [stockcode,' 收益率残差平方序列图'];

title(name);

saveas(gcf,[name,'.jpg']);

figure;

subplot(2,1,1);

autocorr(at.^2);

name = [stockcode,' 收益率残差平方序列自相关函数图'];

title(name);

subplot(2,1,2);

parcorr(at.^2);

name = [stockcode,' 收益率残差平方序列偏自相关函数图'];

title(name);

name = [stockcode,' 收益率残差平方序列自相关和偏自相关函数图'];

saveas(gcf,[name,'.jpg']);

至于为什么可以用偏自相关函数给arch模型定阶，可以看蔡瑞胸的金融时间序列分析。

arch效应检验

我们，取滞后阶数为3，用拉格朗日乘子法进行检验。

% 拉格郎日乘子检验

lags = 3;

alpha=0.05;

disp('%%%%%%%%%%%%%%%%%%%%%%');

disp('arch效应拉格朗日乘子检验');

[h,pValue,stat,cValue] = archtest(at,'Lags',lags,'Alpha',alpha);

fprintf(['在%d%%的置信水平下，滞后阶数为%d时,检验统计量为%f,\n',...

'临界值为%f，概率为%f\n'],alpha*100,lags,stat,cValue,pValue);

1

2

3

4

5

6

7

8

命令行结果输出：

%%%%%%%%%%%%%%%%%%%%%%

arch效应拉格朗日乘子检验

在5%的置信水平下，滞后阶数为3时,检验统计量为27.885162,

临界值为7.814728，概率为0.000004

%%%%%%%%%%%%%%%%%%%%%%

garch模型参数估计

这样，我们确认存在非常显著的arch效应。

建立 garch 模型

我们，此处建立 garch(1，1)模型：

%% garch模型建立

disp('%%%%%%%%%%%%%%%%%%%%%%');

md = garch(1,1);

disp('garch模型参数估计');

estMd = md.estimate(ret,'Display','off');

summarize(estMd);

命令行结果输出：

%%%%%%%%%%%%%%%%%%%%%%

garch模型参数估计

GARCH(1,1) Conditional Variance Model (Gaussian Distribution)

Effective Sample Size: 146

Number of Estimated Parameters: 3

LogLikelihood: 387.87

AIC: -769.74

BIC: -760.79

Value       StandardError    TStatistic      PValue

__________    _____________    __________    __________

Constant    2.8435e-05     2.8106e-05        1.0117         0.31168

GARCH{1}       0.84958        0.14246        5.9639      2.4636e-09

ARCH{1}       0.019211       0.016655        1.1535         0.24872

波动率估计

我们将 garch 模型无条件方差作为估计的波动率：

%% 无条件方差即波动率计算

sigma = estMd.Constant/(1-estMd.GARCH{1}-estMd.ARCH{1});

% 年化

y = 250;

sigma_y = sqrt(sigma*y);

fprintf('假定一年为%d天，年化波动率为%f\n',y,sigma_y);

diary off

1

2

3

4

5

6

7

命令行结果：

假定一年为250天，年化波动率为0.232768

1

代码及文档地址

————————————————