diff --git a/matlab/distributions/inverse_gamma_specification.m b/matlab/distributions/inverse_gamma_specification.m new file mode 100644 index 0000000000000000000000000000000000000000..f2b8585955f849c1b6617e1472d195dcc4d9f31f --- /dev/null +++ b/matlab/distributions/inverse_gamma_specification.m @@ -0,0 +1,76 @@ +function [s,nu] = inverse_gamma_specification(mu,sigma,type) + +% function [s,nu] = inverse_gamma_specification(mu,sigma,type) +% Specification of the inverse Gamma function parameters +% X ~ IG(s,nu) +% +% INPUTS +% mu: expectation +% sigma: standard deviation +% type=1: inverse Gamma 1 +% type=2: inverse Gamma 2 + +% OUTPUTS +% s: shape parameter +% nu: scale parameter +% +% SPECIAL REQUIREMENTS +% none + +% Copyright (C) 2003-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is 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. +% +% Dynare 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. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + +sigma2 = sigma^2; +mu2 = mu^2; + +if type == 2; % Inverse Gamma 2 + nu = 2*(2+mu2/sigma2); + s = 2*mu*(1+mu2/sigma2); +elseif type == 1; % Inverse Gamma 1 + if sigma2 < Inf; + nu = sqrt(2*(2+mu2/sigma2)); + nu2 = 2*nu; + nu1 = 2; + err = 2*mu2*gamma(nu/2)^2-(sigma2+mu2)*(nu-2)*gamma((nu-1)/2)^2; + while abs(nu2-nu1) > 1e-12 + if err > 0 + nu1 = nu; + if nu < nu2 + nu = nu2; + else + nu = 2*nu; + nu2 = nu; + end + else + nu2 = nu; + end + nu = (nu1+nu2)/2; + err = 2*mu2*gamma(nu/2)^2-(sigma2+mu2)*(nu-2)*gamma((nu-1)/2)^2; + end + s = (sigma2+mu2)*(nu-2); + else; + nu = 2; + s = 2*mu2/pi; + end; +else; + s = -1; + nu = -1; +end; + +% 01/18/2004 MJ replaced fsolve with secant +% suppressed chck +% changed order of output parameters \ No newline at end of file diff --git a/matlab/distributions/multivariate_normal_pdf.m b/matlab/distributions/multivariate_normal_pdf.m new file mode 100644 index 0000000000000000000000000000000000000000..66565f19d37e5d9db079c2697e9dc874609d146c --- /dev/null +++ b/matlab/distributions/multivariate_normal_pdf.m @@ -0,0 +1,36 @@ +function density = multivariate_normal_pdf(X,Mean,Sigma_upper_chol,n); +% Evaluates the density of a multivariate gaussian, with expectation Mean +% and variance Sigma_upper_chol'*Sigma_upper_chol, at X. +% +% +% INPUTS +% +% X [double] 1*n vector +% Mean [double] 1*n vector, expectation of the multivariate random variable. +% Sigma_upper_chol [double] n*n matrix, upper triangular Cholesky decomposition of Sigma (the covariance matrix). +% n [integer] dimension. +% +% OUTPUTS +% density [double] density +% +% SPECIAL REQUIREMENTS + +% Copyright (C) 2003-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is 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. +% +% Dynare 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. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + density = (2*pi)^(-.5*n) * ... + prod(diag(Sigma_upper_chol))^(-1) * ... + exp(-.5*(X-Mean)*(Sigma_upper_chol\(transpose(Sigma_upper_chol)\transpose(X-Mean)))); \ No newline at end of file diff --git a/matlab/distributions/multivariate_student_pdf.m b/matlab/distributions/multivariate_student_pdf.m new file mode 100644 index 0000000000000000000000000000000000000000..1865b41aee5860266fa206dd089e3cadafe1239e --- /dev/null +++ b/matlab/distributions/multivariate_student_pdf.m @@ -0,0 +1,36 @@ +function density = multivariate_student_pdf(X,Mean,Sigma_upper_chol,df); +% Evaluates the density of a multivariate student, with expectation Mean, +% variance Sigma_upper_chol'*Sigma_upper_chol and degrees of freedom df, at X. +% +% INPUTS +% +% X [double] 1*n vector +% Mean [double] 1*n vector, expectation of the multivariate random variable. +% Sigma_upper_chol [double] n*n matrix, upper triangular Cholesky decomposition of Sigma (the "covariance matrix"). +% df [integer] degrees of freedom. +% +% OUTPUTS +% density [double] density. +% +% SPECIAL REQUIREMENTS + +% Copyright (C) 2003-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is 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. +% +% Dynare 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. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + nn = length(X); + t1 = gamma( .5*(nn+df) ) / ( gamma( .5*nn ) * (df*pi)^(.5*nn) ) ; + t2 = t1 / prod(diag(Sigma_upper_chol)) ; + density = t2 / ( 1 + (X-Mean)*(Sigma_upper_chol\(transpose(Sigma_upper_chol)\transpose(X-Mean)))/df )^(.5*(nn+df)); \ No newline at end of file diff --git a/matlab/distributions/rand_inverse_wishart.m b/matlab/distributions/rand_inverse_wishart.m new file mode 100644 index 0000000000000000000000000000000000000000..f6de993e00f43edb5c5d4b89c68de4c9083dffcc --- /dev/null +++ b/matlab/distributions/rand_inverse_wishart.m @@ -0,0 +1,53 @@ +function G = rand_inverse_wishart(m, v, H_inv_upper_chol) + +% function G = rand_inverse_wishart(m, v, H_inv_upper_chol) +% rand_inverse_wishart Pseudo random matrices drawn from an +% inverse Wishart distribution +% G = rand_inverse_wishart(m, v, H_inv_upper_chol) +% Returns an m-by-m matrix drawn from an inverse-Wishart distribution. +% +% INPUTS: +% m: dimension of G and H_inv_upper_chol. +% v: degrees of freedom, greater or equal than m. +% H_inv_chol: upper cholesky decomposition of the inverse of the +% matrix parameter. +% The upper cholesky of the inverse is requested here +% in order to avoid to recompute it at every random draw. +% H_inv_upper_chol = chol(inv(H)) +% OUTPUTS: +% G: G ~ IW(m, v, H) where H = inv(H_inv_upper_chol'*H_inv_upper_chol) +% or, equivalently, using the correspondence between Wishart and +% inverse-Wishart: inv(G) ~ W(m, v, S) where +% S = H_inv_upper_chol'*H_inv_upper_chol = inv(H) +% +% SPECIAL REQUIREMENT +% none +% + +% Copyright (C) 2003-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is 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. +% +% Dynare 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. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + + X = randn(v, m) * H_inv_upper_chol; + + + % At this point, X'*X is Wishart distributed + % G = inv(X'*X); + + % Rather compute inv(X'*X) using the SVD + [U,S,V] = svd(X, 0); + SSi = 1 ./ (diag(S) .^ 2); + G = (V .* repmat(SSi', m, 1)) * V'; \ No newline at end of file diff --git a/matlab/distributions/rand_matrix_normal.m b/matlab/distributions/rand_matrix_normal.m new file mode 100644 index 0000000000000000000000000000000000000000..32aa4bafd9ee5d505ce9618e69cc564e5a01d99f --- /dev/null +++ b/matlab/distributions/rand_matrix_normal.m @@ -0,0 +1,43 @@ +function B = rand_matrix_normal(n, p, M, Omega_lower_chol, Sigma_lower_chol) + +% function B = rand_matrix_normal(n, p, M, Omega_lower_chol, Sigma_lower_chol) +% Pseudo random matrices drawn from a matrix-normal distribution +% B ~ MN_n*p(M, Omega, Sigma) +% Equivalent to vec(B) ~ N(vec(Mu), kron(Omega, Sigma)) +% +% INPUTS +% n: row +% p: column +% M: (n*p) matrix, mean +% Omega_lower_chol: (p*p), lower Cholesky decomposition of Omega, +% (Omega_lower_chol = chol(Omega, 'lower')) +% Sigma_lower_chol: (n*n), lower Cholesky decomposition of Sigma, +% (Sigma_lower_chol = chol(Sigma, 'lower')) +% +% OUTPUTS +% B: (n*p) matrix drawn from a Matrix-normal distribution +% +% SPECIAL REQUIREMENTS +% Same notations than: http://en.wikipedia.org/wiki/Matrix_normal_distribution + +% Copyright (C) 2003-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is 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. +% +% Dynare 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. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + + B1 = randn(n * p, 1); + B2 = kron(Omega_lower_chol, Sigma_lower_chol) * B1; + B3 = reshape(B2, n, p); + B = B3 + M; diff --git a/matlab/distributions/rand_multivariate_normal.m b/matlab/distributions/rand_multivariate_normal.m new file mode 100644 index 0000000000000000000000000000000000000000..bd70c51dfd7533efd22b37ad5de12681be6e4e7b --- /dev/null +++ b/matlab/distributions/rand_multivariate_normal.m @@ -0,0 +1,34 @@ +function draw = rand_multivariate_normal(Mean,Sigma_upper_chol,n) +% Pseudo random draws from a multivariate normal distribution, +% \mathcal N_n(Mean,Sigma), with expectation Mean and variance Sigma. +% +% INPUTS +% +% Mean [double] 1*n vector, expectation of the multivariate random variable. +% Sigma_upper_chol [double] n*n matrix, upper triangular Cholesky decomposition of Sigma (the covariance matrix). +% n [integer] dimension. +% +% OUTPUTS +% draw [double] 1*n vector drawn from a multivariate normal distribution with expectation Mean and +% covariance Sigma +% +% SPECIAL REQUIREMENTS + +% Copyright (C) 2003-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is 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. +% +% Dynare 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. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + + draw = Mean + randn(1,n) * Sigma_upper_chol; diff --git a/matlab/distributions/rand_multivariate_student.m b/matlab/distributions/rand_multivariate_student.m new file mode 100644 index 0000000000000000000000000000000000000000..cc4fc0e75a3817f41b19928b5ff59292f386fe9d --- /dev/null +++ b/matlab/distributions/rand_multivariate_student.m @@ -0,0 +1,39 @@ +function draw = rand_multivariate_student(Mean,Sigma_upper_chol,df) +% Pseudo random draws from a multivariate student distribution, +% with expectation Mean, variance Sigma*df/(df-2) and degrees of freedom df>0. +% +% INPUTS +% +% Mean [double] 1*n vector, expectation of the multivariate random variable. +% Sigma_upper_chol [double] n*n matrix, upper triangular Cholesky decomposition of Sigma +% (the covariance matrix up to a factor df/(df-2)). +% df [integer] degrees of freedom. +% +% OUTPUTS +% draw [double] 1*n vector drawn from a multivariate normal distribution with expectation Mean and +% covariance Sigma. +% +% REMARK This is certainly not the most efficient way... +% +% NOTE See Zellner (appendix B.2, 1971) for a definition. +% + +% Copyright (C) 2003-2008 Dynare Team +% +% This file is part of Dynare. +% +% Dynare is 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. +% +% Dynare 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. +% +% You should have received a copy of the GNU General Public License +% along with Dynare. If not, see <http://www.gnu.org/licenses/>. + + n = length(Mean); + draw = Mean + randn(1,n) * Sigma_upper_chol * sqrt(df/sum(randn(df,1).^2));