1function bvar_density(maxnlags)
2% function bvar_density(maxnlags)
3% computes the density of a bayesian var
4%
5% INPUTS
6%    maxnlags:      maximum number of lags in the bvar
7%
8% OUTPUTS
9%    none
10%
11% SPECIAL REQUIREMENTS
12%    none
13
14% Copyright (C) 2003-2007 Christopher Sims
15% Copyright (C) 2007-2017 Dynare Team
16%
17% This file is part of Dynare.
18%
19% Dynare is free software: you can redistribute it and/or modify
20% it under the terms of the GNU General Public License as published by
21% the Free Software Foundation, either version 3 of the License, or
22% (at your option) any later version.
23%
24% Dynare is distributed in the hope that it will be useful,
25% but WITHOUT ANY WARRANTY; without even the implied warranty of
26% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27% GNU General Public License for more details.
28%
29% You should have received a copy of the GNU General Public License
30% along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
31
32global oo_
33
34oo_.bvar.log_marginal_data_density=NaN(maxnlags,1);
35
36for nlags = 1:maxnlags
37    [ny, nx, posterior, prior] = bvar_toolbox(nlags);
38    oo_.bvar.posterior{nlags}=posterior;
39    oo_.bvar.prior{nlags}=prior;
40
41    posterior_int = matrictint(posterior.S, posterior.df, posterior.XXi);
42    prior_int = matrictint(prior.S, prior.df, prior.XXi);
43
44    lik_nobs = posterior.df - prior.df;
45
46    log_dnsty = posterior_int - prior_int - 0.5*ny*lik_nobs*log(2*pi);
47
48    oo_.bvar.log_marginal_data_density(nlags)=log_dnsty;
49
50    skipline()
51    fprintf('The marginal log density of the BVAR(%g) model is equal to %10.4f\n', ...
52            nlags, log_dnsty);
53    skipline()
54end
55
56
57function w = matrictint(S, df, XXi)
58% Computes the log of the integral of the kernel of the PDF of a
59% normal-inverse-Wishart distribution.
60%
61% S:   parameter of inverse-Wishart distribution
62% df:  number of degrees of freedom of inverse-Wishart distribution
63% XXi: first component of VCV matrix of matrix-normal distribution
64%
65% Computes the integral over (Phi, Sigma) of:
66%
67% det(Sigma)^(-k/2)*exp(-0.5*Tr((Phi-PhiHat)'*(XXi)^(-1)*(Phi-PhiHat)*Sigma^(-1)))*
68% det(Sigma)^((df+ny+1)/2)*exp(-0.5*Tr(Sigma^(-1)*S))
69%
70% (where k is the dimension of XXi and ny is the dimension of S and
71% Sigma)
72
73% Original file downloaded from:
74% http://sims.princeton.edu/yftp/VARtools/matlab/matrictint.m
75
76k=size(XXi,1);
77ny=size(S,1);
78[cx,p]=chol(XXi);
79[cs,q]=chol(S);
80
81if any(diag(cx)<100*eps)
82    error('singular XXi')
83end
84if any(diag(cs<100*eps))
85    error('singular S')
86end
87
88% Matrix-normal component
89w1 = 0.5*k*ny*log(2*pi)+ny*sum(log(diag(cx)));
90
91% Inverse-Wishart component
92w2 = -df*sum(log(diag(cs))) + 0.5*df*ny*log(2) + ny*(ny-1)*0.25*log(pi) + ggammaln(ny, df);
93
94w = w1 + w2;
95
96function lgg = ggammaln(m, df)
97if df <= (m-1)
98    error('too few df in ggammaln')
99else
100    garg = 0.5*(df+(0:-1:1-m));
101    lgg = sum(gammaln(garg));
102end
103