clear;
addpath('functions');
T=500;
N=2;
L=1;
B1=[0.2 -0.1 -1; 0.5 -0.1 -1];
B2=[0.5 0.1 1; 0.7 0.1 1];
S1=[3 -0.5;-0.5 3];
S2=[1 0.1;0.1 1];
P=[0.95 0.05;0.05 0.95];
strue=zeros(T,2);
strue(1,1)=1;
strue=simS(strue,P);
e=randn(T,N);
Y=zeros(T,N);
X=zeros(T,N*L+1);
for i=2:T;
X(i,:)=[Y(i-1,:) 1];
if strue(i,1)==1
Y(i,:)=X(i,:)*B1'+e(i,:)*chol(S1);
else
Y(i,:)=X(i,:)*B2'+e(i,:)*chol(S2);
end
end
y=Y;
x=X;
maxtrys=1000;
phiols=x\y;
phi1=vec(phiols);
phi2=vec(phiols);
phi10=phi1;
phi20=phi2;
sig1=eye(N)*3;
sig2=eye(N);
p=0.95;
q=0.95;
pmat=[p 1-q;1-p q];
ncrit=10;
lamdaP = 10;
tauP = 10*lamdaP;
epsilonP= 1/10000;
muP=mean(y)';
sigmaP=[];
deltaP=[];
e0=[];
for i=1:N
ytemp=y(:,i);
xtemp=[lag0(ytemp,1) ones(rows(ytemp),1)];
ytemp=ytemp(2:end,:);
xtemp=xtemp(2:end,:);
btemp=xtemp\ytemp;
etemp=ytemp-xtemp*btemp;
stemp=etemp'*etemp/rows(ytemp);
if abs(btemp(1))>1
btemp(1)=1;
end
deltaP=[deltaP;btemp(1)];
sigmaP=[sigmaP;stemp];
e0=[e0 etemp];
end
[yd,xd] = create_dummies(lamdaP,tauP,deltaP,epsilonP,L,muP,sigmaP,N);
u00=25;
u01=5;
u11=25;
u10=5;
out1=[];
out2=[];
out3=[];
out4=[];
REPS=10000;
BURN=5000;
igibbs=1;
count=1;
while count<REPS-BURN
A = [(eye(2)-pmat);ones(1,2)];
EN=[0;0;1];
ett11= pinv(A'*A)*A'*EN;
iS1=inv(sig1);
iS2=inv(sig2);
lik=0;
filter=zeros(T,2);
for j=1:T
em1=y(j,:)-x(j,:)*reshape(phi1,N*L+1,N);
em2=y(j,:)-x(j,:)*reshape(phi2,N*L+1,N);
neta1=(1/sqrt(det(sig1)))*exp(-0.5*(em1*iS1*em1'));
neta2=(1/sqrt(det(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
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
tranmat=switchg(S+1,[1;2]);
N00=tranmat(1,1);
N01=tranmat(1,2);
N10=tranmat(2,1);
N11=tranmat(2,2);
p0=drchrnd([N00+u00;N01+u01]);
p=p0(1,1);
p0=drchrnd([N10+u10;N11+u11]);
q=p0(2,1);
pmat=[p 1-q;1-p q];
id=find(S==0);
Y1=y(id,:);
X1=x(id,:);
Y0=[Y1;yd];
X0=[X1;xd];
mstar1=vec(X0\Y0);
xx=X0'*X0;
ixx1=xx\eye(cols(xx));
[ phi1,PROBLEM1] = getcoef( mstar1,sig1,ixx1,maxtrys,N,L );
if PROBLEM1
phi1=phi01;
else
phi01=phi1;
end
e=Y0-X0*reshape(phi1,N*L+1,N);
scale=e'*e;
sig1=iwpQ(rows(Y0),inv(scale));
id=find(S==1);
Y2=y(id,:);
X2=x(id,:);
Y0=[Y2;yd];
X0=[X2;xd];
mstar2=vec(X0\Y0);
xx=X0'*X0;
ixx2=xx\eye(cols(xx));
[ phi2,PROBLEM2] = getcoef( mstar2,sig2,ixx2,maxtrys,N,L );
if PROBLEM2
phi2=phi02;
else
phi02=phi2;
end
e=Y0-X0*reshape(phi2,N*L+1,N);
scale=e'*e;
sig2=iwpQ(rows(Y0),inv(scale));
if igibbs>BURN
chck=log(det(sig1))>log(det(sig2));
if chck
out1(count,:)=[phi1' phi2'];
out2(count,:,:)=[sig1 sig2 ];
out3=[out3;S'];
out4=[out4;[p q]];
count=count+1;
end
end
igibbs=igibbs+1;
disp(sprintf(' Replication %s , %s Saved Draws %s. ', ...
num2str(igibbs), num2str(count) ));
end
figure(1)
temp=mean(out3,1);
plot(temp,'c','LineWidth',2);
hold on
plot(strue(:,2),'k','LineWidth',2)
title('Probability of Regime 1');
legend('Estimate','True')
axis tight
figure(2)
tmp=[vec(B1');vec(B2')];
for j=1:rows(tmp);
subplot(6,2,j);
hist(out1(:,j))
vline(tmp(j));
title(strcat('Coefficient:',num2str(j)));
end