clear;
addpath('functions');
T=500;
B1=0.2;
B2=0.9;
C1=1;
C2=-1;
S1=3;
S2=1;
GAMMA0=-1;
GAMMA1=1;
LAMBDA0=10;
strue=zeros(T,1);
strue(1,1)=1;
Z=getar(0.9,T);
e=randn(T,1);
Y=zeros(T,1);
X=zeros(T,1);
SSTAR=zeros(T,1);
ptrue=zeros(T,1);
qtrue=zeros(T,1);
for i=2:T;
X(i,:)=Y(i-1,:);
SSTAR(i,:)=GAMMA0+Z(i,:)*LAMBDA0+GAMMA1*strue(i-1,1)+randn(1,1);
if SSTAR(i,:)>=0
strue(i,1)=1;
end
ptrue(i)=normcdf((-GAMMA0-Z(i,:)*LAMBDA0));
qtrue(i)=1-normcdf((-GAMMA0-Z(i,:)*LAMBDA0-GAMMA1));
if strue(i,1)==0
Y(i)=[X(i,:) 1]*[B1 C1]'+e(i)*sqrt(S1);
else
Y(i)=[X(i,:) 1]*[B2 C2]'+e(i)*sqrt(S2);
end
end
y=Y;
x=[X ones(T,1)];
z=Z;
phi1=[0.5;1];
phi2=[0.8;-1];
sig1=3;
sig2=1;
gamma=[-1 0 1]';
pp=repmat(0.95,T,1);
qq=repmat(0.95,T,1);
ncrit=10;
B0=zeros(2,1);
Sigma0=eye(2);
d0=0.1;
v0=1;
GAMMA00=zeros(3,1);
SGAMMA0=eye(3).*1000;
out1=[];
out2=[];
out3=[];
out4=[];
out5=[];
out6=[];
out7=[];
REPS=20000;
BURN=15000;
igibbs=1;
count=1;
while count<REPS-BURN
pmat=[pp(1) 1-qq(1);
1-pp(1) qq(1)];
A = [(eye(2)-pmat);ones(1,2)];
EN=[0;0;1];
ett11= pinv(A'*A)*A'*EN;
iS1=1/sig1;
iS2=1/sig2;
lik=0;
filter=zeros(T,2);
for j=1:T
pmat=[pp(j) 1-qq(j);
1-pp(j) qq(j)];
em1=y(j)-x(j,:)*phi1;
em2=y(j)-x(j,:)*phi2;
neta1=(1/sqrt(sig1))*exp(-0.5*(em1*iS1*em1'));
neta2=(1/sqrt(sig2))*exp(-0.5*(em2*iS2*em2'));
ett10=pmat*ett11;
ett11=ett10.*[neta1;neta2];
fit=sum(ett11);
ett11=(ett11)/fit;
filter(j,1:2)=ett11';
lik=lik+log(fit);
end
check=-1;
while check<0
S=zeros(T,1);
p1=filter(T,1);
p2=filter(T,2);
p=p1/(p1+p2);
u=rand(1,1);
S(T,1)=(u>=p);
for t=T-1:-1:1
pmat=[pp(t+1) 1-qq(t+1);
1-pp(t+1) qq(t+1)];
if S(t+1)==0
p00=pmat(1,1)*filter(t,1);
p01=pmat(1,2)*filter(t,2);
elseif S(t+1)==1
p00=pmat(2,1)*filter(t,1);
p01=pmat(2,2)*filter(t,2);
end
u=rand(1,1);
p=p00/(p00+p01);
if u<p
S(t)=0;
else
S(t)=1;
end
end
if sum(S==0)>=ncrit && sum(S==1)>=ncrit
check=1;
end
end
sstar=zeros(T,1);
Slag=lag0(S,1);
Slag(1)=Slag(2);
zall=[ones(T,1) z Slag];
mm=zall*gamma;
for t=1:T
if S(t)==1
sstar(t)= normlt_rnd(mm(t),1,0);
elseif S(t)==0;
sstar(t)= normrt_rnd(mm(t),1,0);
end
end
pp=normcdf(-zall(:,1:end-1)*gamma(1:end-1));
qq=1-normcdf(-zall(:,1:end-1)*gamma(1:end-1)-gamma(end));
yy=sstar;
xx=zall;
V=inv(inv(SGAMMA0)+(xx'*xx));
M=V*(inv(SGAMMA0)*GAMMA00+xx'*yy);
gamma=M+(randn(1,3)*chol(V))';
id=find(S==0);
y1=y(id);
x1=x(id,:);
M=inv(inv(Sigma0)+(1/sig1)*(x1'*x1))*(inv(Sigma0)*B0+(1/sig1)*x1'*y1);
V=inv(inv(Sigma0)+(1/sig1)*(x1'*x1));
phi1=M+(randn(1,2)*chol(V))';
id=find(S==1);
y2=y(id);
x2=x(id,:);
M=inv(inv(Sigma0)+(1/sig2)*(x2'*x2))*(inv(Sigma0)*B0+(1/sig2)*x2'*y2);
V=inv(inv(Sigma0)+(1/sig2)*(x2'*x2));
phi2=M+(randn(1,2)*chol(V))';
e1=y1-x1*phi1;
T1=v0+rows(e1);
D1=d0+e1'*e1;
z0=randn(T1,1);
z0z0=z0'*z0;
sig1=D1/z0z0;
e2=y2-x2*phi2;
T2=v0+rows(e2);
D2=d0+e2'*e2;
z0=randn(T2,1);
z0z0=z0'*z0;
sig2=D2/z0z0;
if igibbs>BURN
chck=phi1(2,1)>phi2(2,1);
if chck
out1=[out1;([phi1' phi2'])];
out2=[out2;([sig1 sig2 ])];
out3=[out3;S'];
out4=[out4;pp'];
out5=[out5;qq'];
out6=[out6;gamma'];
out7=[out7;sstar'];
count=count+1;
end
end
igibbs=igibbs+1;
disp(sprintf(' Replication %s , %s Saved Draws %s. ', ...
num2str(igibbs), num2str(count) ));
end
figure(1)
subplot(8,2,1);
hist(out1(:,1),50);
vline(B1)
title('Coefficient regime 1');
axis tight
subplot(8,2,2);
hist(out1(:,3),50);
title('Coefficient regime 2');
vline(B2)
axis tight
subplot(8,2,3);
hist(out1(:,2),50);
vline(C1)
title('Intercept regime 1');
axis tight
subplot(8,2,4);
hist(out1(:,4),50);
title('Intercept regime 2');
vline(C2)
axis tight
subplot(8,2,5);
hist(out2(:,1),50);
vline(S1)
title('\sigma_{1}');
axis tight
subplot(8,2,6);
hist(out2(:,2),50);
vline(S2)
title('\sigma_{2}');
axis tight
subplot(8,2,7);
hist(out6(:,1),50);
title('\gamma_0');
vline(GAMMA0)
axis tight
subplot(8,2,8);
hist(out6(:,2),50);
title('\lambda');
vline(LAMBDA0)
axis tight
subplot(8,2,9);
hist(out6(:,3),50);
title('\gamma_1');
vline(GAMMA1)
axis tight
subplot(8,2,[11 12])
temp=mean(out3,1);
plot(temp,'c','LineWidth',2);
hold on
plot(strue(:,1),'k','LineWidth',2)
title('Probability of Regime 1');
legend('Estimate','True')
axis tight
subplot(8,2,[13 14])
temp=mean(out4,1);
plot(temp,'c','LineWidth',2);
hold on
plot(ptrue(:,1),'k','LineWidth',2)
title('p_00');
legend('Estimate','True')
axis tight
subplot(8,2,[15 16])
temp=mean(out5,1);
plot(temp,'c','LineWidth',2);
hold on
plot(qtrue(:,1),'k','LineWidth',2)
title('p_11');
legend('Estimate','True')
axis tight