clear
addpath('functions');
[ data0 junk ]=xlsread('\data\datain.xls');
[ junk names ]=xlsread('\data\names.xls');
index=xlsread('\data\index.xls');
dindex=index(:,1);
index=index(:,2);
data=[];
for i=1:cols(data0);
if dindex(i)==1
dat=log(data0(:,i));
dat=diff(dat)*100;
elseif dindex(i)==3
dat=diff(data0(:,i));
else
dat=data0(2:end,i);
end
data=[data dat];
end
data=standardise(data);
z=xlsread('\data\baserate.xls');
z=z(2:end);
z=standardise(z);
KK=3;
L=2;
N=KK+1;
NN=cols(data);
T=rows(data)
pmat=extract(data,KK);
beta0=[pmat(1,:) z(1) zeros(1,N)];
ns=cols(beta0);
P00=eye(ns);
rmat=ones(NN,1);
Sigma=eye(N);
reps=5000;
burn=4000;
mm=1;
for m=1:reps;
fload=[];
floadr=[];
error=[];
for i=1:NN
y=data(:,i);
if index(i)==0
x=pmat;
else
x=[pmat z];
end
M=inv(x'*x)*(x'*y);
V=rmat(i)*inv(x'*x);
ff=M+(randn(1,cols(x))*cholx(V))';
if index(i)==0;
fload=[fload;ff'];
floadr=[floadr;0];
else
fload=[fload;ff(1:end-1)'];
floadr=[floadr;ff(end)];
end
error=[error y-x*ff];
end
fload(1:KK,1:KK)=eye(KK);
floadr(24:24+KK-1,1)=zeros(KK,1);
rmat=[];
for i=1:NN
rmati= IG(0,0,error(:,i));
rmat=[rmat rmati];
end
Y=[pmat z];
X=[lag0(Y,1) lag0(Y,2) ones(rows(Y),1)];
Y=Y(2:end,:);
X=X(2:end,:);
M=vec(inv(X'*X)*(X'*Y));
V=kron(Sigma,inv(X'*X));
chck=-1;
while chck<0
beta=M+(randn(1,N*(N*L+1))*cholx(V))';
S=stability(beta,N,L);
if S==0
chck=10;
end
end
beta1=reshape(beta,N*L+1,N);
errorsv=Y-X*beta1;
scale=errorsv'*errorsv;
Sigma=iwpQ(T,inv(scale));
H=zeros(NN,(KK+1)*L);
H(1:rows(fload),1:KK+1)=[fload floadr];
H(rows(floadr)+1,KK+1)=1;
R=diag([rmat 0]);
MU=[beta1(end,:)';zeros(N*(L-1),1)]';
F=[beta1(1:N*L,:)';eye(N*(L-1),N*L)];
Q=zeros(rows(F),rows(F));
Q(1:N,1:N)=Sigma;
beta_tt=[];
ptt=zeros(T,ns,ns);
i=1;
x=H;
beta10=MU+beta0*F';
p10=F*P00*F'+Q;
yhat=(x*(beta10)')';
eta=[data(i,:) z(i,:)]-yhat;
feta=(x*p10*x')+R;
K=(p10*x')*inv(feta);
beta11=(beta10'+K*eta')';
p11=p10-K*(x*p10);
beta_tt=[beta_tt;beta11];
ptt(i,:,:)=p11;
for i=2:T
beta10=MU+beta11*F';
p10=F*p11*F'+Q;
yhat=(x*(beta10)')';
eta=[data(i,:) z(i,:)]-yhat;
feta=(x*p10*x')+R;
K=(p10*x')*inv(feta);
beta11=(beta10'+K*eta')';
p11=p10-K*(x*p10);
ptt(i,:,:)=p11;
beta_tt=[beta_tt;beta11];
end
beta2 = zeros(T,ns);
jv=1:3;
wa=randn(T,ns);
f=F(jv,:);
q=Q(jv,jv);
mu=MU(jv);
i=T;
p00=squeeze(ptt(i,jv,jv));
beta2(i,jv)=beta_tt(i:i,jv)+(wa(i:i,jv)*cholx(p00));
for i=T-1:-1:1
pt=squeeze(ptt(i,:,:));
bm=beta_tt(i:i,:)+(pt*f'*inv(f*pt*f'+q)*(beta2(i+1:i+1,jv)-mu-beta_tt(i,:)*f')')';
pm=pt-pt*f'*inv(f*pt*f'+q)*f*pt;
beta2(i:i,jv)=bm(jv)+(wa(i:i,jv)*cholx(pm(jv,jv)));
end
pmat=beta2(:,1:3);
if m>burn
A0=cholx(Sigma);
yhat=zeros(36,N);
vhat=zeros(36,N);
vhat(3,1:N)=[0 0 0 1];
for i=3:36
yhat(i,:)=[yhat(i-1,:) yhat(i-2,:) 1]*[beta1(1:N*L,:);zeros(1,N)]+vhat(i,:)*A0;
end
yhat1=yhat*H(:,1:KK+1)';
irfmat(mm,1:36,1:NN+1)=(yhat1);
mm=mm+1;
end
end
irf=prctile(irfmat,[50 16 84],1);
figure(1)
j=1
for i=1:size(irf,3)
subplot(4,10,j)
plotx1(squeeze(irf(:,:,i))');
title(strcat('\fontsize{8}', names(i)))
j=j+1
axis tight
end