diff --git a/matlab/PosteriorIRF.m b/matlab/PosteriorIRF.m
index f5098b39a8e60e887e39e40f27828c2407021bff..8923cf8018ba630babdc592ee950f025fde7abe0 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 cdf0c5cdade2e0507d2d9defc3b8b9730d41e2d3..131c15a53d2dde5d775b12ea592c8d580fbdbbf4 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 f61f07f1f9fc0667954cc500b6d28833ec8f9197..e20dc92b88c3df1bbc15042ed6d50b620e80c9cf 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;