人大经济论坛 › 标签 › MATLAB程序

标签: MATLAB程序经管大学堂：名校名师名课

相关帖子	版块	作者	回复/查看	最后发表

【海量分享2】Matlab 2008-2018全套安装软件（2018a 2017b 2017a 等）+视频+书籍+代码	MATLAB等数学软件专版	sdxzfb 2018-2-18	895 48275	zcj111 2024-2-24 15:45:15
matlab教程[从入门到精通]程序	MATLAB等数学软件专版	瓶子02 2019-4-14	27 3430	543913616 2023-12-20 08:39:09
广义帕累托分布的参数估计-极大似然估计 - [悬赏 100 个论坛币]	MATLAB等数学软件专版	yubuxiao 2019-2-16	3 2780	hhdxlll 2022-10-26 16:04:22
求实物期权相关matlab代码！	MATLAB等数学软件专版	miechen1979 2021-11-28	0 343	miechen1979 2021-11-28 11:36:53
随机波动模型 MATLAB00 - [悬赏 100 个论坛币]	MATLAB等数学软件专版	Kingbord 2016-10-6	3 1787	扭曲的虫洞 2021-1-2 17:24:31
MATLAB之入门教程（视频+数据+代码+源文件）	MATLAB等数学软件专版	瓶子02 2019-6-1	14 2482	zhaochuan2 2020-6-25 11:08:34
求教用contour画下面函数的二元等高线图。	MATLAB等数学软件专版	cyy1118 2020-3-7	1 918	402play 2020-3-11 18:29:03
elhost的LMsarsem_panal错误警告	MATLAB等数学软件专版	果果1106 2014-12-29	11 5360	qingqing1993 2018-5-24 11:24:27
求教这两种演化仿真图像的程序，感谢！	MATLAB等数学软件专版	Symphony.Huo1 2017-3-21	0 748	Symphony.Huo1 2017-3-21 23:35:28
MATLAB的一些算法程序	MATLAB等数学软件专版	15951005885 2016-7-14	2 952	15951005885 2016-7-20 19:42:11
个人写的一些Matlab程序集	量化投资	wuyy2007 2015-6-2	89 9260	luotangsang 2016-1-12 16:15:25
Matlab 2012b程序	数据分析与数据挖掘	jamestoneyin 2013-8-24	29 5461	stepward 2013-10-10 09:53:40
MATLAB有用的工具箱	MATLAB等数学软件专版	hssunduan 2010-4-11	2 1626	wangtian0828 2010-4-17 18:46:23

更多...

相关日志

分享中国股市波动率的双重不对称性及其解释 Matlab程序: stanleyjunjun 2013-9-20 13:37; 有很多同学来信要求提供这篇文章以及此前几篇文章的Matlab程序，这里公开《金融研究》上发表论文的程序，这个程序比较成熟。不过程序也还要使用者自己调试并完善，本文程序仅供参考，不确保一定正确，如果使用本文程序所有的风险请自担。 stanleyjunjun@163.com 发表于：中国股市波动率的双重不对称性及其解释——基于MS-TGARCH模型的MCMC估计和分析，金融研究 2011.03，第138-152页。中国股市波动率的双重不对称性及其解释 —— 基于 MS-TGARCH 模型的 MCMC 估计和分析朱钧钧谢识予复旦大学经济学院，上海 200433 摘要：本文应用 MCMC 方法估计了上证综指的 MS-TGARCH 模型，并得出中国股市的波动率存在双重不对称性。其一，每个波动状态中，中国股市的波动率都存在不对称性，其中高波动状态的波动率对正向冲击的反应显著大于负向冲击，而低波动状态则刚好相反；其二，不同状态之间，中国股市波动率的不对称性刚好相反，并且高波动状态伴随着显著大于 0 的平均收益率，而低波动状态的平均收益率显著小于 0 。这种双重不对称性显著区别于国外股市。至于中国股市上述独特特征的形成原因，本文从投资者心理和行为角度，联系中国股市的交易制度设计和做空机制的缺乏进行了解释。关键词：马尔可夫转换模型股市波动率 MCMC 方法 JEL 分类号：文献标识码： A 作者简介：朱钧钧，复旦大学经济学院博士生；谢识予，复旦大学经济学院教授，博士生导师。作者感谢复旦大学（教育部）金融创新研究生开放实验室创新项目基金资助，复旦大学研究生创新基金资助。本文的第一版在厦门的第 10 届中国数量经济学年会上交流过，感谢与会专家的点评，文责自负。 %------------------------------------------------------------------------------ % Main Code %------------------------------------------------------------------------------ clc data = load('SEIweek200910.txt'); y = log(data(312:end,7))-log(data(311:end-1,7)); T = size(y,1); x = ones(T,1); sn = 2; % number of state nk = 1; % regressor number in mean equation np = 1; % p-term in GARCH % prior for parameters pi_pri = ; % prior for pi distrib = 'tdistrib'; indexSwitch =0; m=100; % grid points % choose starting values % P= ; %spec = garchset('distribution', 'T','p',0,'q',3); % = garchfit(spec,y); %garchdisp(coeff,errors) % coeff.C = ; coeff.K = 0.00015; coeff.ARCH = % default model: y(t) = C +e(t); h(t) = K + GARCH*h(t-1) + ARCH*e(t-1)^2 % first try, after checking for GARCH % beta = ; % a0 = ; % % a1 = ; % a2 = ; % after 1000 simulation beta = ; a0 = ; % a1p = ; a1n = ; b1 = ; v = 5; var_err = y.^2; beta_1min = -0.002; beta_1max = 0.022; beta_1step = (beta_1max - beta_1min)/m; beta_2min = -0.02; beta_2max = 0.008; beta_2step = (beta_2max - beta_2min)/m; a0_1min = 8e-06; a0_1max = 3e-04; a0_1step = (a0_1max - a0_1min)/m; a0_2min = 4e-06; a0_2max = 4e-04; a0_2step = (a0_2max - a0_2min)/m; a1p_1min = 0.01; a1p_1max = 0.8; a1p_1step = (a1p_1max - a1p_1min)/m; a1p_2min = 0.01; a1p_2max = 0.6; a1p_2step = (a1p_2max - a1p_2min)/m; a1n_1min = 0.01; a1n_1max = 0.6; a1n_1step = (a1n_1max - a1n_1min)/m; a1n_2min = 0.01; a1n_2max = 0.6; a1n_2step = (a1n_2max - a1n_2min)/m; b1_1min = 0.3; b1_1max = 0.95; b1_1step = (b1_1max - b1_1min)/m; b1_2min = 0.4; b1_2max = 0.95; b1_2step = (b1_2max - b1_2min)/m; v_min =2.3; v_max = 20; v_step = (v_max - v_min)/m; % save result as beta_ = ; a1p_ = ; b1_ = ; P_ = ; prob = zeros(T,1); S_ = ; % switching number, sm: (T-1)*2 tm=zeros(sn,sn); % total number for im=1:sn for jm=1:sn tm(im,jm)=size(sm(sm(:,1)==im sm(:,2)==jm,:),1); % total number of each transition from i state to j state end rr = pi_pri(im,:)+tm(im,:); P_draw = chi2rnd(2*rr); P(im,:)=P_draw./repmat(sum(P_draw,2),1,sn); % P(im,:)=draw_dirichlet(1,sn,pi_pri(im,:)+tm(im,:),0); % rm: prior for number of state transition (ij) end %--------------------------------------------------------------------- %Specify the number of replications nburn=0; %number of burnin replications nrep=30000; %number of retained replications ntotal=nburn+nrep; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Now start Gibbs loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i = 1:ntotal completion_process=i/ntotal %------------------------------------------------------- % Griddy Gibbs a0(1) griddy = 'a0_1'; param = ; % for case of t distribution; plot(a_d, like) = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, a0_1min, a0_1step); a0(1) = griddy_out; % a0(1) = max_out; %------------------------------------------------------- % Griddy Gibbs a0(2) griddy = 'a0_2'; param = ; = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, a0_2min, a0_2step); a0(2) = griddy_out; % a0(2) = max_out; %----------------------------------------------------------------------- % Griddy Gibbs a1(1) griddy = 'a1p_1'; param = ; = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, a1p_1min, a1p_1step); a1p(1) = griddy_out; %----------------------------------------------------------------------- % Griddy Gibbs a1(2) griddy = 'a1p_2'; param = ; = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, a1p_2min, a1p_2step); a1p(2) = griddy_out; %----------------------------------------------------------------------- % Griddy Gibbs a1(1) griddy = 'a1n_1'; param = ; = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, a1n_1min, a1n_1step); a1n(1) = griddy_out; %----------------------------------------------------------------------- % Griddy Gibbs a1(2) griddy = 'a1n_2'; param = ; = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, a1n_2min, a1n_2step); a1n(2) = griddy_out; %----------------------------------------------------------------------- % Griddy Gibbs a2(1) griddy = 'b1_1'; param = ; = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, b1_1min, b1_1step); b1(1) = griddy_out; %----------------------------------------------------------------------- % Griddy Gibbs a2(2) griddy = 'b1_2'; param = ; = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, b1_2min, b1_2step); b1(2) = griddy_out; %------------------------------------------------------ % Griddy Gibbs Sampling beta(1) griddy = 'beta_1'; param = ; = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, beta_1min, beta_1step); beta(1) = griddy_out; % beta(1) = max_out; %------------------------------------------------------ % Griddy Gibbs Sampling beta(2) griddy = 'beta_2'; param = ; = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, beta_2min, beta_2step); beta(2) = griddy_out; griddy = 'v'; param = ; = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, v_min, v_step); v = griddy_out; %-------------------------------------------------------------------- %draw S conditional on rest %-------------------------------------------------------------------- tau = 1; dif_like = MSGARCH_3_S(y, x, S, param, distrib, tau+1); post_S1 = P(1,S(tau+1)); post_S2 = P(2,S(tau+1)); post_S1N = post_S1*dif_like/(post_S1*dif_like + post_S2); %post_S2N = post_S2/(post_S1+post_S2); random_S = rand; if random_S = post_S1N S(1) = 1; else S(1) = 2; end prob(1) = post_S1N; for tau = 2:T-1 dif_like = MSGARCH_3_S(y, x, S, param, distrib, tau); post_S1 = P(S(tau-1),1)*P(1,S(tau+1)); post_S2 = P(S(tau-1),2)*P(2,S(tau+1)); post_S1N = post_S1*dif_like/(post_S1*dif_like + post_S2); % post_S2N = post_S2/(post_S1+post_S2); random_S = rand; if random_S = post_S1N S(tau) = 1; else S(tau) = 2; end prob(tau) = post_S1N; end tau = T; dif_like = MSGARCH_3_S(y, x, S, param, distrib, tau); post_S1 = P(S(tau-1),1); post_S2 = P(S(tau-1),2); post_S1N = post_S1*dif_like/(post_S1*dif_like + post_S2); % post_S2N = post_S2/(post_S1+post_S2); random_S = rand; if random_S = post_S1N S(T) = 1; else S(T)=2; end prob(tau) = post_S1N; % for tau = 2:T-1; = S_Post_mixS2b(y, x, beta, S, a0, a1, a2, tau+1, lamda); end % plot(dif_like_tau) %-------------------------------------------------------------------- %draw P conditional on S %-------------------------------------------------------------------- % P=zeros(sn,sn); sm= ; % switching number, sm: (T-1)*2 tm=zeros(sn,sn); % total number for im=1:sn for jm=1:sn tm(im,jm)=size(sm(sm(:,1)==im sm(:,2)==jm,:),1); % total number of each transition from i state to j state end rr = pi_pri(im,:)+tm(im,:); P_draw = chi2rnd(2*rr); P(im,:)=P_draw./repmat(sum(P_draw,2),1,sn); % P(im,:)=draw_dirichlet(1,sn,pi_pri(im,:)+tm(im,:),0); % rm: prior for number of state transition (ij) end %------------------------------------------------------- %identification constraint %------------------------------------------------------- % desired result: a_1 a_2 a_3 indexconstraint = beta(1) beta(2); % to fix state 1 as low violatility state % indexconstraint = P(1,1)P(2,2); % not good % indexconstraint = a0(1) a0(2); % sounds good % indexconstraint = b(1)b(2); if indexconstraint indexSwitch=indexSwitch+1; ndxv= ; beta=beta(ndxv); a0 = a0(ndxv); a1p = a1p(ndxv); a1n = a1n(ndxv); b1 = b1(ndxv); P=P(ndxv,ndxv); S=3*ones(T,1) - S; end %----------------------------------------------------- % save accepted draw after burn-in replications %----------------------------------------------------- if inburn beta_ = ; a0_ = ; a1p_ = ; a1n_ = ; b1_ = ; v_ = ; Pvec = ; P_ = ; prob_ = ; S_ = ; end end MS9.P=P; MS9.P_=P_; MS9.S=S; MS9.beta_=beta_; MS9.a0_=a0_; MS9.a1p_=a1p_; MS9.a1n_=a1n_; MS9.b1_=b1_; MS9.v_ = v_; MS9.S_ = S_; % disp('indexSwitch') indexSwitch % disp('nswitch1') figure(1) subplot(321); plot(beta_) title ('convergence of beta'); subplot(322); plot (a0_); title('convergence of a0 term') subplot(323); plot(a1p_) title('convergence of ARCH term') subplot(324); plot(a1n_) title('convergence of ARCH term') subplot(325); plot(b1_) title('convergence of GARCH term') subplot(326);plot(v_); title('convergence or v') proc = 0.2; disp('mean and std of beta') for tt = 1:2 end disp('mean and std of volatility param') for tt = 1:2 end %------------------------------------------------------------------------ % function sum_log_like %------------------------------------------------------------------------ function sum_log_like = MSGARCH_3_Post(y, x, S, param, distrib) T = size(y,1); switch distrib case 'normal' beta = param(1,:)'; a0 = param(2,:)'; a1p = param(3,:)'; a1n = param(4,:)'; b1 = param(5,:)'; case 'tdistrib' beta = param(1,:)'; a0 = param(2,:)'; a1p = param(3,:)'; a1n = param(4,:)'; b1 = param(5,:)'; v = param(6,1); end e = y - x.*beta(S); h = zeros(T,1); % h(1) = 0.5*sum(a0)/(1-0.5*sum(a1)-0.5*sum(a2)); % can produce complex value if S(1) == 1 h(1) = 3*a0(1); else h(1) = 3*a0(2); end for t = 2:T h(t) = a0(S(t)) + a1p(S(t))*(e(t-1)0)*(e(t-1)^2) + a1n(S(t))*(e(t-1)0)*(e(t-1)^2) + b1(S(t))*h(t-1); end switch distrib case 'normal' log_like = - 0.5*log(2*pi) - 0.5*log(h) - 0.5*(e.^2)./h; case 'tdistrib' log_like = log(gamma(0.5*(v+1))) - log(gamma(v/2)) - 0.5*log(v*pi*h) - (v+1)/2*log(1+(e.^2)./(v*h)); end sum_log_like = sum(log_like); %---------------------------------------------------------------------------- % function MSGARCH_3_Griddy %---------------------------------------------------------------------------- function = MSGARCH_3_Griddy(y, x, S, param, distrib, griddy, m, a_min, a_step) a_d = zeros(m,1); sum_log_like = zeros(m,1); switch distrib case 'normal' beta = param(1,:)'; a0 = param(2,:)'; a1p = param(3,:)'; a1n = param(4,:)'; b1 = param(5,:)'; case 'tdistrib' beta = param(1,:)'; a0 = param(2,:)'; a1p = param(3,:)'; a1n = param(4,:)'; b1 = param(5,:)'; v = param(6,1); end for i = 1:m a_draw = a_min + i*a_step; switch griddy case 'beta_1' beta(1) = a_draw; case 'beta_2' beta(2) = a_draw; case 'a0_1' a0(1) = a_draw; case 'a0_2' a0(2) = a_draw; case 'a1p_1' a1p(1) = a_draw; case 'a1n_1' a1n(1) = a_draw; case 'a1p_2' a1p(2) = a_draw; case 'a1n_2' a1n(2) = a_draw; case 'b1_1' b1(1) = a_draw; case 'b1_2' b1(2) = a_draw; case 'v' v = a_draw; end a_d(i) = a_draw; switch distrib case 'normal' param = ; sum_log_like(i) = MSGARCH_3_Post(y, x, S, param, distrib); case 'tdistrib' param = ; sum_log_like(i) = MSGARCH_3_Post(y, x, S, param, distrib); end end like = exp(sum_log_like-max(sum_log_like)); a_cdf = zeros(m,1); a_cdf(1) = like(1); for i = 2:m; a_cdf(i) = a_cdf(i-1) + like(i); end; a_cdfn = a_cdf./a_cdf(m); a_c = rand; = min(abs(a_c - a_cdfn)); if index == 1 | index == m; griddy_out = a_d(index); else if a_c a_cdfn(index) | a_c == a_cdfn(index) griddy_out = (a_d(index)-a_d(index-1))*(a_c - a_cdfn(index-1))/(a_cdfn(index) - a_cdfn(index-1)) + a_d(index-1); else index = index +1; griddy_out = (a_d(index)-a_d(index-1))*(a_c - a_cdfn(index-1))/(a_cdfn(index) - a_cdfn(index-1)) + a_d(index-1); end end %----------------------------------------------------------------------- % function MSGARCH_3_S %----------------------------------------------------------------------- function = MSGARCH_3_S(y, x, S, param, distrib, tau) T = size(y,1); switch distrib case 'normal' beta = param(1,:)'; a0 = param(2,:)'; a1p = param(3,:)'; a1n = param(4,:)'; b1 = param(5,:)'; case 'tdistrib' beta = param(1,:)'; a0 = param(2,:)'; a1p = param(3,:)'; a1n = param(4,:)'; b1 = param(5,:)'; v = param(6,1); end S(tau) = 1; e1 = y - x.*beta(S); h1 = zeros(T,1); if S(1) == 1 h1(1) = 3*a0(1); else h1(1) = 3*a0(2); end % switch S_tau % not work well % case 'S_1tau' % % h(1) = a0(1)/(1-a1(1)-a2(1)); % in case of negative (1-a1-a2) % h(1) = 2*a0(1); % case 'S_2tau' % h(1) = a0(2)/(1-a1(2)-a2(2)); % h(1) = 2*a0(2); % end for t = 2:T h1(t) = a0(S(t)) + a1p(S(t))*(e1(t-1)0)*(e1(t-1)^2) + a1n(S(t))*(e1(t-1)0)*(e1(t-1)^2) + b1(S(t))*h1(t-1); end switch distrib case 'normal' log_like1 = - 0.5*log(2*pi) -0.5*log(h1(tau:T)) - 0.5*(e1(tau:T).^2)./h1(tau:T); sum_log_like1 = sum(log_like1); log_like_tau1 = - 0.5*log(2*pi) -0.5*log(h1(tau)) - 0.5*(e1(tau).^2)./h1(tau); case 'tdistrib' log_like1 = log(gamma(0.5*(v+1))) - log(gamma(v/2)) - 0.5*log(v*pi*h1(tau:T)) - (v+1)/2*log(1+(e1(tau:T).^2)./(v*h1(tau:T))); sum_log_like1 = sum(log_like1); log_like_tau1 = log(gamma(0.5*(v+1))) - log(gamma(v/2)) - 0.5*log(v*pi*h1(tau)) - (v+1)/2*log(1+(e1(tau).^2)./(v*h1(tau))); end S(tau) = 2; e2 = y - x.*beta(S); h2 = zeros(T,1); if S(1) == 1 h2(1) = 3*a0(1); else h2(1) = 3*a0(2); end for t = 2:T h2(t) = a0(S(t)) + a1p(S(t))*(e2(t-1)0)*(e2(t-1)^2) + a1n(S(t))*(e2(t-1)0)*(e2(t-1)^2) + b1(S(t))*h2(t-1); end switch distrib case 'normal' log_like2 = - 0.5*log(2*pi) -0.5*log(h2(tau:T)) - 0.5*(e2(tau:T).^2)./h2(tau:T); sum_log_like2 = sum(log_like2); log_like_tau2 = - 0.5*log(2*pi) -0.5*log(h2(tau)) - 0.5*(e2(tau).^2)./h2(tau); case 'tdistrib' log_like2 = log(gamma(0.5*(v+1))) - log(gamma(v/2)) - 0.5*log(v*pi*h2(tau:T)) - (v+1)/2*log(1+(e2(tau:T).^2)./(v*h2(tau:T))); sum_log_like2 = sum(log_like2); log_like_tau2 = log(gamma(0.5*(v+1))) - log(gamma(v/2)) - 0.5*log(v*pi*h2(tau)) - (v+1)/2*log(1+(e2(tau).^2)./(v*h2(tau))); end dif_log_like = sum_log_like1 - sum_log_like2; dif_like = exp(dif_log_like); dif_like_tau = exp(log_like_tau1 - log_like_tau2); %------------------------------------------------------------------------------------ % function MSGARCH_3_PlotSurf %------------------------------------------------------------------------------------ function = MSGARCH_3_PlotSurf(y,x,S,param,griddy, distrib,m,hor_min, hor_step, ver_min, ver_step) T = size(y,1); switch distrib case 'normal' beta = param(1,:)'; a0 = param(2,:)'; a1p = param(3,:)'; a1n = param(4,:)'; b1 = param(5,:)'; case 'tdistrib' beta = param(1,:)'; a0 = param(2,:)'; a1p = param(3,:)'; a1n = param(4,:)'; b1 = param(5,:)'; v = param(6,1); end sum_log_like = zeros(m,m); for i = 1:m ver_draw = ver_min + (i-1)*ver_step; for j = 1:m hor_draw = hor_min + (j-1)*hor_step; switch griddy case 'beta' beta(1) = hor_draw; beta(2) = ver_draw; case 'a0' a0(1) = hor_draw; a0(2) = ver_draw; case 'a1p' a1p(1) = hor_draw; a1p(2) = ver_draw; case 'a1n' a1n(1) = hor_draw; a1n(2) = ver_draw; case 'b1' b1(1) = hor_draw; b1(2) = ver_draw; end e = y - x.*beta(S); h = zeros(T,1); % h(1) = 0.5*sum(a0)/(1-0.5*sum(a1)-0.5*sum(a2)); % can produce complex value if S(1) == 1 h(1) = 3*a0(1); else h(1) = 3*a0(2); end for t = 2:T h(t) = a0(S(t)) + a1p(S(t))*(e(t-1)0)*(e(t-1)^2) + a1n(S(t))*(e(t-1)0)*(e(t-1)^2) + b1(S(t))*h(t-1); end switch distrib case 'normal' log_like = - 0.5*log(2*pi) - 0.5*log(h) - 0.5*(e.^2)./h; case 'tdistrib' log_like = log(gamma(0.5*(v+1))) - log(gamma(v/2)) - 0.5*log(v*pi*h) - (v+1)/2*log(1+(e.^2)./(v*h)); end sum_log_like(i,j) = sum(log_like); end end ver = zeros(m,1); hor = zeros(m,1); for ii = 1:m ver(ii) = ver_min + (ii-1)*ver_step; hor(ii) = hor_min + (ii-1)*hor_step; end like = exp(sum_log_like - max(max(sum_log_like))*ones(m,m));; 个人分类: Bayesian Econometrics|31 次阅读|1 个评论