【转自半美人】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
代码及文档地址
————————————————
|
|
转播
