From b619eef664aada54786c049ec4cd7ad295b344d9 Mon Sep 17 00:00:00 2001
From: sebastien <sebastien@ac1d8469-bf42-47a9-8791-bf33cf982152>
Date: Fri, 28 Sep 2007 13:43:21 +0000
Subject: [PATCH] v4 rand_matrix_normal.m: inverted order of last two arguments
for more consistency with matrix-normal notation
git-svn-id: https://www.dynare.org/svn/dynare/dynare_v4@1400 ac1d8469-bf42-47a9-8791-bf33cf982152
---
matlab/PosteriorIRF.m | 2 +-
matlab/bvar_forecast.m | 2 +-
matlab/rand_matrix_normal.m | 16 ++++++++++------
3 files changed, 12 insertions(+), 8 deletions(-)
diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m
index f5098b39a8..8923cf8018 100644
--- a/matlab/PosteriorIRF.m
+++ b/matlab/PosteriorIRF.m
@@ -187,7 +187,7 @@ while b<=B
SIGMA_inv_upper_chol);
% draw from the conditional posterior of PHI
PHI_draw = rand_matrix_normal(NumberOfLagsTimesNvobs,nvobs, PHI, ...
- chol(iXX)', chol(SIGMAu_draw)');
+ chol(SIGMAu_draw)', chol(iXX)');
Companion_matrix(1:nvobs,:) = transpose(PHI_draw);
% Check for stationarity
explosive_var = any(abs(eig(Companion_matrix))>1.000000001);
diff --git a/matlab/bvar_forecast.m b/matlab/bvar_forecast.m
index cdf0c5cdad..131c15a53d 100644
--- a/matlab/bvar_forecast.m
+++ b/matlab/bvar_forecast.m
@@ -33,7 +33,7 @@ function bvar_forecast(nlags)
% Matlab, so using transpose
Sigma_lower_chol = chol(Sigma)';
- Phi = rand_matrix_normal(k, ny, posterior.PhiHat, XXi_lower_chol, Sigma_lower_chol);
+ Phi = rand_matrix_normal(k, ny, posterior.PhiHat, Sigma_lower_chol, XXi_lower_chol);
% All the eigenvalues of the companion matrix have to be on or inside the unit circle
Companion_matrix(1:ny,:) = Phi(1:ny*nlags,:)';
diff --git a/matlab/rand_matrix_normal.m b/matlab/rand_matrix_normal.m
index f61f07f1f9..e20dc92b88 100644
--- a/matlab/rand_matrix_normal.m
+++ b/matlab/rand_matrix_normal.m
@@ -5,16 +5,20 @@ function B = rand_matrix_normal(n, p, M, Omega_lower_chol, Sigma_lower_chol)
% B = rand_matrix_normal(n, p, M, Omega_lower_chol, Sigma_lower_chol)
%
% Returns an n-by-p matrix drawn from a Matrix-normal distribution
+%
+% B ~ MN_n*p(M, Omega, Sigma)
+%
+% Equivalent to vec(B) ~ N(vec(Mu), kron(Omega, Sigma))
+%
% Same notations than: http://en.wikipedia.org/wiki/Matrix_normal_distribution
+%
% M is the mean, n-by-p matrix
-% Omega_lower_chol is n-by-n, lower Cholesky decomposition of Omega
+% Omega_lower_chol is p-by-p, lower Cholesky decomposition of Omega
% (Omega_lower_chol = chol(Omega, 'lower'))
-% Sigma_lower_chol is p-by-p, lower Cholesky decomposition of Sigma
+% Sigma_lower_chol is n-by-n, lower Cholesky decomposition of Sigma
% (Sigma_lower_chol = chol(Sigma, 'lower'))
-%
-% Equivalent to vec(B) ~ N(vec(Mu), kron(Sigma, Omega))
B1 = randn(n * p, 1);
- B2 = kron(Sigma_lower_chol, Omega_lower_chol) * B1;
+ B2 = kron(Omega_lower_chol, Sigma_lower_chol) * B1;
B3 = reshape(B2, n, p);
- B = B3 + M;
\ No newline at end of file
+ B = B3 + M;
--
GitLab