bvar_density

 function bvar_density(maxnlags)
 computes the density of a bayesian var

 INPUTS
    maxnlags:      maximum number of lags in the bvar

 OUTPUTS
    none

 SPECIAL REQUIREMENTS
    none

0001 function bvar_density(maxnlags)
0002 % function bvar_density(maxnlags)
0003 % computes the density of a bayesian var
0004 %
0005 % INPUTS
0006 %    maxnlags:      maximum number of lags in the bvar
0007 %
0008 % OUTPUTS
0009 %    none
0010 %
0011 % SPECIAL REQUIREMENTS
0012 %    none
0013 
0014 % Copyright (C) 2003-2007 Christopher Sims
0015 % Copyright (C) 2007-2009 Dynare Team
0016 %
0017 % This file is part of Dynare.
0018 %
0019 % Dynare is free software: you can redistribute it and/or modify
0020 % it under the terms of the GNU General Public License as published by
0021 % the Free Software Foundation, either version 3 of the License, or
0022 % (at your option) any later version.
0023 %
0024 % Dynare is distributed in the hope that it will be useful,
0025 % but WITHOUT ANY WARRANTY; without even the implied warranty of
0026 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0027 % GNU General Public License for more details.
0028 %
0029 % You should have received a copy of the GNU General Public License
0030 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
0031 
0032 for nlags = 1:maxnlags
0033     [ny, nx, posterior, prior] = bvar_toolbox(nlags);
0034     
0035     posterior_int = matrictint(posterior.S, posterior.df, posterior.XXi);
0036     prior_int = matrictint(prior.S, prior.df, prior.XXi);
0037     
0038     lik_nobs = posterior.df - prior.df;
0039     
0040     log_dnsty = posterior_int - prior_int - 0.5*ny*lik_nobs*log(2*pi);
0041     
0042     disp(' ')
0043     fprintf('The marginal log density of the BVAR(%g) model is equal to %10.4f\n', ...
0044             nlags, log_dnsty);
0045     disp(' ')
0046 end
0047 
0048 
0049 function w = matrictint(S, df, XXi)
0050 % Computes the log of the integral of the kernel of the PDF of a
0051 % normal-inverse-Wishart distribution.
0052 %
0053 % S:   parameter of inverse-Wishart distribution
0054 % df:  number of degrees of freedom of inverse-Wishart distribution
0055 % XXi: first component of VCV matrix of matrix-normal distribution
0056 %
0057 % Computes the integral over (Phi, Sigma) of:
0058 %
0059 % det(Sigma)^(-k/2)*exp(-0.5*Tr((Phi-PhiHat)'*(XXi)^(-1)*(Phi-PhiHat)*Sigma^(-1)))*
0060 % det(Sigma)^((df+ny+1)/2)*exp(-0.5*Tr(Sigma^(-1)*S))
0061 %
0062 % (where k is the dimension of XXi and ny is the dimension of S and
0063 % Sigma)
0064 
0065 % Original file downloaded from:
0066 % http://sims.princeton.edu/yftp/VARtools/matlab/matrictint.m
0067 
0068 k=size(XXi,1);
0069 ny=size(S,1);
0070 [cx,p]=chol(XXi);
0071 [cs,q]=chol(S);
0072 
0073 if any(diag(cx)<100*eps)
0074     error('singular XXi')
0075 end
0076 if any(diag(cs<100*eps))
0077     error('singular S')
0078 end
0079 
0080 % Matrix-normal component
0081 w1 = 0.5*k*ny*log(2*pi)+ny*sum(log(diag(cx)));
0082 
0083 % Inverse-Wishart component
0084 w2 = -df*sum(log(diag(cs))) + 0.5*df*ny*log(2) + ny*(ny-1)*0.25*log(pi) + ggammaln(ny, df);
0085 
0086 w = w1 + w2;
0087 
0088 function lgg = ggammaln(m, df)
0089 if df <= (m-1)
0090     error('too few df in ggammaln')
0091 else
0092     garg = 0.5*(df+(0:-1:1-m));
0093     lgg = sum(gammaln(garg));
0094 end