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