function varmcmc = postdraw(vout,n,nosigprior)
% function varmcmc = postdraw(vout,n,nosigprior)
%  Copyright Christopher Sims, 2007
%  This program is distributed under the GPL3 license.  See <http://www.gnu.org/licenses/>.
% Provides n draws from the joint posterior on residual covariance and coefficients of a
% reduced form VAR.
% vout:	return value from rfvar3()
% n:		Number of draws
% nosigprior:	If true, don't use the Jeffreys-motivated improper prior.
% --------------------
  xxi = chol(vout.xxi);
  % ncf = size(xxi)(1); % this is Octave dialect
  ncf = size(xxi,1);
  S = vout.u' * vout.u;
  % df = size(vout.u)(1) - ncf;     % This effectively uses a |Sigma|^{-(nvar+1)/2} prior;
  df = size(vout.u,1) - ncf;
  %neq = size(vout.u)(2);
  neq = size(vout.u,2);
  lags = size(vout.By,3); %(ncf - size(vout.Bx)(2)) / neq;
  if nosigprior 
     df = df-size(S,1)-1;          % This undoes the |Sigma|^{-(nvar+1)/2} prior; 
  endif	
  wmat = rwwish(df,inv(S),n);
  for (it = 1:n) 
     wmat( :, :, it) = inv(wmat( :, :, it));
  endfor
  nmat = randn(ncf, neq, n);
  % cfmat = t(cbind(matrix(vout$By,neq,neq*lags),vout$Bx))
  cfmat = [reshape(vout.By, [neq, neq * lags]), vout.Bx]';
  for (ir = 1:n)
    nmat( :, :, ir) = xxi' * nmat(:, :, ir) * wmat(:, :, ir) + cfmat;
  endfor
  Byx = permute(nmat, [2,1,3]);
  By = Byx( :, 1:(neq*lags), :);
  By = reshape(By, [neq,neq,lags,n]);
  % Bx = as.vector(vout$Bx)+aperm(nmat,c(2,1,3))[,(neq*lags+1):ncf,]  # Bx added in both here
  % and in cfmat. Bug caught by A.Zawadwoski;
  Bx = Byx(:, (neq*lags+1):ncf, :);
  % Note that wmat[,,i] is ready for use as input to impulsdtrf;
  varmcmc.By=By;
  varmcmc.Bx=Bx;
  varmcmc.smat=wmat;