1 clear;
2 addpath('functions')
3 %generate artificial data
4 nobs=996; %996 months 332 quarters
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  %assume first variable is subject to temporal aggregation
21  dataQ=zeros(nobs/3,1); %quarterly data Y
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); %monthly data X
29 %arrange data
30 %put missing observations
31 dataN=[  nan(rows(dataQ),2) dataQ(:,1) ];  %puts NANs for missing obs
32 dataN=vecr(dataN);
33 data0=[ zeros(rows(dataQ),2) dataQ(:,1) ];  %same as above but zeros for missing
34 data0=vecr(data0);
35 %initial value of data just repeated observations
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);  %id for missing obs
42 N=cols(data);
43 REPS=11000;
44 BURN=10500;
45 L=3;  %lags
46 Y=data;
47 X=prepare(data,L); %X=[Y(-1),Y(-2)...constant]
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 %initial values for VAR coefficients
54 b0=X\Y;  %ols
55 e0=Y-X*b0;
56 sigma=eye(N);
57 %priors for VAR coefficients (Banbura et.al)
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 %dummy data to implement priors see http://ideas.repec.org/p/ecb/ecbwps/20080966.html
81 [yd,xd] = create_dummies(lamdaP,tauP,deltaP,epsilonP,L,muP,sigmaP,N);
82 %Initial values for the Kalman filter B0/0
83 beta0=[];
84 for j=0:L-1
85     beta0=[beta0 Y(L-j,:)];
86 end
87 P00=eye(cols(beta0))*0.1;  %P[0/0]
88 % Gibbs sampler
89 gibbs1=1;
90 for gibbs=1:REPS
91 %step 1 Draw VAR coefficients
92 X0=[X;xd]; %add dummy obs
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))'; %draw but keep stable
99 ee=stability(varcoef,N,L);
100 if ee==0;
101     chck=1;
102 end
103 end
104 %step 2 Draw VAR covariance
105  resids=Y0-X0*reshape(varcoef,N*L+1,N);
106 scaleS=(resids'*resids);
107 sigma=iwpQ(T,invpd(scaleS)); %draw for inverse Wishart
108 %step 3 Carter Kohn algorithm  to draw monthly data
109 ns=cols(P00);
110 [F,MUx]=comp(varcoef,N,L,1); %companion form for coefficients
111 Q=zeros(ns,ns);
112 Q(1:N,1:N)=sigma; %companion form for covariance
113 %Carter and Kohn algorithm to draw the factor
114 beta_tt=zeros(T,ns);          %will hold the filtered state variable
115 ptt=zeros(T,ns,ns);    % will hold its variance
116 % %%%%%%%%%%%Step 6a run Kalman Filter
117 beta11=beta0;
118 p11=P00;
119 for i=1:T
120 nanid=mid(i,1); %checks if data on GDP is missing
121 if nanid==1 %missing
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;  %big variance so missing data ignored
127     R=diag(rr);
128 else  %valid  observation for first variable every 3rd month
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     %Prediction
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 %updating
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 % Backward recursion to calculate the mean and variance of the distribution of the state
153 %vector
154 beta2 = zeros(T,ns);   %this will hold the draw of the state variable
155 bm2=beta2;
156 jv=1:2; %index of non singular block
157 jv1=[1 3 5]; %state variables to draw, 3, 5 are lagged states
158 wa=randn(T,ns);
159 i=T;  %period 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);%beta_tt(i:i,jv1)+(wa(i:i,jv1)*cholx(p00));   %draw for beta in period t from N(beta_tt,ptt)
163 q=Q(jv,jv);
164 mu=MUx(jv);
165 f=F(jv,:);
166 %periods t-1..to .1
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);     %bm(jv1)+(wa(i:i,jv1)*cholx(pm(jv1,jv1)));
174 bm2(i,:)=bm;
175 end
176 out=beta2(:,1); %draw of monthly data
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