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
|
function [cT,vR,kdf] = fn_gibbsglb(Sbd,idmat0,nvar,fss)
% [cT,vR,kdf] = gibbsglb(Sbd,idmat0,nvar,fss)
% Global setup outside the Gibbs loop -- c.f. gibbsvar
% Ref.: D.F. Waggoner and T.A. Zha: "Does Normalization Matter for Inference?"
% See Note Forecast (2) pp. 44-51
%
% Sbd: cell(nvar,1). Sbd=diag(S(1), ..., S(m)). Already divided by 'fss' for the
% posterior of a0 or A0(:) in Waggoner and Zha when one has asymmetric prior.
% Note,"bd" stands for block diagonal.
% idmat0: identification matrix for A0 with asymmetric prior; column -- equation.
% nvar: rank of A0 or # of variables
% fss: effective sample size (in the exponential term) --
% # of observations + # of dummy observations
% (or nSample - lags + # of dummy observations)
%-------------
% 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)
% kdf: the polynomial power in the Gamma or 1d Wishart distribution, used in
% "gibbsvar.m"
%
%
% 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/>.
%
kdf = fss; %2; %fss;
cT = cell(nvar,1);
for k=1:nvar
cT{k} = chol(Sbd{k}); % upper triangular but lower triangular Choleski
end
Q = cell(nvar,1);
% Q{i}: nvar-by-nvar with rank nvar-q(i) with ith equation a and
% q(i) which is # of non-zero elements in a. Note: Q*a=0.
vR = cell(nvar,1);
for k=1:nvar
Q{k} = diag(idmat0(:,k)==0); % constructing Q{k}
%
[jnk1,d1,v1] = svd(Q{k}/cT{k});
d1max=max(diag(d1));
if d1max==0
Idxk=1:nvar;
else
Idxk = find(diag(d1)<eps*d1max);
end
lenk = length(find(idmat0(:,k)));
if ( length(Idxk)<lenk )
warning('Dimension of non-zero A0(:,k) is different from svd(Q*inv(T))')
disp('Press ctrl-c to abort')
pause
else
jnk1 = v1(:,Idxk);
vR{k} = jnk1(:,1:lenk); % orthonormal basis for T*R
end
end
|
