Markov regime switching Model

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

最近朋友找我帮忙做一个计量模型的 Markov regime swithing 的模型，我想了一周，也没改造成功，看到论坛那么多有才的人，比如epoh，我这里发一个贴，期待各位的回应。看看能解决这个问题不。这里发一个arch_3_state的状态空间模型，期待各位的回复。

1. %% coed_MS_ARCH_3S.txt
   
2. %-------------------------------------------
   
3. % Written by Zhu, Junjun,
   
4. % Ph.D candidate in School of Economics
   
5. % Fudan University, Shanghai 200433
   
6. % Nov. 2009
   
7. 
   
8. % please notice, we can not guarantee that these codes are without mistakes.
   
9. % You may use these at your own risks.
   
10. %-------------------------------------------
    
11. 
    
12. clc;
    
13. clear;
    
14. load data_200910weekex;
    
15. y = data_200910weekex;
    
16. n = size(y,1);
    
17. x = ones(n,1);
    
18. sn = 3; % number of state
    
19. nk = 1; % regressor number in mean equation
    
20. np = 3; % p-term in GARCH
    
21. y_p = y(np+1:n);
    
22. x_p = x(np+1:n);
    
23. pi_pri = [18 2 2; 2 28 2; 2 2 16];  % prior for pi
    
24. bols = inv(x'*x)*x'*y;
    
25. s2ols = (y-x*bols)'*(y-x*bols)/(n-1);
    
26. xsquare=x'*x;
    
27. indexSwitch =0;
    
28. switch_1 = 0;
    
29. switch_2 = 0;
    
30. switch_3 = 0;
    
31. P=[0.8 0.1 0.1;0.1 0.8 0.1; 0.1 0.1 0.8];
    
32. %spec = garchset('distribution', 'T','p',0,'q',3);
    
33. %[coeff, errors,LLF,innovations,sigmas,summary] = garchfit(spec,y);
    
34. %garchdisp(coeff,errors)
    
35. % coeff.C = [0.0005;]; coeff.K = 0.00015; coeff.ARCH = [0.197,0.26,0.292;]
    
36. % default model: y(t) = C +e(t); h(t) = K + GARCH*h(t-1) + ARCH*e(t-1)^2
    
37. %   beta = [0.0016; 0.0016; 0.0016];  
    
38. %   a_1 = [0.00074; 0.2; 0.26; 0.29]; % ARCH coeff. a0, a1, a2, a3 for state 1
    
39. %   a_2 = [0.00074; 0.2; 0.26; 0.29];
    
40. %   a_3 = [0.00074; 0.2; 0.26; 0.29];
    
41. beta = [-0.0011; 0.02; -0.0061];  
    
42. a_1 = [0.0005; 0.11; 0.15; 0.08];
    
43. a_2 = [0.0014; 0.22; 0.17; 0.13];
    
44. a_3 = [0.0022; 0.14; 0.10; 0.15];
    
45. var_err = (y-x*beta(1)).^2;
    
46. 
    
47. beta_ = [];
    
48. a1_ = [];
    
49. a2_ = [];
    
50. a3_ = [];
    
51. P_ = [];
    
52. nburn=0; %number of burnin replications
    
53. nrep=30000; %number of retained replications
    
54. ntotal=nburn+nrep;
    
55. 
    
56. 
    
57. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    
58. %-----------------------------------------------------------------------
    
59. %starting value for S, P
    
60. %-----------------------------------------------------------------------
    
61.     S=zeros(n,1);
    
62.    
    
63.     var_err_mean = mean(var_err);
    
64.     for t=1:n
    
65.         if var_err(t)<=var_err_mean-0.00005              % state 1 as low volatility state
    
66.             S(t)=1;
    
67.         else
    
68.             if var_err(t)<=var_err_mean+0.00015
    
69.             S(t)=2;
    
70.             else
    
71.             S(t) = 3;
    
72.             end
    
73.         end
    
74.     end
    
75.    
    
76. %----------------------------------------------------------------------
    
77. % starting value for P
    
78. %----------------------------------------------------------------------
    
79.     T = n;
    
80.     sm=[S(1:T-1) S(2:T)];  % switching number, sm: (T-1)*2
    
81.     tm=zeros(sn,sn);   % total number
    
82.     for im=1:sn
    
83.        for jm=1:sn
    
84.         tm(im,jm)=size(sm(sm(:,1)==im & sm(:,2)==jm,:),1); % total number of each transition from i state to j state              
    
85.        end
    
86.     rr = pi_pri(im,:)+tm(im,:);
    
87.     P_draw = chi2rnd(2*rr);
    
88.     P(im,:)=P_draw./repmat(sum(P_draw,2),1,sn);
    
89.     % P(im,:)=draw_dirichlet(1,sn,pi_pri(im,:)+tm(im,:),0); % rm: prior for number of state transition (ij)  
    
90.     end
    
91. % 第一部分结束，

复制代码

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏6 回帖

关键词：switching Markov regime switch Ching model Markov regime switching

1. [code]%紧接第一部分代码
   
2. 
   
3. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Now start Gibbs loop%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
4. for i = 1:ntotal
   
5. completion_process=i/ntotal
   
6. %-------------------------------------------------------
   
7. % caculate draw probability
   
8. %---------------S1---------------------------------------
   
9. %log likelihood for the regression model with ARCH errors -- old draw
   
10. y_1 = y(S==1); x_1 = x(S==1);
    
11. err_1 = y_1 - x_1*beta(1);
    
12. n1 = length(err_1);
    
13. h1 =a_1(1)*ones(n1-np,1); % vols: 524*1
    
14. for j=1:np
    
15. h1 = h1 + a_1(j+1,1)*err_1(np+1-j:n1-j,1).^2; % conditional variance h
    
16. end
    
17. % alternative: vols = vols + adraw(2,1)*errs(p:t-1,1) + adraw(3,1)*errs(p-1:t-2) + ...
    
18. % adraw(4,1)*errs(p-2:t-3,1) + adraw(5,1)*errs(p-3:t-4,1)
    
19. lterm=(err_1(np+1:n1,1).^2)./h1;
    
20. llike1 = -.5*sum(log(h1)) -.5*sum(lterm);
    
21. 
    
22. %---------------S2---------------------------------------
    
23. y_2 = y(S==2); x_2 = x(S==2);
    
24. err_2 = y_2 - x_2*beta(2);
    
25. n2 = length(err_2);
    
26. h2 =a_2(1)*ones(n2-np,1);% vols: 524*1
    
27. for j=1:np
    
28. h2 = h2 + a_2(j+1,1)*err_2(np+1-j:n2-j,1).^2; % conditional variance h
    
29. end
    
30. lterm=(err_2(np+1:n2,1).^2)./h2;
    
31. llike2 = -.5*sum(log(h2)) -.5*sum(lterm);
    
32. 
    
33. %---------------S3---------------------------------------
    
34. y_3 = y(S==3); x_3 = x(S==3);
    
35. err_3 = y_3 - x_3*beta(3);
    
36. n3 = length(err_3);
    
37. h3 =a_3(1)*ones(n3-np,1);% vols: 524*1
    
38. for j=1:np
    
39. h3 = h3 + a_3(j+1,1)*err_3(np+1-j:n3-j,1).^2; % conditional variance h
    
40. end
    
41. lterm=(err_3(np+1:n3,1).^2)./h3;
    
42. llike3 = -.5*sum(log(h3)) -.5*sum(lterm);
    
43. 
    
44. %---------------------------------------------------------
    
45. % take candidate draw, caculate candidate probability
    
46. %---------------------------------------------------------
    
47. %take candidate draw
    
48. beta_sig = [(0.003/4)^2 0 0; 0 (0.013/4)^2 0; 0 0 (0.018/4)^2 ];
    
49. beta_can = beta + norm_rnd(beta_sig);
    
50. % a_1_sig = [(0.00074/6)^2 0 0 0;0 (0.1/6)^2 0 0; 0 0 (0.1/6)^2 0; 0 0 0 (0.2/6)^2];
    
51. a_1_sig = [(0.0001/4)^2 0 0 0;0 (0.075/4)^2 0 0; 0 0 (0.07/4)^2 0; 0 0 0 (0.647/4)^2];
    
52. a_1_can = a_1 + norm_rnd(a_1_sig);
    
53. % a_2_sig = [(0.00074/4)^2 0 0 0;0 (0.1/4)^2 0 0; 0 0 (0.1/4)^2 0; 0 0 0 (0.2/4)^2];
    
54. a_2_sig = [(0.0007/4)^2 0 0 0;0 (0.1406/4)^2 0 0; 0 0 (0.13/4)^2 0; 0 0 0 (0.12/4)^2];
    
55. a_2_can = a_2 + norm_rnd(a_2_sig);
    
56. % a_3_sig = [(0.00074/4)^2 0 0 0;0 (0.1/4)^2 0 0; 0 0 (0.1/4)^2 0; 0 0 0 (0.2/4)^2];
    
57. a_3_sig = [(0.0011/4)^2 0 0 0;0 (0.1473/4)^2 0 0; 0 0 (0.0855/4)^2 0; 0 0 0 (0.1824/4)^2];
    
58. a_3_can = a_3 + norm_rnd(a_3_sig);
    
59. %if this draw violates positivity or stationarity restrictions do not accept it
    
60. violat_1 =0; violat_2 =0; violat_3 =0;
    
61. %first check positivity of arch coefficients
    
62. for j=1:np+1
    
63. if a_1_can(j,1)<=0
    
64. violat_1=1;
    
65. end
    
66. if a_2_can(j,1)<=0
    
67. violat_2=1;
    
68. end
    
69. if a_3_can(j,1)<=0
    
70. violat_3=1;
    
71. end
    
72. end
    
73. if a_1_can(2) + a_1_can(3) +a_1_can(4) >1; violat_1=1; end
    
74. if a_2_can(2) + a_2_can(3) +a_2_can(4) >1; violat_2=1; end
    
75. if a_3_can(2) + a_3_can(3) +a_3_can(4) >1; violat_3=1; end
    
76. 
    
77. %-------------------------------------------------------
    
78. % Accept the candidate with probability
    
79. %-------------------------------------------------------
    
80. if violat_1 ==0
    
81. %log likelihood for the regression model with ARCH errors -- candidate draw
    
82. err_1_can = y_1 - x_1*beta_can(1);
    
83. h1_can=a_1_can(1,1)*ones(n1-np,1);
    
84. for j=1:np
    
85. h1_can = h1_can + a_1_can(j+1,1)*err_1_can(np+1-j:n1-j,1).^2;
    
86. end
    
87. lterm=(err_1_can(np+1:n1,1).^2)./h1_can;
    
88. llike1_p = -.5*sum(log(h1_can)) -.5*sum(lterm);
    
89. %log of acceptance probability
    
90. laccprob_1 = llike1_p - llike1;
    
91. %now accept candidate with appropriate probability
    
92. if laccprob_1 > log(rand)
    
93. beta(1)=beta_can(1);
    
94. a_1 = a_1_can;
    
95. switch_1 = switch_1+1;
    
96. end
    
97. end
    
98. %-------------------------- S2--------------------------------------------
    
99. if violat_2 ==0
    
100. %log likelihood for the regression model with ARCH errors -- candidate draw
     
101. err_2 = y_2 - x_2*beta_can(2);
     
102. h2=a_2_can(1,1)*ones(n2-np,1);
     
103. for j=1:np
     
104. h2 = h2 + a_2_can(j+1,1)*err_2(np+1-j:n2-j,1).^2;
     
105. end
     
106. lterm=(err_2(np+1:n2,1).^2)./h2;
     
107. llike2_p = -.5*sum(log(h2)) -.5*sum(lterm);
     
108. %log of acceptance probability
     
109. laccprob_2 = llike2_p - llike2;
     
110. %now accept candidate with appropriate probability
     
111. if laccprob_2 > log(rand)
     
112. beta(2)=beta_can(2);
     
113. a_2 = a_2_can;
     
114. switch_2 = switch_2+1;
     
115. end
     
116. end
     
117. 
     
118. %-------------------------- S2--------------------------------------------
     
119. if violat_3 ==0
     
120. %log likelihood for the regression model with ARCH errors -- candidate draw
     
121. err_3 = y_3 - x_3*beta_can(3);
     
122. h3=a_3_can(1,1)*ones(n3-np,1);
     
123. for j=1:np
     
124. h3 = h3 + a_3_can(j+1,1)*err_3(np+1-j:n3-j,1).^2;
     
125. end
     
126. lterm=(err_3(np+1:n3,1).^2)./h3;
     
127. llike3_p = -.5*sum(log(h3)) -.5*sum(lterm);
     
128. %log of acceptance probability
     
129. laccprob_3 = llike3_p - llike3;
     
130. %now accept candidate with appropriate probability
     
131. if laccprob_3 > log(rand)
     
132. beta(3)=beta_can(3);
     
133. a_3 = a_3_can;
     
134. switch_3 = switch_3+1;
     
135. end
     
136. end
     
137. 
     
138. %--------------------------------------------------------------------
     
139. %draw S conditional on [lnl,filprob,simstate] = MSLinRegDrawState(T,k1,k2,nm,indSigma,p,P,yv,X1,X2,Beta,hv);
     
140. %--------------------------------------------------------------------
     
141. err_S = zeros(n,sn);
     
142. h_S= zeros(n,sn);
     
143. for j=1:sn
     
144. for i = 1:n
     
145. err_S(i,j) = y(i) - x(i)*beta(j);
     
146. end
     
147. end
     
148. h_S(1,1) = a_1(1)/(1-a_1(2)-a_1(3)-a_1(4));
     
149. h_S(2,1) = a_1(1) + a_1(2)*err_S(1,1)^2;
     
150. h_S(3,1) = a_1(1) + a_1(2)*err_S(2,1)^2+a_1(3)*err_S(1,1)^2;
     
151. h_S(1,2) = a_2(1)/(1-a_2(2)-a_2(3)-a_2(4));
     
152. h_S(2,2) = a_2(1) + a_2(2)*err_S(1,2)^2;
     
153. h_S(3,2) = a_2(1) + a_2(2)*err_S(2,2)^2+a_2(3)*err_S(1,2)^2;
     
154. h_S(1,3) = a_3(1)/(1-a_3(2)-a_3(3)-a_3(4));
     
155. h_S(2,3) = a_3(1) + a_3(2)*err_S(1,3)^2;
     
156. h_S(3,3) = a_3(1) + a_3(2)*err_S(2,3)^2+a_3(3)*err_S(1,3)^2;
     
157. for j = 4:n
     
158. h_S(j,1) = a_1(1) + a_1(2)*err_S(j-1,1)^2 + a_1(3)*err_S(j-2,1)^2 + a_1(4)*err_S(j-3,1)^2;
     
159. h_S(j,2) = a_2(1) + a_2(2)*err_S(j-1,2)^2 + a_2(3)*err_S(j-2,2)^2 + a_2(4)*err_S(j-3,2)^2;
     
160. h_S(j,3) = a_3(1) + a_3(2)*err_S(j-1,3)^2 + a_3(3)*err_S(j-2,3)^2 + a_3(4)*err_S(j-3,3)^2;
     
161. end
     
162. % (normpdf) lnlike = exp(-0.5 * ((x - mu)./sigma).^2) ./ (sqrt(2*pi) .* sigma);
     
163. lnpdat(:,1) = -0.5*log(2*pi*h_S(:,1))-0.5*(err_S(:,1).^2)./h_S(:,1);
     
164. lnpdat(:,2) = -0.5*log(2*pi*h_S(:,2))-0.5*(err_S(:,2).^2)./h_S(:,2);
     
165. lnpdat(:,3) = -0.5*log(2*pi*h_S(:,3))-0.5*(err_S(:,3).^2)./h_S(:,3);
     
166. lnpdat(:,1) = exp(lnpdat(:,1));
     
167. lnpdat(:,2) = exp(lnpdat(:,2));
     
168. lnpdat(:,3) = exp(lnpdat(:,3));
     
169. filprob=zeros(n,3);
     
170. k=S(2);
     
171. piSti1(1) = P(1,k)*lnpdat(1,1)/(P(1,k)*lnpdat(1,1)+P(2,k)*lnpdat(1,2)+P(3,k)*lnpdat(1,3));
     
172. piSti2(1) = P(2,k)*lnpdat(1,2)/(P(1,k)*lnpdat(1,1)+P(2,k)*lnpdat(1,2)+P(3,k)*lnpdat(1,3));
     
173. piSti3(1) = P(3,k)*lnpdat(1,3)/(P(1,k)*lnpdat(1,1)+P(2,k)*lnpdat(1,2)+P(3,k)*lnpdat(1,3));
     
174. random_S = rand;
     
175. if random_S<=piSti1(1)
     
176. S(1) = 1;
     
177. else
     
178. if random_S <=piSti1(1)+piSti2(1)
     
179. S(1)=2;
     
180. else
     
181. S(1) = 3;
     
182. end
     
183. end
     
184. 
     
185. filprob(1,1)=piSti1(1);
     
186. filprob(1,2)=piSti2(1);
     
187. filprob(1,3)=piSti3(1);
     
188. 
     
189. for t =2:n-1
     
190. l=S(t-1);
     
191. k=S(t+1);
     
192. piSti1(t) = P(l,1)*P(1,k)*lnpdat(t,1)/(P(l,1)*P(1,k)*lnpdat(t,1)+P(l,2)*P(2,k)*lnpdat(t,2)+P(l,3)*P(3,k)*lnpdat(t,3));
     
193. piSti2(t) = P(l,2)*P(2,k)*lnpdat(t,2)/(P(l,1)*P(1,k)*lnpdat(t,1)+P(l,2)*P(2,k)*lnpdat(t,2)+P(l,3)*P(3,k)*lnpdat(t,3));
     
194. piSti3(t) = P(l,3)*P(3,k)*lnpdat(t,3)/(P(l,1)*P(1,k)*lnpdat(t,1)+P(l,2)*P(2,k)*lnpdat(t,2)+P(l,3)*P(3,k)*lnpdat(t,3));
     
195. random_St = rand;
     
196. if random_St<=piSti1(t)
     
197. S(t) = 1;
     
198. else
     
199. if random_St<=piSti1(t)+piSti2(t)
     
200. S(t)=2;
     
201. else
     
202. S(t)=3;
     
203. end
     
204. end
     
205. filprob(t,1)=piSti1(t);
     
206. filprob(t,2)=piSti2(t);
     
207. filprob(t,3)=piSti3(t);
     
208. end
     
209. 
     
210. l=S(n-1);
     
211. piSti1(n) = P(l,1)*lnpdat(n,1)/(P(l,1)*lnpdat(n,1)+P(l,2)*lnpdat(n,2)+P(l,3)*lnpdat(n,3));
     
212. piSti2(n) = P(l,2)*lnpdat(n,2)/(P(l,1)*lnpdat(n,1)+P(l,2)*lnpdat(n,2)+P(l,3)*lnpdat(n,3));
     
213. piSti3(n) = P(l,3)*lnpdat(n,3)/(P(l,1)*lnpdat(n,1)+P(l,2)*lnpdat(n,2)+P(l,3)*lnpdat(n,3));
     
214. random_Sn = rand;
     
215. if random_Sn<=piSti1(n)
     
216. S(n) = 1;
     
217. else
     
218. if random_Sn<=piSti1(n)+piSti2(n)
     
219. S(n)=2;
     
220. else
     
221. S(n) =3;
     
222. end
     
223. end
     
224. 
     
225. filprob(n,1)=piSti1(n);
     
226. filprob(n,2)=piSti2(n);
     
227. filprob(n,3)=piSti3(n);
     
228. %--------------------------------------------------------------------
     
229. %draw P conditional on S
     
230. %--------------------------------------------------------------------
     
231. sm=[S(1:n-1) S(2:n)]; % switching number, sm: (T-1)*2
     
232. tm=zeros(sn,sn); % total number
     
233. for im=1:sn
     
234. for jm=1:sn
     
235. tm(im,jm)=size(sm(sm(:,1)==im & sm(:,2)==jm,:),1); % total number of each transition from i state to j state
     
236. end
     
237. rr = pi_pri(im,:)+tm(im,:);
     
238. P_draw = chi2rnd(2*rr);
     
239. P(im,:)=P_draw./repmat(sum(P_draw,2),1,sn);
     
240. % P(im,:)=draw_dirichlet(1,sn,pi_pri(im,:)+tm(im,:),0); % rm: prior for number of state transition (ij)
     
241. end
     
242. 
     
243. 
     
244. %-------------------------------------------------------
     
245. %identification constraint
     
246. %-------------------------------------------------------
     
247. indexconstraint1=a_1(1)>a_2(1);
     
248. %if indexconstraint1 & abs(beta(1)-beta(2))<0.002
     
249. if a_1(1)-a_2(1)>0.0003 && beta(1)-beta(2)>0.005
     
250. indexSwitch=indexSwitch+1;
     
251. ndxv = [2 1 3];
     
252. ABSETZ = a_2;
     
253. a_2 = a_1;
     
254. a_1 = ABSETZ;
     
255. beta=beta(ndxv);
     
256. P=P(:,ndxv);
     
257. P=P(ndxv,:);
     
258. filprob=filprob(:,ndxv);
     
259. end
     
260. 
     
261. % indexconstraint2=beta(1)>beta(3)>beta(2);
     
262. indexconstraint2=a_1(1)>a_3(1);
     
263. %if indexconstraint2 & abs(beta(1)-beta(3))<0.002
     
264. if a_1(1)-a_3(1)>0.0003 && abs(beta(1)-beta(3))>0.002
     
265. indexSwitch=indexSwitch+1;
     
266. ndxv = [3 2 1];
     
267. ABSETZ = a_3;
     
268. a_3 = a_1;
     
269. 
     
270. a_1 = ABSETZ;
     
271. beta=beta(ndxv);
     
272. P=P(:,ndxv);
     
273. P=P(ndxv,:);
     
274. filprob=filprob(:,ndxv);
     
275. end
     
276. 
     
277. indexconstraint3=a_2(1)>a_3(1);
     
278. if indexconstraint3 && beta(3)-beta(2)>0.003
     
279. indexSwitch=indexSwitch+1;
     
280. ndxv = [1 3 2];
     
281. ABSETZ = a_3;
     
282. a_3 = a_2;
     
283. a_2 = ABSETZ;
     
284. beta=beta(ndxv);
     
285. P=P(:,ndxv);
     
286. P=P(ndxv,:);
     
287. filprob=filprob(:,ndxv);
     
288. end
     
289. 
     
290. %-----------------------------------------------------
     
291. % save accepted draw after burn-in replications
     
292. %-----------------------------------------------------
     
293. 
     
294. if i>nburn
     
295. beta_ = [beta_ beta];
     
296. a1_ = [a1_ a_1];
     
297. a2_ = [a2_ a_2];
     
298. a3_ = [a3_ a_3];
     
299. Pvec = [P(1,:)' ; P(2,:)'; P(3,:)'];
     
300. P_ = [P_ Pvec];
     
301. end
     
302. end

复制代码

1. MS3.P=P;
   
2. MS3.P_=P_;
   
3. MS3.S=S;
   
4. MS3.beta_=beta_;
   
5. MS3.a1_=a1_;
   
6. MS3.a2_=a2_;
   
7. MS3.a3_=a3_;
   
8. MS3.err_S=err_S;
   
9. MS3.h_S=h_S;
   
10. MS3.filprob=filprob;
    
11. 
    
12. switch_1/ntotal
    
13. switch_2/ntotal
    
14. switch_3/ntotal
    
15. indexSwitch
    
16. 
    
17. figure(1)
    
18. subplot (3,2,1); plot(beta_')
    
19. title ('convergence of omega');
    
20. subplot(3,2,2); plot ([a1_(1,:)' a2_(1,:)' a3_(1,:)']);
    
21. title('convergence of alpha');
    
22. subplot(323); plot([a1_(2,:)' a2_(2,:)' a3_(2,:)'])
    
23. title('convergence of beta1')
    
24. subplot(324); plot([a1_(3,:)' a2_(3,:)'  a3_(3,:)'])
    
25. title('convergence of beta2')
    
26. subplot(325); plot([a1_(4,:)' a2_(4,:)'  a3_(4,:)'])
    
27. title('convergence of beta3')
    
28. 
    
29. %%%% calculate the conditional variance
    
30. 
    
31. ha = h_S.*[(S==1) (S==2) (S==3)];
    
32. ha = ha(:,1) + ha(:,2)+ ha(:,3);
    
33.    
    
34. figure(2)
    
35. plot(ha)
    
36. 
    
37. 
    
38. 
    
39. %------------------FUNCTION-------------------------------------
    
40. function y = norm_rnd(sig);
    
41. % PURPOSE: random multivariate random vector based on
    
42. %          var-cov matrix sig
    
43. %---------------------------------------------------
    
44. % USAGE:   y = norm_rnd(sig)
    
45. % where:   sig = a square-symmetric covariance matrix
    
46. % NOTE: for mean b, var-cov sig use: b +  norm_rnd(sig)
    
47. %---------------------------------------------------      
    
48. % RETURNS: y = random vector normal draw mean 0, var-cov(sig)
    
49. %---------------------------------------------------
    
50. 
    
51. % written by:
    
52. % James P. LeSage, Dept of Economics
    
53. % University of Toledo
    
54. % 2801 W. Bancroft St,
    
55. % Toledo, OH 43606
    
56. % jlesage@spatial-econometrics.com
    
57. 
    
58. if nargin ~= 1
    
59. error('Wrong # of arguments to norm_rnd');
    
60. end;
    
61. 
    
62. h = chol(sig);
    
63. [nrow, ncol] = size(sig);
    
64. x = randn(nrow,1);
    
65. y = h'*x;

复制代码

不会也要顶一下

太深奥了

糊里糊涂

先mark一下  最近正在学MRS-GARCH模型  现在似然函数编出来的，却不知道该如何做参数估计，还好有楼主的帖子做个参考，实在感谢

你好，又要打扰楼主了，你发的程序我现在看了1/4多吧，但是又几个疑惑，是在等不急，就先问你一下吧：
首先是帖子第二楼的“start gibbs toop中，求了三个似然函数（09行——42行），ARCH效应中本期的条件方差是用前几期（这程序里是滞后阶数是3）的残差平方表示的，但是在程序中却是先把残差按波动的大小分为3个状态，而在3个状态中分别求条件方差h,那没在时间上，显然t期的h不是用t-1,t-2,t-3期的残差平方表示了。所以这里求h的方法让我无法理解，楼主能不能把原理讲一下，谢谢。
另一个问题仍是在GIBBS TOOP中，第49行有一个norm_rnd函数，这个是不是从正态分布中抽取随机数？因为我的matlab中现实没有这个函数，所以我就改成了normrnd（），但是又显示
Error using normrnd (line 28)
Requires at least two input arguments.

Error in code_ms_3s (line 147)
   beta_can = beta + normrnd(beta_sig);
我不知道问题出在哪里，还是希望楼主能指点一下，谢谢

49行可不可以改成
beta_can=beta+[normrnd(0,(0.003/4)^2);normrnd(0,(0.013/4)^2);normrnd(0,0.018/4)^2)]

syasd 发表于 2012-6-9 15:13
你好，又要打扰楼主了，你发的程序我现在看了1/4多吧，但是又几个疑惑，是在等不急，就先问你一下吧：
首先 ...

不好意思，过去几个月在工作。很少来论坛交流，也有7个月的时间做客户营销，没时间运行和修改matlab程序；
你问的问你很深入，也比较难以解决。毕竟MR-GARCH不是简单的。Griddy gibbs 抽样技术，也许你们自己有书，我从论坛找到国外的一个作者的贝叶斯金融统计建模，里面也提到过这个。不过很多我还没理解。
关于你提到的 normrnd()函数，这个是出自 lesage 的 JPL7 计量经济学工具箱里，是他本人自己写的一个多元正态分布随机函数，你可以下载JPL7，将这里的代码放进去，就可以正常运行。
如果你改用系统的多元正态分布随机数函数，因为输入参数和输出，也就是我们说的函数调用不匹配问题，导致出错。

明天要考试。所以不深入聊下去。我很高兴你这样深入的看这个程序和提这么多问题。有时间的时候多交流。
我的Q：280201722 潇湘品茗 30岁，大学毕业5年了。