1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
|
function A0gbs = gibbsvar(A0gbs,cT,vR,nvar,fss,kdf,Indxcol)
% A0gbs = gibbsvar(A0gbs,cT,vR,nvar,fss,kdf)
% One-step Gibbs sampler for structural VARs -- simultaneous equations approach
% Ref.: D.F. Waggoner and T. Zha: "Does Normalization Matter for Inference?"
% See Note Forecast (2) pp. 44-51
%
% A0gbs: the last draw of A0 matrix
% cT{i}: nvar-by-nvar where T'*T=Sbd{i} which is kind of covariance martrix
% divided by fss already
% vR{i}: nvar-by-q{i} -- orthonormral basis for T*R, which is obtained through
% single value decomposition of Q*inv(T). See gibbsglb.m
% nvar: rank of A0 or # of variables
% fss: effective sample size == nSample (T)-lags+# of dummy observations
% kdf: polynomial power in the Gamma or 1d Wishart distribution
% Indxcol: a vector of random columns this Gibbs draws.
% When this input is not supplied, the Gibbs draws all columns
%------------------
% A0bgs: new draw of A0 matrix in this Gibbs step
%
% NOTE: Added Indxcol on 2/13/00 so that previous programs may not be compatible.
%
%
% Copyright (C) 1997-2012 Daniel Waggoner and Tao Zha
%
% This free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% It is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% If you did not received a copy of the GNU General Public License
% with this software, see <http://www.gnu.org/licenses/>.
%
if (nargin==6), Indxcol=[1:nvar]; end
%---------------- Local loop for Gibbs given last A0gbs ----------
%* uR{i}: nvar-by-q{i} -- orthonormal with first q(i)-1 vectors lies in the
% span(T*a(j)|j~=i)
%*** Constructing u(1),...,u(q{i}) at each Gibbs step
%
sw0 = zeros(nvar,1);
for k=Indxcol % given last A0gbs and general new A0bgs
X = cT{k}*A0gbs; % given the latest updated A0gbs
X(:,k) = 0; % want to find non-zero sw s.t., X'*sw=0
[jL,Ux] = lu(X');
jIx0 = min(find(abs(diag(Ux))<eps)); % if isempty(jIx0), then something is wrong here
%
sw0(jIx0+1:end) = 0;
sw0(jIx0) = 1;
jA = Ux(1:jIx0-1,1:jIx0-1);
jb = Ux(1:jIx0-1,jIx0);
jy = -jA\jb;
sw0(1:jIx0-1) = jy;
sw = sw0/sqrt(sum(sw0.^2));
%
lenk = length(vR{k}(1,:));
gkb = zeros(lenk,1); % greek beta's
uR = zeros(nvar,lenk);
sx = zeros(nvar,lenk);
sx(:,1) = vR{k}(:,1);
for ki = 1:lenk-1
wxv = [sw'*sx(:,ki);sw'*vR{k}(:,ki+1)]; % w'*x and w'*v(+1)
dwxv = sqrt(sum(wxv.^2));
if (dwxv<eps)
uR(:,ki)=sx(:,ki); sx(:,ki+1)=vR{k}(:,ki+1);
else
wxv = wxv/dwxv;
uR(:,ki) = wxv(1)*vR{k}(:,ki+1) - wxv(2)*sx(:,ki);
sx(:,ki+1) = wxv(2)*vR{k}(:,ki+1) + wxv(1)*sx(:,ki);
end
end
uR(:,lenk) = sx(:,lenk); % uR now constructed
%
%--------- Gibbs loop ----------
%*** draw independently beta's that combine uR to form a's (columns of A0)
jcon = sqrt(1/fss);
gkb(1:lenk-1) = jcon*randn(lenk-1,1);
%* gamma or 1-d Wishart draw
jnk = jcon*randn(kdf+1,1);
if rand(1)<0.5
gkb(lenk) = sqrt(jnk'*jnk);
else
gkb(lenk) = -sqrt(jnk'*jnk);
end
%
%*** form new a(i) - ith column of A0
A0gbs(:,k) = (cT{k}\uR)*gkb;
end
|
