% returns a matrix (N,number of shocks) of multinormal random variables 
% with mean mu and variance Sigma.
% if S=1, shocks are standardized
% [z]=multinormrnd(mu,Sigma,N,S)
% nsh=length(mu);
% R = chol(Sigma);
% z = repmat(mu(:)',N,1) + randn(N,nsh)*R;

function [z]=multinormrnd(mu,Sigma,N,S);
% mu = [1 2];
% Sigma = [1 .5; .5 2]; 
if nargin<4;
    S=0;
end
nsh=length(mu);
R = chol(Sigma);
z =randn(nsh,N);
if S && N>1
   z=z-(repmat(mean(z,2),1,size(z,2)));
   z=z./repmat(std(z,[],2),1,size(z,2));
end
 z=repmat(mu(:),1,N) + R'*z;
end