clear;
addpath('functions')
nobs=996;
btrue=[0.95 0.1;
0.1 0.95;
-0.1 0;
0 -0.1;
-0.05 0;
0 -0.05;
0 0];
sigmatrue=[2 1;
1 2];
datatrue=zeros(nobs,2);
for j=4:nobs
datatrue(j,:)=[datatrue(j-1,:) datatrue(j-2,:) datatrue(j-3,:) 1]*btrue+randn(1,2)*chol(sigmatrue);
end
dataQ=zeros(nobs/3,1);
jj=1;
for j=1:3:nobs
tmp=datatrue(j:j+2,1);
dataQ(jj,:)=mean(tmp);
jj=jj+1;
end
dataM=datatrue(:,2);
dataN=[ nan(rows(dataQ),2) dataQ(:,1) ];
dataN=vecr(dataN);
data0=[ zeros(rows(dataQ),2) dataQ(:,1) ];
data0=vecr(data0);
dataX=repmat(dataQ(:,1),1,3);
dataX=vecr(dataX);
data=[dataX dataM];
dataid=[ dataN dataM];
dataid0=[ data0 dataM];
mid=isnan(dataid);
N=cols(data);
REPS=11000;
BURN=10500;
L=3;
Y=data;
X=prepare(data,L);
Y=Y(L+1:end,:);
X=X(L+1:end,:);
dataid0=dataid0(L+1:end,:);
dataM=dataM(L+1:end,:);
T=rows(X);
b0=X\Y;
e0=Y-X*b0;
sigma=eye(N);
lamdaP = 1;
tauP = 10*lamdaP;
epsilonP= 1;
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);
beta0=[];
for j=0:L-1
beta0=[beta0 Y(L-j,:)];
end
P00=eye(cols(beta0))*0.1;
gibbs1=1;
for gibbs=1:REPS
X0=[X;xd];
Y0=[Y;yd];
mstar=vec(X0\Y0);
vstar=kron(sigma,invpd(X0'*X0));
chck=-1;
while chck<0
varcoef=mstar+(randn(1,N*(N*L+1))*chol(vstar))';
ee=stability(varcoef,N,L);
if ee==0;
chck=1;
end
end
resids=Y0-X0*reshape(varcoef,N*L+1,N);
scaleS=(resids'*resids);
sigma=iwpQ(T,invpd(scaleS));
ns=cols(P00);
[F,MUx]=comp(varcoef,N,L,1);
Q=zeros(ns,ns);
Q(1:N,1:N)=sigma;
beta_tt=zeros(T,ns);
ptt=zeros(T,ns,ns);
beta11=beta0;
p11=P00;
for i=1:T
nanid=mid(i,1);
if nanid==1
H=[0 0 0 0 0 0;
0 1 0 0 0 0];
rr=zeros(1,N);
rr(1)=1e10;
R=diag(rr);
else
H=[1/3 0 1/3 0 1/3 0;
0 1 0 0 0 0];
rr=zeros(1,N);
R=diag(rr);
end
x=H;
beta10=MUx+beta11*F';
p10=F*p11*F'+Q;
yhat=(x*(beta10)')';
eta=dataid0(i,:)-yhat;
feta=(x*p10*x')+R;
K=(p10*x')*invpd(feta);
beta11=(beta10'+K*eta')';
p11=p10-K*(x*p10);
ptt(i,:,:)=p11;
beta_tt(i,:)=beta11;
end
beta2 = zeros(T,ns);
bm2=beta2;
jv=1:2;
jv1=[1 3 5];
wa=randn(T,ns);
i=T;
p00=squeeze(ptt(i,jv1,jv1));
beta2(i,:)=beta_tt(i,:);
beta2(i,jv1)=mvnrnd(beta_tt(i:i,jv1),p00,1);
q=Q(jv,jv);
mu=MUx(jv);
f=F(jv,:);
for i=T-1:-1:1
pt=squeeze(ptt(i,:,:));
bm=beta_tt(i:i,:)+(pt*f'*invpd(f*pt*f'+q)*(beta2(i+1:i+1,jv)-mu-beta_tt(i,:)*f')')';
pm=pt-pt*f'*invpd(f*pt*f'+q)*f*pt;
beta2(i,:)=bm;
beta2(i:i,jv1)=mvnrnd(bm(jv1),pm(jv1,jv1),1);
bm2(i,:)=bm;
end
out=beta2(:,1);
datax=[out dataM];
Y=datax;
X=prepare(Y,L);
Y=Y(L+1:end,:);
X=X(L+1:end,:);
disp(sprintf('Iteration Number= %s ', num2str(gibbs)));
if gibbs>=BURN
dmat(:,gibbs1)=out;
bmat(gibbs1,:)=varcoef;
smat(gibbs1,:,:)=sigma;
gibbs1=gibbs1+1;
end
end
figure(1)
tmp=prctile(dmat,[50 ],2);
plot(tmp,'r');hold on;
plot(datatrue(4:end,1),'k')
legend('Posterior Median','True Data');