1 clear;
2 addpath('functions')
3
4 nobs=996;
5 btrue=[0.95 0.1;
6 0.1 0.95;
7 -0.1 0;
8 0 -0.1;
9 -0.05 0;
10 0 -0.05;
11 0 0];
12
13 sigmatrue=[2 1;
14 1 2];
15
16 datatrue=zeros(nobs,2);
17 for j=4:nobs
18 datatrue(j,:)=[datatrue(j-1,:) datatrue(j-2,:) datatrue(j-3,:) 1]*btrue+randn(1,2)*chol(sigmatrue);
19 end
20
21 dataQ=zeros(nobs/3,1);
22 jj=1;
23 for j=1:3:nobs
24 tmp=datatrue(j:j+2,1);
25 dataQ(jj,:)=mean(tmp);
26 jj=jj+1;
27 end
28 dataM=datatrue(:,2);
29
30
31 dataN=[ nan(rows(dataQ),2) dataQ(:,1) ];
32 dataN=vecr(dataN);
33 data0=[ zeros(rows(dataQ),2) dataQ(:,1) ];
34 data0=vecr(data0);
35
36 dataX=repmat(dataQ(:,1),1,3);
37 dataX=vecr(dataX);
38 data=[dataX dataM];
39 dataid=[ dataN dataM];
40 dataid0=[ data0 dataM];
41 mid=isnan(dataid);
42 N=cols(data);
43 REPS=11000;
44 BURN=10500;
45 L=3;
46 Y=data;
47 X=prepare(data,L);
48 Y=Y(L+1:end,:);
49 X=X(L+1:end,:);
50 dataid0=dataid0(L+1:end,:);
51 dataM=dataM(L+1:end,:);
52 T=rows(X);
53
54 b0=X\Y;
55 e0=Y-X*b0;
56 sigma=eye(N);
57
58 lamdaP = 1;
59 tauP = 10*lamdaP;
60 epsilonP= 1;
61 muP=mean(Y)';
62 sigmaP=[];
63 deltaP=[];
64 e0=[];
65 for i=1:N
66 ytemp=Y(:,i);
67 xtemp=[lag0(ytemp,1) ones(rows(ytemp),1)];
68 ytemp=ytemp(2:end,:);
69 xtemp=xtemp(2:end,:);
70 btemp=xtemp\ytemp;
71 etemp=ytemp-xtemp*btemp;
72 stemp=etemp'*etemp/rows(ytemp);
73 if abs(btemp(1))>1
74 btemp(1)=1;
75 end
76 deltaP=[deltaP;btemp(1)];
77 sigmaP=[sigmaP;stemp];
78 e0=[e0 etemp];
79 end
80
81 [yd,xd] = create_dummies(lamdaP,tauP,deltaP,epsilonP,L,muP,sigmaP,N);
82
83 beta0=[];
84 for j=0:L-1
85 beta0=[beta0 Y(L-j,:)];
86 end
87 P00=eye(cols(beta0))*0.1;
88
89 gibbs1=1;
90 for gibbs=1:REPS
91
92 X0=[X;xd];
93 Y0=[Y;yd];
94 mstar=vec(X0\Y0);
95 vstar=kron(sigma,invpd(X0'*X0));
96 chck=-1;
97 while chck<0
98 varcoef=mstar+(randn(1,N*(N*L+1))*chol(vstar))';
99 ee=stability(varcoef,N,L);
100 if ee==0;
101 chck=1;
102 end
103 end
104
105 resids=Y0-X0*reshape(varcoef,N*L+1,N);
106 scaleS=(resids'*resids);
107 sigma=iwpQ(T,invpd(scaleS));
108
109 ns=cols(P00);
110 [F,MUx]=comp(varcoef,N,L,1);
111 Q=zeros(ns,ns);
112 Q(1:N,1:N)=sigma;
113
114 beta_tt=zeros(T,ns);
115 ptt=zeros(T,ns,ns);
116
117 beta11=beta0;
118 p11=P00;
119 for i=1:T
120 nanid=mid(i,1);
121 if nanid==1
122 H=[0 0 0 0 0 0;
123 0 1 0 0 0 0];
124
125 rr=zeros(1,N);
126 rr(1)=1e10;
127 R=diag(rr);
128 else
129 H=[1/3 0 1/3 0 1/3 0;
130 0 1 0 0 0 0];
131
132 rr=zeros(1,N);
133 R=diag(rr);
134
135
136 end
137
138 x=H;
139
140 beta10=MUx+beta11*F';
141 p10=F*p11*F'+Q;
142 yhat=(x*(beta10)')';
143 eta=dataid0(i,:)-yhat;
144 feta=(x*p10*x')+R;
145
146 K=(p10*x')*invpd(feta);
147 beta11=(beta10'+K*eta')';
148 p11=p10-K*(x*p10);
149 ptt(i,:,:)=p11;
150 beta_tt(i,:)=beta11;
151 end
152
153
154 beta2 = zeros(T,ns);
155 bm2=beta2;
156 jv=1:2;
157 jv1=[1 3 5];
158 wa=randn(T,ns);
159 i=T;
160 p00=squeeze(ptt(i,jv1,jv1));
161 beta2(i,:)=beta_tt(i,:);
162 beta2(i,jv1)=mvnrnd(beta_tt(i:i,jv1),p00,1);
163 q=Q(jv,jv);
164 mu=MUx(jv);
165 f=F(jv,:);
166
167 for i=T-1:-1:1
168
169 pt=squeeze(ptt(i,:,:));
170 bm=beta_tt(i:i,:)+(pt*f'*invpd(f*pt*f'+q)*(beta2(i+1:i+1,jv)-mu-beta_tt(i,:)*f')')';
171 pm=pt-pt*f'*invpd(f*pt*f'+q)*f*pt;
172 beta2(i,:)=bm;
173 beta2(i:i,jv1)=mvnrnd(bm(jv1),pm(jv1,jv1),1);
174 bm2(i,:)=bm;
175 end
176 out=beta2(:,1);
177
178 datax=[out dataM];
179 Y=datax;
180 X=prepare(Y,L);
181 Y=Y(L+1:end,:);
182 X=X(L+1:end,:);
183
184 disp(sprintf('Iteration Number= %s ', num2str(gibbs)));
185 if gibbs>=BURN
186 dmat(:,gibbs1)=out;
187 bmat(gibbs1,:)=varcoef;
188 smat(gibbs1,:,:)=sigma;
189 gibbs1=gibbs1+1;
190 end
191 end
192 figure(1)
193 tmp=prctile(dmat,[50 ],2);
194 plot(tmp,'r');hold on;
195 plot(datatrue(4:end,1),'k')
196 legend('Posterior Median','True Data');
197