1 clear;
2 addpath('functions');
3 %generate artificial data
4 T=500;
5 B1=0.2;
6 B2=0.9;
7 C1=1;
8 C2=-1;
9 S1=10;
10 S2=1;
11 P=[0.95 0.05;0.05 0.95];
12 Q=[0.97 0;0.03 1];
13 strue=zeros(T,2);
14 strue(1,1)=1;
15 strue=simS(strue,P);
16 vtrue=zeros(T,2);
17 vtrue(1,1)=1;
18 check=-1;
19 while check<0
20 vtrue=simS(vtrue,Q);
21 if sum(vtrue(:,1))>20
22     check=1;
23 end
24 end
25 e=randn(T,1);
26 Y=zeros(T,1);
27 X=zeros(T,1);
28 for i=2:T;
29     X(i,:)=Y(i-1,:);
30     if strue(i,1)==1
31         if vtrue(i,1)==1
32     Y(i)=[X(i,:) 1]*[B1 C1]'+e(i)*sqrt(S1);
33         elseif vtrue(i,2)==1
34
35     Y(i)=[X(i,:) 1]*[B1 C1]'+e(i)*sqrt(S2);
36         end
37     elseif strue(i,2)==1
38    if vtrue(i,1)==1
39     Y(i)=[X(i,:) 1]*[B2 C2]'+e(i)*sqrt(S1);
40     elseif vtrue(i,2)==1
41     Y(i)=[X(i,:) 1]*[B2 C2]'+e(i)*sqrt(S2);
42    end
43     end
44 end
45 %data
46 y=Y;
47 x=[X ones(T,1)];
48 %specify starting values
49 phi1=[0.5;1];   %regime 1 coefficients
50 phi2=[0.8;-1];   %regime 2 coefficients
51 sig1=3;       %regime 1 variance
52 sig2=1;       %regime 2 variance
53 p=0.95;
54 q=0.95;
55 px=0.98;
56 pmat=[p 1-q;1-p q];
57 qmat=[px 0; 1-px 1];
58 VV=zeros(T,1);
59 ncrit=5; %each regime should have ncrit obs
60 %set Priors
61 %coefficients
62 B0=zeros(2,1); %prior mean
63 Sigma0=eye(2); %prior variance
64 %variances
65 d0=0.1; %prior scale
66 v0=1;   %prior df
67 %transition probabilities
68 u00=25; %p00~D(u11,u22)
69 u01=5;
70 u11=25; %p11~D(u22,u21)
71 u10=5;
72 v00=25; %q00~D(v00,v01)
73 v01=5;
74 out1=[];  %save coefficients
75 out2=[];  %save variances
76 out3=[];  %save S
77 out4=[]; %save p
78 out5=[]; %save VV
79 REPS=10000;
80 BURN=5000;
81 igibbs=1;
82 count=1;
83 while count<REPS-BURN
84 %step 1: sample S[t]
85  %%%%%%%%%%%%%%%%Run Hamilton Filter%%%%%%%%%%%%%%%%
86    %unconditional probabilities
87 A = [(eye(2)-pmat);ones(1,2)];
88            EN=[0;0;1];
89            ett11= pinv(A'*A)*A'*EN;
90
91     lik=0;
92     filter=zeros(T,2);
93     for j=1:T
94         if VV(j)==0
95     iS1=1/sig1;
96     iS2=1/sig1;
97     dsig1=sqrt(sig1);
98     dsig2=sqrt(sig1);
99         else
100          iS1=1/sig2;
101     iS2=1/sig2;
102     dsig1=sqrt(sig2);
103     dsig2=sqrt(sig2);
104         end
105         em1=y(j)-x(j,:)*phi1;
106         em2=y(j)-x(j,:)*phi2;
107         neta1=(1/dsig1)*exp(-0.5*(em1*iS1*em1'));%F(Y\S=0)
108         neta2=(1/dsig2)*exp(-0.5*(em2*iS2*em2'));%F(Y\S=1)
109         %%%Prediction Step%%%%
110         ett10=pmat*ett11;
111         %%%%Update Step%%%%
112         ett11=ett10.*[neta1;neta2]; %joint density F(Y,S)
113         fit=sum(ett11);           %Marginal density F(Y)
114         ett11=(ett11)/fit;    %conditional density F(S\Y) the weights of the likelihood
115         filter(j,1:2)=ett11';      %save filter probability ett
116         lik=lik+log(fit);      %save log likelihood
117
118     end
119
120
121
122  check=-1;
123  while check<0
124    %backward recursion to sample from H(S[t]\S[t+1],y)
125    S=zeros(T,1);
126    %time T
127    p1=filter(T,1);
128    p2=filter(T,2);
129    p=p1/(p1+p2);
130    u=rand(1,1);
131    S(T,1)=(u>=p);
132
133    for t=T-1:-1:1
134    if S(t+1)==0
135 p00=pmat(1,1)*filter(t,1);
136 p01=pmat(1,2)*filter(t,2);
137 elseif S(t+1)==1
138 p00=pmat(2,1)*filter(t,1);
139 p01=pmat(2,2)*filter(t,2);
140    end
141   u=rand(1,1);
142   p=p00/(p00+p01);
143   if u<p
144       S(t)=0;
145   else
146       S(t)=1;
147   end
148    end
149
150 if sum(S==0)>=ncrit && sum(S==1)>=ncrit
151     check=1;
152 end
153  end
154
155
156
157
158  %step 1: sample S[t]
159  %%%%%%%%%%%%%%%%Run Hamilton Filter%%%%%%%%%%%%%%%%
160    %unconditional probabilities
161      ett11= [1;0]; %start in regime 0 by assumption
162     iS1=1/sig1;
163     iS2=1/sig2;
164     dsig1=sqrt(sig1);
165     dsig2=sqrt(sig2);
166
167
168     filterx=zeros(T,2);
169     filterx(1,:)=ett11';
170     for j=2:T
171         if S(j)==0
172         phi1x=phi1;
173         phi2x=phi1;
174         else
175         phi1x=phi2;
176         phi2x=phi2;
177         end
178
179         em1=y(j)-x(j,:)*phi1x;
180         em2=y(j)-x(j,:)*phi2x;
181         neta1=(1/dsig1)*exp(-0.5*(em1*iS1*em1'));%F(Y\S=0)
182         neta2=(1/dsig2)*exp(-0.5*(em2*iS2*em2'));%F(Y\S=1)
183         %%%Prediction Step%%%%
184         ett10=qmat*ett11;
185         %%%%Update Step%%%%
186         ett11=ett10.*[neta1;neta2]; %joint density F(Y,S)
187         fit=sum(ett11);           %Marginal density F(Y)
188         ett11=(ett11)/fit;    %conditional density F(S\Y) the weights of the likelihood
189         filterx(j,1:2)=ett11';      %save filter probability ett
190
191     end
192
193
194
195
196  checkx=-1;
197  while checkx<0
198    %backward recursion to sample from H(S[t]\S[t+1],y)
199    VV=zeros(T,1);
200     %time T
201    p1=filterx(T,1);
202    p2=filterx(T,2);
203    p=p1/(p1+p2);
204    u=rand(1,1);
205    VV(T,1)=(u>=p);
206    for t=T-1:-1:1
207    if VV(t+1)==0
208 p00=qmat(1,1)*filterx(t,1);
209 p01=qmat(1,2)*filterx(t,2);
210 elseif VV(t+1)==1
211 p00=qmat(2,1)*filterx(t,1);
212 p01=qmat(2,2)*filterx(t,2);
213    end
214   u=rand(1,1);
215   p=p00/(p00+p01);
216   if u<p
217       VV(t)=0;
218   else
219       VV(t)=1;
220   end
221    end
222
223 if sum(VV==0)>=ncrit && sum(VV==1)>=ncrit
224     checkx=1;
225 end
226  end
227
228
229
230  %step 3 sample the transition matrix P
231
232     tranmat=switchg(S+1,[1;2]); %calculate the number of regime switches
233     N00=tranmat(1,1); %S(t-1)=0 S(t)=0
234     N01=tranmat(1,2); %S(t-1)=0 S(t)=1
235     N10=tranmat(2,1); %S(t-1)=1 S(t)=0
236     N11=tranmat(2,2); %S(t-1)=1 S(t)=1
237     %draw from the dirichlet density
238     p0=drchrnd([N00+u00;N01+u01]);
239     p=p0(1,1); %p00
240     p0=drchrnd([N10+u10;N11+u11]);
241     q=p0(2,1); %p11
242     pmat=[p 1-q;1-p q]; %transition prob matrix
243
244     %step 4 sample the transition matrix Q
245
246     tranmatx=switchg(VV+1,[1;2]); %calculate the number of regime switches
247     N00x=tranmat(1,1); %V(t-1)=0 V(t)=0
248     N01x=tranmat(1,2); %V(t-1)=0 V(t)=1
249
250     %draw from the dirichlet density
251     p0=drchrnd([N00x+v00;N01x+v01]);
252     px=p0(1,1); %p00
253
254     qmat=[px 0;1-px 1]; %transition prob matrix
255
256     %step 3 sample beta
257     %calculate time series of sigma[t]
258     sigmat=(VV==0).*sig1+(VV==1).*sig2;
259     sigmat1=sigmat(S==0);
260     sigmat2=sigmat(S==1);
261     % Select data in regime 1
262     id=find(S==0);
263     y1=y(id)./sqrt(sigmat1); %remove heteroscedasticity
264     x1=x(id,:)./(repmat(sqrt(sigmat1),1,2));
265     M=inv(inv(Sigma0)+(x1'*x1))*(inv(Sigma0)*B0+x1'*y1);
266     V=inv(inv(Sigma0)+(x1'*x1));
267     phi1=M+(randn(1,2)*chol(V))';
268     %Select data in regime 2
269     id=find(S==1);
270     y2=y(id)./sqrt(sigmat2);
271     x2=x(id,:)./(repmat(sqrt(sigmat2),1,2));
272     M=inv(inv(Sigma0)+(x2'*x2))*(inv(Sigma0)*B0+x2'*y2);
273     V=inv(inv(Sigma0)+(x2'*x2));
274     phi2=M+(randn(1,2)*chol(V))';
275
276
277     %step 4 sample sigma
278     residuals=(y-x*phi1).*(S==0)+(y-x*phi2).*(S==1);
279     %residuals regime 1
280     e1=residuals(VV==0);
281     T1=v0+rows(e1);
282     D1=d0+e1'*e1;
283     %draw from IG
284    z0=randn(T1,1);
285     z0z0=z0'*z0;
286    sig1=D1/z0z0;
287    %residuals regime 2
288     e2=residuals(VV==1);
289     T2=v0+rows(e2);
290     D2=d0+e2'*e2;
291     %draw from IG
292    z0=randn(T2,1);
293     z0z0=z0'*z0;
294    sig2=D2/z0z0;
295
296
297
298    %save and impose regime identification
299    if igibbs>BURN
300        chck=phi1(2,1)>phi2(2,1); %constant bigger in regime 1
301        if chck
302            out1=[out1;([phi1' phi2'])];
303            out2=[out2;([sig1 sig2 ])];
304            out3=[out3;S'];
305            out4=[out4;[p q px]];
306            out5=[out5;VV'];
307            count=count+1;
308        end
309
310    end
311    igibbs=igibbs+1;
312      disp(sprintf(' Replication %s , %s Saved Draws %s. ', ...
313              num2str(igibbs), num2str(count) ));
314 end
315
316 figure(1)
317 subplot(6,2,1);
318 hist(out1(:,1),50);
319 vline(B1)
320 title('Coefficient regime 1');
321 axis tight
322 subplot(6,2,2);
323 hist(out1(:,3),50);
324 title('Coefficient regime 2');
325 vline(B2)
326 axis tight
327 subplot(6,2,3);
328 hist(out1(:,2),50);
329 vline(C1)
330 title('Intercept regime 1');
331 axis tight
332 subplot(6,2,4);
333 hist(out1(:,4),50);
334 title('Intercept regime 2');
335 vline(C2)
336 axis tight
337 subplot(6,2,5);
338 hist(out2(:,1),50);
339 vline(S1)
340 title('\sigma_{1}');
341 axis tight
342 subplot(6,2,6);
343 hist(out2(:,2),50);
344 vline(S2)
345 title('\sigma_{2}');
346 axis tight
347 subplot(6,2,7);
348 hist(out4(:,1),50);
349 title('P_{00}');
350 vline(P(1,1))
351 axis tight
352 subplot(6,2,8);
353 hist(out4(:,2),50);
354 title('p_{11}');
355 vline(P(2,2))
356 axis tight
357 subplot(6,2,9);
358 hist(out4(:,3),50);
359 title('q_{00}');
360 vline(Q(1,1))
361 axis tight
362 subplot(6,2,10)
363 temp=mean(out3,1);
364 plot(temp,'c','LineWidth',2);
365 hold on
366 plot(strue(:,2),'k','LineWidth',2)
367 title('Probability of Regime 1');
368 legend('Estimate','True')
369 axis tight
370 subplot(6,2,11)
371 temp=mean(out5,1);
372 plot(temp,'c','LineWidth',2);
373 hold on
374 plot(vtrue(:,2),'k','LineWidth',2)
375 title('Probability of Regime 1 (Variance)');
376 legend('Estimate','True')
377 axis tight