teaching

% Estimating Markov Switching Models by MCMC

%

% yt=St'mu+et

%

% THIS PROGRAM WAS WRITTEN FOR MATLAB BY

% B. HARK YOO

% Bank of Korea

% harkyoo@bok.or.kr

randn('state',2);

rand('state',2);

% Generating data

kbar=2;

n=200;

mu=cumsum(0.5*ones(kbar,1));

sig=0.1;sig2=sig^2;

pr00=0.98;

pr11=0.8;

tranp=[pr00 1-pr00

1-pr11 pr11];

keye=eye(kbar);

st=markov(tranp,keye(:,1),n);

y=st'*mu+randn(n,1)*sig;

%%%%%%%%%% MCMC %%%%%%%%%%%%

disp(['Starting time=' num2str(clock)])

begtime=clock;

intvl=2;

burn=200;

nn=200;

bnn=burn+nn;

disp(['intvl=' num2str(intvl)])

disp(['burn=' num2str(burn)])

disp(['nn=' num2str(nn)])

% Initial Values

kbar=2;

keye=eye(kbar);

tp=rand(kbar,kbar);

tpsum=sum(tp');

for i=1:kbar

tp(i,:)=tp(i,:)/tpsum(i);

end

tp=[0.9 0.1

0.1 0.9];

zt=markov(tp,keye(:,1),n);

x=zt';

k=size(x,2);

m=inv(x'*x)*x'*y;

s=0.1;s2=s^2;

theta=[m;s];

mcmc=theta';

mcmctp=tp(:)';

ztsum=zeros(kbar,n);

% Priors

mm0=zeros(k,1);vm0=diag(1*ones(k,1));ivm0=inv(vm0);

ig0=3;ig1=5;

u0=keye*8+ones(kbar,kbar);

for ii=1:bnn

% Draw St by Multi-move Sampling

zt=mmove(y,x,theta,tp,kbar);

x=zt';

if ii>burn;

ztsum=ztsum+zt;

end

% Draw Transition Probabilities from Beta

np=nswitch(zt);

np=np+u0;

for i=1:kbar

if kbar==2;

p00=betarnd(np(i,1),np(i,2));

tp(i,:)=[p00 1-p00];

elseif kbar==3;

ns=np(i,(kbar-2):kbar);

p=betagen(ns');

tp(i,:)=p';

elseif kbar>=4;

j=4;

ns=np(i,(kbar-2):kbar);

p=betagen(ns');

for j=4:kbar

ns=np(i, (kbar-j+1):kbar);

p=betatran(ns',p);

end

tp(i,:)=p';

end

end

% Draw mu

ystar=y./s;

xstar=x./s;

mvar=inv(xstar'*xstar+ivm0);

mhat=mvar*(xstar'*ystar+ivm0*mm0);

mtemp=mhat+chol(mvar)'*randn(size(m,1),1);

if (mtemp(2:size(m,1))-mtemp(1:size(m,1)-1))>0;

m=mtemp;

end

% Draw s

e=y-x*m;e2=e.^2;

ss2=gamrnd(n/2+ig0,inv(1/ig1+1/2*sum(e2)));

s2=1/ss2;

s=sqrt(s2);

theta=[m;s];

mcmc=[mcmc

theta'];

mcmctp=[mcmctp

tp(:)'];

end

mcmc=[mcmc mcmctp];

runtime=etime(clock,begtime)/60;

disp(['Running time(in minute)=' num2str(runtime)])

% MCMC Results

[pos_m,pos_me,pos_sd,autocorr,hpi]=pstats(mcmc,intvl,burn,0.05);

disp('mean hpdi autocorr')

disp([ num2str([pos_m hpi autocorr])])

figure(1)

title('MCMC Draws')

for i=1:7;

subplot(4,2,i)

plot(mcmc(:,i))

end

pzt=ztsum/nn;

figure(2)

title('st and zt')

plot([st(2,:)' pzt(2,:)'])