From 2074327cd979cd54fa9807b3bdb974b30fa180b9 Mon Sep 17 00:00:00 2001 From: sebastien Date: Thu, 5 Jul 2007 12:49:26 +0000 Subject: [PATCH] v4 bvar_forecast.m: no more use 'lower' option of chol(), to ensure compatibility with older Matlab versions git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1340 ac1d8469-bf42-47a9-8791-bf33cf982152 --- matlab/bvar_forecast.m | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/matlab/bvar_forecast.m b/matlab/bvar_forecast.m index a6d75d97a..6c52e5cb8 100644 --- a/matlab/bvar_forecast.m +++ b/matlab/bvar_forecast.m @@ -12,15 +12,20 @@ function bvar_forecast(nlags) sims_no_shock = NaN(options_.forecast, ny, options_.bvar_replic); sims_with_shocks = NaN(options_.forecast, ny, options_.bvar_replic); - S_inv_chol = chol(inv(posterior.S)); - XXi_lower_chol = chol(posterior.XXi, 'lower'); + S_inv_upper_chol = chol(inv(posterior.S)); + + % Option 'lower' of chol() not available in old versions of + % Matlab, so using transpose + XXi_lower_chol = chol(posterior.XXi)'; k = ny*nlags+nx; for d = 1:options_.bvar_replic - Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_chol); + Sigma = rand_inverse_wishart(ny, posterior.df, S_inv_upper_chol); - Sigma_lower_chol = chol(Sigma, 'lower'); + % Option 'lower' of chol() not available in old versions of + % Matlab, so using transpose + Sigma_lower_chol = chol(Sigma)'; Phi = rand_matrix_normal(k, ny, posterior.PhiHat, XXi_lower_chol, Sigma_lower_chol); @@ -93,31 +98,31 @@ function bvar_forecast(nlags) end -function G = rand_inverse_wishart(m, v, H_inv_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_chol) +% G = rand_inverse_wishart(m, v, H_inv_upper_chol) % % Returns an m-by-m matrix drawn from an inverse-Wishart distribution. % -% m: dimension of G and H_inv_chol. +% 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 +% 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_chol = chol(inv(H)) +% H_inv_upper_chol = chol(inv(H)) % % In other words: -% G ~ IW(m, v, H) where H = inv(H_inv_chol'*H_inv_chol) +% 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_chol'*H_inv_chol = inv(H) +% inv(G) ~ W(m, v, S) where S = H_inv_upper_chol'*H_inv_upper_chol = inv(H) X = NaN(v, m); for i = 1:v - X(i, :) = randn(1, m) * H_inv_chol; + X(i, :) = randn(1, m) * H_inv_upper_chol; end % At this point, X'*X is Wishart distributed