MRS DCC GARCH 模型的MATLAB 程序修改

8# zhangtao 朱钧钧的没有错误的。也许我上传的时候，删掉了 Lesage 的一个函数，那个 norm_rnd()，这个函数是用来做正态分布的随机数。因为我装了 Lesage  的计量经济学程序文件夹，那个 JPL7,所以爱随便改程序的我删掉了那个。其他朋友有的说有问题，有的说就是看到数字乱跳，然后19分钟后，有结果。朱钧钧的程序没问题的。其实除了我们在做的MRS DCC MVGARCH  。 KEVEN 的工具箱里有 DCC MVGARCH  的。调用都差不多。你问的那个DEMO 里的 MVGARCH ，也是在 USCD_GARCH TOOLBOX 里。这个人是 ENGEL 的学生，在牛津大学教学吧。他还有另一个工具箱，那个还涉及比较前沿的 “已实现波动模型（RV）”，有程序代码了。其实ENGEL 的 ARCH 是1980s推出的，后来有几个学生，一个是 波洛索夫，这个推出了 GARCH ，他的另一个学生 KEVEN 推出了 DCC GARCH 。后来，engel还将注意力放在另外的前沿领域，这个你可以看看经典的书籍，汉密尔顿的《时间序列分析》。据说这本书在国外，以英镑计价的话，会比较贵。因为中国是发展中国家嘛。所以很多国外的书籍都比国外的便宜。

老兄,刚看了marcucci和你的 source code
发现有一个问题,值得思考与讨论
你的ms_dcc_mvgarch.m是
for i=1:k
    fprintf(1,'Estimating GARCH model for Series %d\n',i)
    [univariate{i}.parameters, univariate{i}.likelihood, univariate{i}.stderrors, univariate{i}.robustSE, univariate{i}.ht, univariate{i}.scores] ...
        = fattailed_garch(data(:,i) , archP(i) , garchQ(i) , 'NORMAL',[], options);
    stdresid(:,i)=data(:,i)./sqrt(univariate{i}.ht);
end

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

还是应该像marcucci,mrsgarchestfor_con_all.m
[parameters, LLF, EXITFLAG, OUTPUT, Lambda, GRAD, HESSIAN] =  fmincon('swgarchlik', x0 ,A, b, [],[],lowerBounds, upperBounds, [], options, data, errortype, flag1, constr);
[LLF,likelihoods,e,p1,p1t,p1sm,h1,h2] = swgarchlik(parameters,data,errortype,flag1,constr);
..........

10# tulipsliu
恩，哥们儿。你说的这个人的似然函数我也看过。不过，我也仔细看过它的代码，大概写在5年前左右吧。它会打开一个文件句柄，fid=fopen('xxx.txt');我不知道会影响程序效率不，不过个人对他的程序在循环里多次这样调用，感觉会影响效率，所以非常反感。

DCC GACH 的估计比较复杂。用的是 2SQMLE 估计方法，也就是所谓的 "两阶段伪极大似然法"。
其中，你提到的这个部分，是它生成 e_t*η_t,这个我说起来也复杂。
而 Ht=Dt*Rt*Dt,
你在提到的这个用 fmincon 寻找的，是 Et 的模拟结果。
实际上，还有两个似然函数，一个simulation，所以DCC估计会异常的复杂。
一个是用 fmincon 寻找的dccparameters 的优化值（以logL最大化为最优），似然函数是 'dcc_mvgarch_likelihood'。
而且这一步可以有一个 海塞矩阵。可是依然不适合于DCC，参见后面参数估计时，还有一个专用的 dcc_hessian 函数。
另一个是第二阶段的极大化，用第一阶段估计得到的参数 dccparameters ，带入另一个似然函数‘dcc_mvgarch_full_likelihood’。其中，这个函数里还包含一个，'dcc_mvgarch_simulation'。这个用于生成没有dcc 效用的 Ht,而且返回这个 Ht，同时这个函数返回的logL,才是最后的。所以这个非常的复杂。

至于你提到的：
还是应该像marcucci,mrsgarchestfor_con_all.m
[parameters, LLF, EXITFLAG, OUTPUT, Lambda, GRAD, HESSIAN] =  fmincon('swgarchlik', x0 ,A, b, [],[],lowerBounds, upperBounds, [], options, data, errortype, flag1, constr);
[LLF,likelihoods,e,p1,p1t,p1sm,h1,h2] = swgarchlik(parameters,data,errortype,flag1,constr);

这个我也参考过，不过因为是否是 mvgarch，一个是 errortipe（T分布，还是GED）等。这个swgarchlik 对 mrs_dcc_mvgarch_likelihood的设计是没帮助的吧。毕竟方程组的设置就很不一样。

1. function [parameters, loglikelihood, Ht, Qt,  stdresid, likelihoods, stderrors, A,B, jointscores]=dcc_mvgarch(data,dccP,dccQ,archP,garchQ)
   
2. % PURPOSE:
   
3. %        Estimates a multivariate GARCH model using the DCC estimator of Engle and Sheppard
   
4. %
   
5. % USAGE:
   
6. %        [parameters, loglikelihood, Ht, Qt,  likelihoods, stdresid, stderrors, A,B, jointscores]...
   
7. %                =dcc_mvgarch(data,dccP,dccQ,archP,garchQ)
   
8. %
   
9. % INPUTS:
   
10. %      data          = A zero mean t by k vector of residuals from some filtration
    
11. %      dccP          = The lag length of the innovation term in the DCC estimator
    
12. %      dccQ          = The lag length of the lagged correlation matrices in the DCC estimator
    
13. %      archP         = One of three things:    Empty   in which case a 1 innovation model is estimated for each series
    
14. %                                      A scalar, p     in which case a p innovation model is estimated for each series
    
15. %                                      A k by 1 vector in which case the ith series has innovation terms p=archP(i)
    
16. %      garchQ        = One of three things:    Empty   in which case a 1 GARCH lag is used in estimation for each series
    
17. %                                      A scalar, q     in which case a q GARCH lags is used in estimation for each series
    
18. %                                      A k by 1 vector in which case the ith series has lagged variance terms q=archQ(i)
    
19. %
    
20. % OUTPUTS:
    
21. %      parameters    = A vector of parameters estimated form the model of the form
    
22. %                          [GarchParams(1) GarchParams(2) ... GarchParams(k) DCCParams]
    
23. %                          where the garch parameters from each estimation are of the form
    
24. %                          [omega(i) alpha(i1) alpha(i2) ... alpha(ip(i)) beta(i1) beta(i2) ... beta(iq(i))]
    
25. %      loglikelihood = The log likelihood evaluated at the optimum
    
26. %      Ht            = A k by k by t array of conditional variances
    
27. %      Qt            = A k by k by t array of Qt elements
    
28. %      likelihoods   = the estimated likelihoods t by 1
    
29. %      stderrors     = A length(parameters)^2 matrix of estimated correct standard errors
    
30. %      A             = The estimated A form the rebust standard errors
    
31. %      B             = The estimated B from the standard errors
    
32. %      scores        = The estimated scores of the likelihood t by length(parameters)
    
33. %
    
34. %
    
35. % COMMENTS:
    
36. %
    
37. %
    
38. % Author: Kevin Sheppard
    
39. % kevin.sheppard@economics.ox.ac.uk
    
40. % Revision: 2    Date: 12/31/2001
    
41. 
    
42. % Lets do some error checking and clean up
    
43. 
    
44. [t,k]=size(data);
    
45. 
    
46. if k<2
    
47.     error('Must have at least 2 data series')
    
48. end
    
49. 
    
50. if length(dccP)~=length(dccQ) || length(dccP)~=1
    
51.     error('dccP and dccQ must be scalars')
    
52. end
    
53. 
    
54. if ~(isempty(archP) || length(archP)==1 || length(archP)==k)
    
55.     error('Wrong size for archP')
    
56. end
    
57. 
    
58. if ~(isempty(garchQ) || length(garchQ)==1 || length(garchQ)==k)
    
59.     error('Wrong size for garchQ')
    
60. end
    
61. 
    
62. if isempty(archP)
    
63.     archP=ones(1,k);
    
64. elseif length(archP)==1
    
65.     archP=ones(1,k)*archP;
    
66. end
    
67. 
    
68. if isempty(garchQ)
    
69.     garchQ=ones(1,k);
    
70. elseif length(garchQ)==1
    
71.     garchQ=ones(1,k)*garchQ;
    
72. end
    
73. 
    
74. 
    
75. % Now lest do the univariate garching using fattailed_garch as it's faster then garchpq
    
76. stdresid=data;
    
77. options=optimset('fmincon');
    
78. options=optimset(options,'Display','off','Diagnostics','off','MaxFunEvals',1000*max(archP+garchQ+1),'MaxIter',1000*max(archP+garchQ+1),'LargeScale','off','MaxSQPIter',1000);
    
79. options  =  optimset(options , 'MaxSQPIter' , 1000);
    
80. for i=1:k
    
81.     fprintf(1,'Estimating GARCH model for Series %d\n',i)
    
82.     [univariate{i}.parameters, univariate{i}.likelihood, univariate{i}.stderrors, univariate{i}.robustSE, univariate{i}.ht, univariate{i}.scores] ...
    
83.         = fattailed_garch(data(:,i) , archP(i) , garchQ(i) , 'NORMAL',[], options);
    
84.     stdresid(:,i)=data(:,i)./sqrt(univariate{i}.ht);
    
85. end
    
86. 
    
87. 
    
88. options=optimset('fmincon');
    
89. options  =  optimset(options , 'Display'     , 'iter');
    
90. options  =  optimset(options , 'Diagnostics' , 'on');
    
91. options  =  optimset(options , 'LevenbergMarquardt' , 'on');
    
92. options  =  optimset(options , 'LargeScale'  , 'off');
    
93. 
    
94. 
    
95. dccstarting=[ones(1,dccP)*.01/dccP ones(1,dccQ)*.97/dccQ];
    
96. fprintf(1,'\n\nEstimating the DCC model\n')
    
97. 
    
98. [dccparameters,~,~,~,~,GRAD]=fmincon('dcc_mvgarch_likelihood',dccstarting,ones(size(dccstarting)),1-2*options.TolCon,[],[],zeros(size(dccstarting))+2*options.TolCon,[],[],options,stdresid,dccP,dccQ);
    
99. 
    
100. 
     
101. 
     
102. % We now have all of the estimated parameters
     
103. parameters=[];
     
104. H=zeros(t,k);
     
105. for i=1:k
     
106.     parameters=[parameters;univariate{i}.parameters];
     
107.     H(:,i)=univariate{i}.ht;
     
108. end
     
109. parameters=[parameters;dccparameters'];
     
110. 
     
111. 
     
112. 
     
113. %We now have Ht and the likelihood
     
114. [loglikelihood, Rt, likelihoods, Qt]=dcc_mvgarch_full_likelihood(parameters, data, archP,garchQ,dccP,dccQ);
     
115. likelihoods=-likelihoods;
     
116. loglikelihood=-loglikelihood;
     
117. Ht=zeros(k,k,t);
     
118. stdresid=zeros(t,k);
     
119. Hstd=H.^(0.5);
     
120. for i=1:t
     
121.     Ht(:,:,i)=diag(Hstd(i,:))*Rt(:,:,i)*diag(Hstd(i,:));
     
122.     stdresid(i,:)=data(i,:)*Ht(:,:,i)^(-0.5);
     
123. end
     
124. save tempHt Ht
     
125. clear Ht
     
126. 
     
127. if nargout >=7
     
128.     %How ar we going to get STD errors?  Partitioned invers probably.  Well, we need to get the scores form the dcc model, the joint likelihood.
     
129.     %We then need to get A12 and A22 so we can have it all.  We also need to get A11 in the correct form.
     
130.     A=zeros(length(parameters),length(parameters));
     
131.     index=1;
     
132.     for i=1:k
     
133.         workingsize=size(univariate{i}.stderrors);
     
134.         A(index:index+workingsize-1,index:index+workingsize-1)=univariate{i}.stderrors^(-1);
     
135.         index=index+workingsize;
     
136.     end
     
137.    
     
138.     % Ok so much for a All and A12 and A22, as we have them all between whats above
     
139.     fprintf(1,'\n\nCalculating Standard Errors, this can take a while\n');
     
140.     otherA=dcc_hessian('dcc_mvgarch_full_likelihood',parameters, dccP+dccQ, data, archP,garchQ,dccP,dccQ);
     
141.     A(length(parameters)-dccP-dccQ+1:length(parameters),:)=otherA;
     
142.     % tempA=hessian('dcc_garch_full_likelihood',parameters, data, archP,garchQ,dccP,dccQ);
     
143.     % A(length(parameters)-1:length(parameters),:)=tempA(length(parameters)-1:length(parameters),:);
     
144.     %That finishes A
     
145.    
     
146.     % We now need to get the scores for the DCC estimator so we can finish B
     
147.     jointscores=zeros(t,length(parameters));
     
148.     index=1;
     
149.     for i=1:k
     
150.         workingsize=size(univariate{i}.scores,2);
     
151.         jointscores(:,index:index+workingsize-1)=univariate{i}.scores;
     
152.         index=index+workingsize;
     
153.     end
     
154.    
     
155.     %Now all we need to do is calculate the scores form teh dcc estimator and we have everything
     
156.     h=max(abs(parameters/2),1e-2)*eps^(1/3);
     
157.     hplus=parameters+h;
     
158.     hminus=parameters-h;
     
159.     likelihoodsplus=zeros(t,length(parameters));
     
160.     likelihoodsminus=zeros(t,length(parameters));
     
161.     for i=length(parameters)-dccP-dccQ+1:length(parameters)
     
162.         hparameters=parameters;
     
163.         hparameters(i)=hplus(i);
     
164.         [~, ~, indivlike] = dcc_mvgarch_full_likelihood(hparameters, data, archP,garchQ,dccP,dccQ);
     
165.         likelihoodsplus(:,i)=indivlike;
     
166.     end
     
167.     for i=length(parameters)-dccP-dccQ+1:length(parameters)
     
168.         hparameters=parameters;
     
169.         hparameters(i)=hminus(i);
     
170.         [~, ~, indivlike] = dcc_mvgarch_full_likelihood(hparameters, data, archP,garchQ,dccP,dccQ);
     
171.         likelihoodsminus(:,i)=indivlike;
     
172.     end
     
173.     DCCscores=(likelihoodsplus(:,length(parameters)-dccP-dccQ+1:length(parameters))-likelihoodsminus(:,length(parameters)-dccP-dccQ+1:length(parameters)))...
     
174.         ./(2*repmat(h(length(parameters)-dccP-dccQ+1:length(parameters))',t,1));
     
175.     jointscores(:,length(parameters)-dccP-dccQ+1:length(parameters))=DCCscores;
     
176.     B=cov(jointscores);
     
177.     A=A/t;
     
178.     stderrors=A^(-1)*B*A'^(-1)*t^(-1);
     
179. end
     
180. %Done!
     
181. load tempHt

复制代码

那几个函数我用不同颜色标注的，怎么发出来就没了。郁闷。
这个就是 dcc_mvgarch 的主程序。

??? Error: File: C:\MarkovSwitchingDCC\ms_dcc_mvgarch.m Line: 100 Column: 17
Incomplete or misformed expression or statement.

Error in ==> Demo_MS_DCC at 26
[parameters, loglikelihood, Ht, Qt,  stdresid, likelihoods, smooth,Pr]=...


??? Error: File: C:\MarkovSwitchingDCC\ms_dcc_mvgarch.m Line: 100 Column: 17
Incomplete or misformed expression or statement.

??? Error: File: C:\MarkovSwitchingDCC\ms_dcc_mvgarch_full_likelihood.m Line: 21 Column: 7
Incomplete or misformed expression or statement.

??? Error: File: C:\MarkovSwitchingDCC\ms_dcc_mvgarch_full_likelihood.m Line: 21 Column: 7
Incomplete or misformed expression or statement.

??? Input argument "stdresid" is undefined.

Error in ==> ms_dcc_mvgarch_likelihood at 10
[t,k]=size(stdresid);

??? Input argument "x" is undefined.

Error in ==> tvp_markov_filter at 53
[n k] = size(x);
tulipsliu 朋友，你的程序在matlab7上运行时经常提示错误，修改完一个后，又是一个，
比Gauss程序麻烦多了，不知你在运行时有以上错误吗？如何修改？
我因为给学生上《数学建模》，所以懂一些matlab，但你的程序我打印出阅读，非常难懂，
不知你是否有什么建议如何能容易读懂？
我认为你是论坛中matlab方面少有的高手，希望能多帮助和指教，谢谢！
如果涉及到知识产权，希望你发到邮箱：zhangws027@163.com,
我将严格遵守知识产权法，谢谢！

15# zhangtao
朋友，我运行的时候肯定很多错误，不只时候语法错误，还包括逻辑错误（状态空间转移概率竟然为1.4，我想，概率论说的是 0~1 ，我的程序算出来有到1.4的，你应该知道是错的）。
当然，上传论坛的程序，有逻辑错误，应该没语法错误。我再检验试试，可以运行的。在我的电脑里，算是最近修改的程序里最恰当的一个。至于你说的 stdresid 没有定义，x 没有定义。我可以这样猜测，因为你的版本是 7.0 版本的。里面没有计量经济学工具箱。我调用的数据不是我整理的，你应该看到我的程序包里没有数据。而是 MATLAB 系统自己带的demo 里的数据。是西方主要工业国的股票证券指数，包括 FTS100(英国）、CAC(法国），S&P500（标普500）。Econometrics toolbox ，我没记错的话，是在 2008A 版本后才有，7.0 应该比较早，是2005还是哪一年的版本呢。所以你的程序里没有这个数据的。这个是数据问题，怪你的版本太低。我现在用的是 MATLAB R2011a 版本。呵呵，计量经济学工具箱里还扩增了 “协整检验”。

要解决这个问题，你自己重新调入数据。这个程序应该可以运行的，当时结果不对。好几个地方不对。fattail_garch()这里，依然用 KEVEN 的函数，这个没有问题。只是系统会提示 fmincon 里的算法可能不恰当。
后面的 mrs_dcc_mvgarch_likelihood,这个会被 fmincon 调用，这个是第一步的对 dccparameters 的优化（因为DCC MVGARCH 用的是2SQMLE）。fmincon 只对这个优化。它的问题是，优化的数值是复数的。

另一个似然函数是第二阶段的，叫： mrs_dcc_mvgarhc_full_likelihood。这个是另一个函数，这个函数基本就可以估计出最重要的参数了。然后才是主估计程序：mrs_dcc_mvgarch，这个部分，因为Ht无法得到正确的估计值。所以还没想KEVEN那样，把后面其他参数给估计出来。比如： A、B、 jiontscores 等参数。这些还没写。因为第一步的数值不对，或者叫，有bug。所以在论坛发帖，求大家帮忙。把程序修正。

MATLAB 算是比较友好的语言，如果是用C语言来做算法。我看过用C做的金融建模程序的，那个量很大。而且，如果还是可以用类库的“面向对象的程序设计”，C++ 这一类，就更加难以理解。
至于这个程序难以理解。不只是这个。UCSD_GARHC 里，还有KEVEN 自己的 dcc_mvgach。你可以参阅KEVEN写的。因为DCC GARCH 算是比较难的模型了。所以写的程序代码比较难懂。调用的函数太多。

KEVEN 是engel的学生，好像现在在牛津大学上课，他还写了另一个工具箱 MFE toolbox。这个工具箱2011年有更新，我下载了这个工具箱了。。里面还包括前沿的“高频数据建模”。金融计量经济学里的以每一笔证券成交的高频数据，来建模的模型，有 ACD,RV（已实现波动模型）。国内的教材，我看的是张世英的。其他人的我还没看。
其中，ACD模型的估计，可以看 Perlin的程序，他也做Markov regime switchng  的模型，不过他的程序不包括GARCH 类型的，所以无法直接用。他的程序有问题，就是有些地方写得不好。也不好阅读。这个是一个诟病。

KEVEN 的程序写得非常漂亮，其实算是最好阅读，也是最好理解和修改的。KEVEN 的最新MFE 里，添加了对 RV 模型的估计。这个你可以看看，非常的精彩。其实不只KEVEN 。UCSD，应该是 芝加哥大学，圣迭戈。这里还有一个 Lesage，我和朋友讨论空间计量经济学时用到他的程序，写得非常好，格林的高级教程《计量经济学》里的很多，他都写出来了，包括那些特殊的分布，比如 维希特分布，当然，现在的matlab高版本里，系统也提供这个分布类型的函数。Lesage 的一个很重要的程序，我看人的名字，好像2005年，也有一个中国人帮他写了程序。非常不错。

KEVEN  的2011版本的工具箱已经包含这个，下载这个，开个玩笑，已经和国际接轨，和牛津大学的同学同步了。呵呵。

原来您是老师啊，您好。谢谢您的参与。我最近加了几个朋友，他们还是在读博士和硕士。其中有4个和我一起也在看这个 MRS DCC GARCH 。不过他们也没找到解决的办法。

容我说白一点
MRS-DCC-GARCH model
第一步是估计一般garch model 的系数 ?
还是应该估计mrs-garch model 的系数 ?

请看page 10/46
The estimation are done in a two-step procedure.
    First, for each of the two time series, a univariate
           two-regime MS-GARCH model (2) is estimated.
    Second,Then with the parameters fixed at the estimates from the previous stage,
           the rest of the parameters of model(3) are estimated.

Regime Switching in Volatilities and Correlation between Stock and Bond markets.pdf
Regime_Switching_in_Volatilities.pdf (1.46 MB)

2011-6-27 09:39:01 上传

至于Marcucci MRS-GARCH
大家都说跑得慢,
最主要是它采rolling window
只要把N设为T,N=T
他的文献 page 31/42,MRS-GARCH-N
约十来分钟就可跑完.
跑出的结果,跟他的文献差距不大
0.0534 -1.7698 0.0071 0.5304 0.0283 0.0265 0.9081 0.9587 0.9815 0.1658

11# epoh
老兄,刚看了marcucci和你的 source code
发现有一个问题,值得思考与讨论
你的ms_dcc_mvgarch.m是
for i=1:k
     fprintf(1,'Estimating GARCH model for Series %d\n',i)
     [univariate{i}.parameters, univariate{i}.likelihood, univariate{i}.stderrors, univariate{i}.robustSE, univariate{i}.ht, univariate{i}.scores] ...
         = fattailed_garch(data(:,i) , archP(i) , garchQ(i) , 'NORMAL',[], options);
     stdresid(:,i)=data(:,i)./sqrt(univariate{i}.ht);
end

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

还是应该像marcucci,mrsgarchestfor_con_all.m
[parameters, LLF, EXITFLAG, OUTPUT, Lambda, GRAD, HESSIAN] =  fmincon('swgarchlik', x0 ,A, b, [],[],lowerBounds, upperBounds, [], options, data, errortype, flag1, constr);
[LLF,likelihoods,e,p1,p1t,p1sm,h1,h2] = swgarchlik(parameters,data,errortype,flag1,constr);

是的，你看到这里的问题了，不过这里也要修改的话，程序会更复杂了。
其中：stdresid=univariate{i}.stderrors;这个会被后面的参数优化带入，  [dccpara,~,~]=fmincon(@mrs_dcc_mvgarch_liklihoods(para,stdresid,……）……）；
也就是这个也得考虑用 swgarch，替代这里的 fattailed_garch（）函数的问题。不知道这样理解得对不。