function BVARPosteriorSampling(DSGEModel,CurrINI,controls)
% BVARPosteriorSampling: Runs the Gibbs sampler for a Bayesian VAR model.
%
% USAGE:
%
%       BVARPosteriorSampling(DSGEModel,CurrINI,controls)
%
% REQUIRED INPUT:  DSGEModel (structure) whose fields contain all required information
%                            about the model, such as input files, variable and sample
%                            data. Actual observations are not stored in this file.
%
%                  CurrINI (structure), with initialization information.
%
%                  controls (structure), whose fields provide handles to all the UI
%                            controls in YADA's main dialog.
%
% REQUIRED OUTPUT: None currently.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2008 European Central Bank.
%
%                       First version: January 25, 2007.
%                        This version: December 19, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% LICENSE INFORMATION:
%
%      YADA is free software: you can redistribute it and/or modify
%      it under the terms of the GNU General Public License as published by
%      the Free Software Foundation, either version 3 of the License, or
%      (at your option) any later version.
%
%      This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty of
%      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%      GNU General Public License for more details.
%
%      You should have received a copy of the GNU General Public License
%      along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%      YADA is released under the GNU General Public License, Version 3,
%      29 June 2007 <http://www.gnu.org/licenses/>. The current release of
%      the program was last modified by the ECB on the "This version" date
%      above.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% CHANGELOG:
%
% * 26-1-2007:  Continued work on the function.
%
% * 30-1-2007:  Continued work on the function.
%
% * 31-1-2007:  Made sure that the posterior mode is always used to initialize the function.
%
% * 1-2-2007:   Continued work on the Gibbs sampler needed to estimate the marginal likelihood with
%               Chib's marginal likelihood identity approach.
%
% * 9-2-2007:   Made sure that the function reads sufficient data from the posterior mode estimation
%               output file that exactly the same data and prior is used.
%
% * 28-2-2007:  The function now takes into account that the control "controls.posterior.burnin" displays the number
%               of burn-in draws rather than the percentage relative to the full posterior sample.
%
% * 13-6-2007:  The calculation of the inputs for "logMultiNormalPDF" now take numerical issues with the log-determinant
%               and the inverse of a covariance matrix into account by an eigenvalue decomposition.
%
% * 13-11-2007: Updated the documentation.
%
% * 2-1-2008:   Changed code using "DSGEModel.InitializeKalman".
%
% * 10-1-2008:  Changed the last drawnow call to pause(0.02).
%
% * 2-5-2008:   Missing "KalmanFirstObservation" is now not updated. The clear function is now called directly after
%               the last time a parameter group has been used.
%
% * 23-5-2008:  Updated the documentation.
%
% * 28-5-2008:  Added a wait dialog when the Gibbs sampler is running without displaying a progress dialog.
%
% * 29-5-2008:  Added the call to the GetCurrentDate function.
%
% * 19-12-2008: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

SetLogText(strvcat(' ','Initializing Gibbs sampling of Bayesian VAR...',' * Setting up the data matrices'));
%
% Load posterior mode data
%
ModeFile = [DSGEModel.OutputDirectory '\mode\PosteriorMode-BVAR-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '.mat'];
ModeData = load(ModeFile);
%
% check for additional fields
%
if isfield(ModeData,'DSGEModel')==0;
   SetLogText(strvcat('   STATUS: Insufficient data in posterior mode estimation file...',' ','Posterior mode estimation of Bayesian VAR was aborted.'));
   txt = 'The posterior mode estimation output file does not include all the necessary data for running the Gibbs sampler. Please run the posterior mode estimation routine once more.';
   About(txt,'information','YADA - BVAR Mode Estimation Issue',200,500,CurrINI);
   return;
end;
%
% the prior
%
OverallTightness = ModeData.OverallTightness;
CrossEqTightness = ModeData.CrossEqTightness;
HarmonicLagDecay = ModeData.HarmonicLagDecay;
VarianceTightness = ModeData.VarianceTightness;
thetaPi = ModeData.thetaPi;
OmegaPi = ModeData.OmegaPi;
Amatrix = ModeData.Amatrix;
qDF = ModeData.qDF;
thetaPsi = ModeData.thetaPsi;
OmegaPsi = ModeData.OmegaPsi;
%
% the data
%
x = ModeData.x;
T = size(x,2);
X = ModeData.X;
d = ModeData.d;
dLag = ModeData.dLag;
%
% original data in the DSGEModel structure
%
DSGEModel.BVARY = ModeData.DSGEModel.BVARY;
DSGEModel.BVARX = ModeData.DSGEModel.BVARX;
DSGEModel.BVARLags = ModeData.DSGEModel.BVARLags;
DSGEModel.SubBeginYear = ModeData.DSGEModel.SubBeginYear;
DSGEModel.SubBeginPeriod = ModeData.DSGEModel.SubBeginPeriod;
DSGEModel.SubEndYear = ModeData.DSGEModel.SubEndYear;
DSGEModel.SubEndPeriod = ModeData.DSGEModel.SubEndPeriod;
DSGEModel.DataFrequency = ModeData.DSGEModel.DataFrequency;
%
% the posterior mode estimates
%
PsiInit = ModeData.Psi;
PiInit = ModeData.Pi;
OmegaInit = ModeData.Omega;
PsiMode = PsiInit;
PiMode = PiInit;
OmegaMode = OmegaInit;
PsiPartialInit = PsiInit;
OmegaPartialInit = OmegaInit;
%
SetLogText(strvcat('   STATUS: OK',' * Using posterior mode estimates for initial values'));
%
% Collect the number of lags and number of variables
%
p = length(DSGEModel.BVARY);
k = DSGEModel.BVARLags;
q = length(DSGEModel.BVARX);
%
% Determine the sample
%
[FirstPeriod,LastPeriod] = CreateSubSample(DSGEModel);
KalmanFirstObservation = DSGEModel.KalmanFirstObservation;
%
% Check if we should adjust this based on k, the lag order
%
if FirstPeriod+KalmanFirstObservation-1<k+1;
   KalmanFirstObservation = k+2-FirstPeriod;
end;
%
% sequential marginal likelihood start value and step length
%
SeqStartStr = get(controls.dsge.sequentialstart,'String');
DSGEModel.SequentialStartIterationValue = str2num(strrep(StringTrim(SeqStartStr(DSGEModel.SequentialStartIteration,:)),',',''));
SeqStepStr = get(controls.dsge.sequentialstep,'String');
DSGEModel.SequentialStepLengthValue = str2num(strrep(StringTrim(SeqStepStr(DSGEModel.SequentialStepLength,:)),',',''));
%
% actual sample that we will use
%
[NewStartYear,NewStartPeriod] = AdjustSampleStart(DSGEModel.SubBeginYear,DSGEModel.SubBeginPeriod,DSGEModel.DataFrequency,KalmanFirstObservation-1);
%
if DSGEModel.PriorType==1;
   %
   % Minnesota-style prior on Pi
   %
   SetLogText(strvcat('   PRIOR: Minnesota-style',['          - Overall tightness:        ' num2str(OverallTightness)],['          - Cross-equation tightness: ' num2str(CrossEqTightness)],['          - Harmonic lag decay:       ' num2str(HarmonicLagDecay)]));
elseif DSGEModel.PriorType==2;
   %
   % Normal conditional on covariance matrix on Pi
   %
   SetLogText(strvcat('   PRIOR: Normal conditional on covariance',['          - Overall tightness:        ' num2str(OverallTightness)],['          - Harmonic lag decay:       ' num2str(HarmonicLagDecay)]));
else;
   SetLogText('   PRIOR: Diffuse');
end;
SetLogText(strvcat('   STATUS: OK',' * Collecting prior data for covariance parameters'));
%
% Covariance matrix parameters
%
if DSGEModel.OmegaPriorType==1;
   %
   % diffuse prior
   %
   SetLogText(strvcat('   PRIOR: Diffuse','   STATUS: OK'));
else;
   %
   % inverted Wishart
   %
   SetLogText('   PRIOR: Inverse Wishart');
   if DSGEModel.WishartType==1;
      %
      % We use the ML estimate of Omega for the A matrix
      %
      SetLogText(strvcat('          - ML estimate of Omega for A',['          - Degrees of freedom:       ' num2str(qDF,'%0.0f')],'   STATUS: OK'));
   else;
      SetLogText(strvcat(['          - variance tightness:       ' num2str(VarianceTightness)],['          - Degrees of freedom:       ' num2str(qDF,'%0.0f')],'   STATUS: OK'));
   end;
end;
%
% We're OK. Collect data about the number of random samples, etc.
%
%
% 1. Number of posterior draws
%
DrawsStr = get(controls.posterior.draws,'String');
NumPostDraws = str2num(StringTrim(strrep(DrawsStr(DSGEModel.PosteriorDrawsValue,:),',','')));
%
% 2. Number of posterior sample batches to save per chain
%
BatchStr = get(controls.posterior.batch,'String');
NumBatches = str2num(StringTrim(strrep(BatchStr(DSGEModel.SampleBatchValue,:),',','')));
if NumBatches>NumPostDraws;
   NumBatches = NumPostDraws;
end;
NumBatchDraws = NumPostDraws/NumBatches;
%
% 3. Number of burnin draws per chain
%
BurninStr = get(controls.posterior.burnin,'String');
NumBurnin = str2num(strrep(StringTrim(BurninStr(DSGEModel.BurnInValue,:)),',',''));
if NumBurnin>=NumPostDraws;
   SetLogText(strvcat(' ','Too many burn-in draws relation to total draws','Posterior sampling was aborted.'));
   txt = ['The number of burn-in draws (' StringTrim(BurninStr(DSGEModel.BurnInValue,:)) ') leaves no posterior draws (' StringTrim(DrawsStr(DSGEModel.PosteriorDrawsValue,:)) ') for inference. Please select a lower value or increase the number of posterior draws.'];
   About(txt,'information','YADA - Burnin Draws Problem',180,500,CurrINI);
   return;
end;
%
% make sure the output directory is exists
%
stat = MakeDir(DSGEModel.OutputDirectory,'gibbs');
if stat~=1;
   SetLogText(strvcat(' ','Couldn''t create the directory:',['  ' DSGEModel.OutputDirectory '\gibbs'],'Posterior sampling was aborted.'));
   txt = ['YADA was for some reason unable to create the directory "' DSGEModel.OutputDirectory '\gibbs". Posterior sampling has therefore been aborted.'];
   About(txt,'information','YADA - Directory Creation Problem',180,500,CurrINI);
   return;
end;
%
% Check if we should keep old draws
%
StartChain = 1;
NBatchDraws = 1;
BatchIndex = (1:NumBatchDraws);
StartBatch = 1;
StartPostDraws = 1;
%
% Initialize the output
%
if q>0;
   PsiPostSample = zeros(NumPostDraws,p*q);
   PsiPartialPostSample = PsiPostSample;
   LogPsiDensity = zeros(NumPostDraws,1);
end;
PiPostSample = zeros(NumPostDraws,p*p*k);
OmegaPostSample = zeros(NumPostDraws,p*(p+1)/2);
OmegaPartialPostSample = OmegaPostSample;
LogPiDensity = zeros(NumPostDraws,1);
NState = [];
UState = [];
%
% show in log window
%
SetLogText(strvcat(' ',[' * Time:                              ' StringTrim(datestr(now,14)) ', ' GetCurrentDate],[' * Number of draws:                   ' num2str(NumPostDraws,'%0.0f')],[' * Number of sample batches:          ' num2str(NumBatches,'%0.0f')],[' * Number of burn in draws:           ' num2str(NumBurnin,'%0.0f')],[' * Number of parameters to estimate:  ' num2str(p*(q+(p*k)+((p+1)/2)),'%0.0f')],[' * Effective sample size:             ' num2str(T,'%0.0f')],[' * Estimation sample:                 ' NewStartYear ':' NewStartPeriod ' - ' DSGEModel.SubEndYear ':' DSGEModel.SubEndPeriod],[' * Number of endogenous variables:    ' num2str(p,'%0.0f')],[' * Number of exogenous variables:     ' num2str(q,'%0.0f')],[' * Number of lags:                    ' num2str(k,'%0.0f')],' '));
%
% check if there are old draws that we could collect
%
if DSGEModel.OverwriteDraws==0;
   %
   % check which files may already have been stored.
   %
   ExitLoop = 0;
   BatchIndex = (1:NumBatchDraws);
   for NBatch=1:NumBatches;
      PostSampleFile = [DSGEModel.OutputDirectory '\gibbs\BVAR-Gibbs-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NBatch,'%0.0f') '.' num2str(NumBatches,'%0.0f') '.' num2str(NumPostDraws,'%0.0f') '.mat'];
      if FileExist(PostSampleFile)==1;
         %
         % load data and add to matrices
         %
         GibbsData = load(PostSampleFile);
         OmegaPostSample(BatchIndex,:) = GibbsData.BatchOmegaPostSample;
         OmegaPartialPostSample(BatchIndex,:) = GibbsData.BatchOmegaPartialPostSample;
         PiPostSample(BatchIndex,:) = GibbsData.BatchPiPostSample;
         LogPiDensity(BatchIndex,:) = GibbsData.BatchLogPiDensity;
         if q>0;
            PsiPostSample(BatchIndex,:) = GibbsData.BatchPsiPostSample;
            PsiPartialPostSample(BatchIndex,:) = GibbsData.BatchPsiPartialPostSample;
            LogPsiDensity(BatchIndex,:) = GibbsData.BatchLogPsiDensity;
         end;
         NState = GibbsData.BatchNState;
         UState = GibbsData.BatchUState;
         %
         PiInit = reshape(GibbsData.BatchPiPostSample(length(BatchIndex),:),p,p*k);
         OmegaInit = reshape(DuplicationMatrix(p)*GibbsData.BatchOmegaPostSample(length(BatchIndex),:)',p,p);
         OmegaPartialInit = reshape(DuplicationMatrix(p)*GibbsData.BatchOmegaPartialPostSample(length(BatchIndex),:)',p,p);
         if q>0;
            PsiInit = reshape(GibbsData.BatchPsiPostSample(length(BatchIndex),:),p,q);
            PsiPartialInit = reshape(GibbsData.BatchPsiPartialPostSample(length(BatchIndex),:),p,q);
         end;
         %
         BatchIndex = BatchIndex+NumBatchDraws;
      else;
         ExitLoop = 1;
         break;
      end;
   end;
   if (ExitLoop==1)|(NBatch==NumBatches);
      StartBatch = NBatch;
      StartPostDraws = max(BatchIndex)-NumBatchDraws+1;
   end;
end;

if StartPostDraws>NumPostDraws;
   %
   % Everything has already been gone through
   %
   SetLogText('Posterior sampling via the Gibbs sampler has already been completed!');
   txt = 'The Gibbs sampling algorithm has already been run and stored to disk. If you wish to re-run the posterior sampling, please check mark the "Overwrite old draws" checkbox in the "DSGE Posterior Sampling" frame on the "Settings" tab.';
   About(txt,'information','YADA - Posterior Sampling',180,500,CurrINI);
   return;
end;
%
% check if there are draws we can use
%
if (DSGEModel.OverwriteDraws==1)&(DSGEModel.PosteriorDrawsValue>1);
   %
   % The same number of batches is required
   %
   NPDrawsValue = DSGEModel.PosteriorDrawsValue-1;
   PostSampleFile = '';
   while NPDrawsValue>=1;
      NPDraws = str2num(StringTrim(strrep(DrawsStr(NPDrawsValue,:),',','')));
      NBatches = NumBatches;
      if NBatches>NPDraws;
         NBatches = NPDraws;
      end;
      %
      % check if output exists
      %
      for NBatch=1:NBatches;
         PostSampleFileI = [DSGEModel.OutputDirectory '\gibbs\BVAR-Gibbs-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NBatch,'%0.0f') '.' num2str(NumBatches,'%0.0f') '.' num2str(NPDraws,'%0.0f') '.mat'];
         if FileExist(PostSampleFileI)==1;
            PostSampleFile = PostSampleFileI;
         else;
            NBatch = NBatch-1;
            break;
         end;
      end;
      %
      % Go to the next entry
      %
      if isempty(PostSampleFile)==1;
         NPDrawsValue = NPDrawsValue-1;
      else;
         break;
      end;
   end;
   %
   % check if the PostSampleFile is empty. if not, disp NBatch, NBatches, NPDraws
   %
   if isempty(PostSampleFile)==0;
      txt = ['YADA has located ' num2str(ceil(NPDraws*(NBatch/NBatches)),'%0.0f') ' previous posterior draws. Would you like to include these in your ' num2str(NumPostDraws,'%0.0f') ' draws?'];
      answer = Query(txt,'question',150,'Posterior Sampling',500,'no',CurrINI);
      if strcmp(lower(answer),'yes')==1;
         %
         % first we load the data
         %
         NTempBatchDraws = NPDraws/NBatches;
         %
         BatchInd = (1:NTempBatchDraws);
         %
         for i=1:NBatch;
            PostSampleFile = [DSGEModel.OutputDirectory '\gibbs\BVAR-Gibbs-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NBatch,'%0.0f') '.' num2str(NumBatches,'%0.0f') '.' num2str(NPDraws,'%0.0f') '.mat'];
            %
            GibbsData = load(PostSampleFile);
            OmegaPostSample(BatchInd,:) = GibbsData.BatchOmegaPostSample;
            OmegaPartialPostSample(BatchInd,:) = GibbsData.BatchOmegaPartialPostSample;
            PiPostSample(BatchInd,:) = GibbsData.BatchPiPostSample;
            LogPiDensity(BatchInd,:) = GibbsData.BatchLogPiDensity;
            if q>0;
               PsiPostSample(BatchInd,:) = GibbsData.BatchPsiPostSample;
               PsiPartialPostSample(BatchInd,:) = GibbsData.BatchPsiPartialPostSample;
               LogPsiDensity(BatchInd,:) = GibbsData.BatchLogPsiDensity;
            end;
            NState = GibbsData.BatchNState;
            UState = GibbsData.BatchUState;
            BatchInd = BatchInd+NTempBatchDraws;
         end;
         StartPostDraws = ceil(NPDraws*(NBatch/NBatches))+1;
         PiInit = reshape(PiPostSample(StartPostDraws-1,:),p,p*k);
         OmegaInit = reshape(DuplicationMatrix(p)*OmegaPostSample(StartPostDraws-1,:)',p,p);
         OmegaPartialInit = reshape(DuplicationMatrix(p)*OmegaPartialPostSample(StartPostDraws-1,:)',p,p);
         if q>0;
            PsiInit = reshape(PsiPostSample(StartPostDraws-1,:),p,q);
            PsiPartialInit = reshape(PsiPartialPostSample(StartPostDraws-1,:),p,q);
         end;
         %
         %
         % save to file if NumBatches>1
         %
         if NumBatches>1;
            BatchIndex = (1:NumBatchDraws);
            %
            % how many batches can we store
            %
            for i=1:NumBatches;
               %
               % check if we should store data to file
               %
               %
               BatchLogPost = PiPostSample(BatchIndex,1);
               if sum(BatchLogPost~=0)==NumBatchDraws;
                  %
                  % we have a full sample
                  %
                  PostSampleFile = [DSGEModel.OutputDirectory '\gibbs\BVAR-Gibbs-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(i,'%0.0f') '.' num2str(NumBatches,'%0.0f') '.' num2str(NumPostDraws,'%0.0f') '.mat'];
                  %
                  BatchOmegaPostSample = OmegaPostSample(BatchIndex,:);
                  BatchOmegaPartialPostSample = OmegaPartialPostSample(BatchIndex,:);
                  BatchPiPostSample = PiPostSample(BatchIndex,:);
                  BatchLogPiDensity = LogPiDensity(BatchIndex,:);
                  if q>0;
                     BatchPsiPostSample = PsiPostSample(BatchIndex,:);
                     BatchPsiPartialPostSample = PsiPartialPostSample(BatchIndex,:);
                     BatchLogPsiDensity = LogPsiDensity(BatchIndex,:);
                  else;
                     BatchPsiPostSample = [];
                     BatchPsiPartialPostSample = [];
                     BatchLogPsiDensity = [];
                  end;
                  %
                  BatchNState = NState;
                  BatchUState = UState;
                  SetLogText(['   Saving: ' GetFilename(PostSampleFile)]);
                  %
                  % save to file
                  %
                  save(PostSampleFile,'BatchOmegaPostSample','BatchPiPostSample','BatchPsiPostSample','BatchOmegaPartialPostSample','BatchPsiPartialPostSample','BatchLogPiDensity','BatchLogPsiDensity','BatchNState','BatchUState');
                  %
                  BatchIndex = BatchIndex+NumBatchDraws;
               else;
                  NBatchDraws = sum(BatchLogPost~=0)+1;
                  StartBatch = i;
                  break;
               end;
            end;
         end;
      end;
   end;
end;
%
% deal with the random number generator
%
if DSGEModel.RandomNumberValue==1;
   rand('state',0);
   randn('state',0);
   %
   % make sure the number of draws equals the previous number of draws
   %
   if (isempty(NState)==0)&(isempty(UState)==0);
      rand('state',UState);
      randn('state',NState);
   end;
else;
   %
   % we set the state based on the current time
   %
   rand('state',sum(100*clock));
   randn('state',sum(100*clock));
end;
%
%
%
abort = '';
NBatch = StartBatch;
if DSGEModel.PriorType==1;
   %
   % Minnesota-style prior on Pi
   %
   % Set up initial values.
   % NOTE: Lower case means that we use vectorized values
   %
   thetaBarPi = vec(PiInit);
   ThetaBarPi = PiInit;
   if q>0;
      thetaBarPsi = vec(PsiInit);
      ThetaBarPsi = PsiInit;
      ThetaPartialPsi = PsiPartialInit;
   else;
      thetaBarPsi = [];
      ThetaBarPsi = [];
      ThetaPartialPsi = [];
   end;
   OmegaBar = OmegaInit;
   OmegaPartial = OmegaPartialInit;
   invOmegaPi = inv(OmegaPi);
   invOmegaPsi = inv(OmegaPsi);
   df = T+qDF;
   %
   % compute the value of the log likelihood + the log prior at the posterior mode
   %
   LogPosterior = BVARLogPosteriorMinnesota(OmegaMode,PiMode,PsiMode,Amatrix,qDF,thetaPsi,OmegaPsi,thetaPi,OmegaPi,x,X,d,dLag);
   %
   % Calculate the log of the full conditional posterior of Omega at the mode
   %
   LogOmegaDensity = logInverseWishartPDF(OmegaMode,(df+p+1)*BVAROmegaMinnesota(PiMode,PsiMode,Amatrix,qDF,x,X,d,dLag),df);
   %
   % check if we should show a progress dialog
   %
   if DSGEModel.ShowProgress==1;
      ProgressStructure.title = ['Progress for ' num2str(NumPostDraws,'%0.0f') ' posterior draws'];
      ProgressStructure.facecolor = CurrINI.progress_facecolor;
      ProgressStructure.startfacecolor = CurrINI.progress_startfacecolor;
      ProgressStructure.edgecolor = CurrINI.progress_edgecolor;
      ProgressStructure.bgcolor = CurrINI.progress_bgcolor;
      ProgressStructure.stop = 0;
      ProgressStructure.clock = DSGEModel.ShowProgressClock;
      ProgressStructure.label = 'Number of draws:';
      %
      ProgressStructure.name = 'Gibbs Sampling of Bayesian VAR';
      ProgressStructure.CurrINI = CurrINI;
      mywait = ProgressDLG(0,ProgressStructure);
      set(mywait,'Color',get(0,'defaultuicontrolbackgroundcolor'));
      drawnow;
   else;
      txt = ['Please wait while YADA runs the Gibbs sampler, where total number of posterior draws is ' int2str(NumPostDraws) ', while the number of burn-in draws ' int2str(NumBurnin) '. The current chain was started at: ' StringTrim(datestr(now,14)) ' on ' GetCurrentDate '.'];
      WaitHandle = WaitDLG(txt,'information','Gibbs Sampling of Bayesian VAR',500,200,CurrINI,0);
      WaitControls = get(WaitHandle,'UserData');
   end;
   %
   tic;
   for NPostDraws=StartPostDraws:NumPostDraws;
      %
      % update the progress dialog
      %
      if DSGEModel.ShowProgress==1;
         abort = get(mywait,'UserData');
         if (strcmp(abort,'cancel')==1)|(strcmp(abort,'done')==1);
            break;
         else;
            if NPostDraws<=NumPostDraws;
               ProgressDLG([NPostDraws/NumPostDraws NPostDraws]);
            end;
         end;
      end;
      %
      % we start by drawing OmegaBar conditional on Psi and Pi
      %
      OmegaBar = InvWishartRndFcn((df+p+1)*BVAROmegaMinnesota(ThetaBarPi,ThetaBarPsi,Amatrix,qDF,x,X,d,dLag),df);
      %
      % draw OmegaPartial with fixed value of Pi
      %
      OmegaPartial = InvWishartRndFcn((df+p+1)*BVAROmegaMinnesota(PiMode,ThetaPartialPsi,Amatrix,qDF,x,X,d,dLag),df);
      %
      % Next we draw vec(Pi) conditional on Psi and Omega
      %
      [MeanThetaBarPi,SigmaBarPi] = BVARPiMeanMinnesota(OmegaBar,ThetaBarPsi,thetaPi,invOmegaPi,x,X,d,dLag);
      thetaBarPi = MultiNormalRndFcn(vec(MeanThetaBarPi),SigmaBarPi,1);
      ThetaBarPi = reshape(thetaBarPi,p,p*k);
      %
      % now we draw vec(Psi) conditional on Pi and Omega
      %
      if q>0;
         [MeanThetaBarPsi,SigmaBarPsi] = BVARPsiMean(OmegaBar,ThetaBarPi,thetaPsi,invOmegaPsi,x,X,d,dLag);
         thetaBarPsi = MultiNormalRndFcn(vec(MeanThetaBarPsi),SigmaBarPsi,1);
         ThetaBarPsi = reshape(thetaBarPsi,p,q);
         %
         % draw ThetaPartialPsi
         %
         [MeanThetaPartialPsi,SigmaPartialPsi] = BVARPsiMean(OmegaPartial,PiMode,thetaPsi,invOmegaPsi,x,X,d,dLag);
         thetaPartialPsi = MultiNormalRndFcn(vec(MeanThetaPartialPsi),SigmaPartialPsi,1);
         ThetaPartialPsi = reshape(thetaPartialPsi,p,q);
         %
         % compute the height of the log density of PsiMode given OmegaPartial and PiMode
         %
         [V,D] = eig(SigmaPartialPsi);
         LogPsiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PsiMode),vec(MeanThetaPartialPsi),sum(log(diag(D))),(V/D)*V');
%         LogPsiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PsiMode),vec(MeanThetaPartialPsi),log(det(SigmaPartialPsi)),SigmaPartialPsi\eye(p*q));
      end;
      %
      % compute the height of the log density of PiMode given ThetaBarPsi and OmegaBar
      %
      [MeanThetaBarPi,SigmaBarPi] = BVARPiMeanMinnesota(OmegaBar,ThetaBarPsi,thetaPi,invOmegaPi,x,X,d,dLag);
      [V,D] = eig(SigmaBarPi);
      LogPiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PiMode),vec(MeanThetaBarPi),sum(log(diag(D))),(V/D)*V');
%      LogPiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PiMode),vec(MeanThetaBarPi),log(det(SigmaBarPi)),SigmaBarPi\eye(p*p*k));
      %
      % Store draws
      %
      OmegaPostSample(NPostDraws,:) = vech(OmegaBar)';
      OmegaPartialPostSample(NPostDraws,:) = vech(OmegaPartial)';
      PiPostSample(NPostDraws,:) = thetaBarPi';
      if q>0;
         PsiPostSample(NPostDraws,:) = thetaBarPsi';
         PsiPartialPostSample(NPostDraws,:) = thetaPartialPsi';
      end;
      %
      % Check the batch counter
      %
      if NBatchDraws<NumBatchDraws;
         %
         % increase
         %
         NBatchDraws = NBatchDraws+1;
      else;
         %
         % save a batch
         %
         PostSampleFile = [DSGEModel.OutputDirectory '\gibbs\BVAR-Gibbs-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NBatch,'%0.0f') '.' num2str(NumBatches,'%0.0f') '.' num2str(NumPostDraws,'%0.0f') '.mat'];
         BatchOmegaPostSample = OmegaPostSample(BatchIndex,:);
         BatchOmegaPartialPostSample = OmegaPartialPostSample(BatchIndex,:);
         BatchPiPostSample = PiPostSample(BatchIndex,:);
         BatchLogPiDensity = LogPiDensity(BatchIndex,:);
         if q>0;
            BatchPsiPostSample = PsiPostSample(BatchIndex,:);
            BatchPsiPartialPostSample = PsiPartialPostSample(BatchIndex,:);
            BatchLogPsiDensity = LogPsiDensity(BatchIndex,:);
         else;
            BatchPsiPostSample = [];
            BatchPsiPartialPostSample = [];
            BatchLogPsiDensity = [];
         end;
         %
         BatchNState = randn('state');
         BatchUState = rand('state');
         SetLogText(['   Saving: ' GetFilename(PostSampleFile)]);
         save(PostSampleFile,'BatchOmegaPostSample','BatchPiPostSample','BatchPsiPostSample','BatchOmegaPartialPostSample','BatchPsiPartialPostSample','BatchLogPiDensity','BatchLogPsiDensity','BatchNState','BatchUState');
         %
         NBatch = NBatch+1;
         NBatchDraws = 1;
         BatchIndex = BatchIndex+NumBatchDraws;
      end;
      drawnow;
   end;
   GibbsDrawTime = toc;
   %
   % close the current waitbar
   %
   if DSGEModel.ShowProgress==1;
      if ishandle(mywait)==1;
         set(mywait,'UserData','done');
         close(mywait);
         drawnow;
         pause(0.02);
      end;
   else;
      if ishandle(WaitHandle)==1;
         set(WaitControls.text,'UserData','done');
         delete(WaitHandle);
         drawnow;
         pause(0.02);
      end;
   end;
elseif DSGEModel.PriorType==2;
   %
   % Normal conditional on covariance matrix on Pi
   %
   %
   % Set up initial values.
   % NOTE: Lower case means that we use vectorized values
   %
   ThetaBarPi = PiInit;
   ThetaPi = reshape(thetaPi,p,p*k);
   if q>0;
      thetaBarPsi = vec(PsiInit);
      ThetaBarPsi = PsiInit;
      ThetaPartialPsi = PsiPartialInit;
   else;
      thetaBarPsi = [];
      ThetaBarPsi = [];
      ThetaPartialPsi = [];
   end;
   OmegaBar = OmegaInit;
   OmegaPartial = OmegaPartialInit;
   invOmegaPi = inv(OmegaPi);
   invOmegaPsi = inv(OmegaPsi);
   df = T+qDF+(p*k);
   %
   % compute the value of the log likelihood + the log prior at the posterior mode
   %
   LogPosterior = BVARLogPosteriorNormalCond(OmegaMode,PiMode,PsiMode,Amatrix,qDF,thetaPsi,OmegaPsi,ThetaPi,OmegaPi,x,X,d,dLag);
   %
   % Calculate the log of the full conditional posterior of Omega at the mode
   %
   LogOmegaDensity = logInverseWishartPDF(OmegaMode,(df+p+1)*BVAROmegaNormal(PiMode,PsiMode,Amatrix,qDF,ThetaPi,invOmegaPi,x,X,d,dLag),df);
   %
   % check if we should show a progress dialog
   %
   if DSGEModel.ShowProgress==1;
      ProgressStructure.title = ['Progress for ' num2str(NumPostDraws,'%0.0f') ' posterior draws'];
      ProgressStructure.facecolor = CurrINI.progress_facecolor;
      ProgressStructure.startfacecolor = CurrINI.progress_startfacecolor;
      ProgressStructure.edgecolor = CurrINI.progress_edgecolor;
      ProgressStructure.bgcolor = CurrINI.progress_bgcolor;
      ProgressStructure.stop = 0;
      ProgressStructure.clock = DSGEModel.ShowProgressClock;
      ProgressStructure.label = 'Number of draws:';
      %
      ProgressStructure.name = 'Gibbs Sampling of Bayesian VAR';
      ProgressStructure.CurrINI = CurrINI;
      mywait = ProgressDLG(0,ProgressStructure);
      set(mywait,'Color',get(0,'defaultuicontrolbackgroundcolor'));
      drawnow;
   end;
   %
   tic;
   for NPostDraws=StartPostDraws:NumPostDraws;
      %
      % update the progress dialog
      %
      if DSGEModel.ShowProgress==1;
         abort = get(mywait,'UserData');
         if (strcmp(abort,'cancel')==1)|(strcmp(abort,'done')==1);
            break;
         else;
            if NPostDraws<=NumPostDraws;
               ProgressDLG([NPostDraws/NumPostDraws NPostDraws]);
            end;
         end;
      end;
      %
      % we start by drawing OmegaBar conditional on Psi and Pi
      %
      OmegaBar = InvWishartRndFcn((df+p+1)*BVAROmegaNormal(ThetaBarPi,ThetaBarPsi,Amatrix,qDF,ThetaPi,invOmegaPi,x,X,d,dLag),df);
      %
      % draw OmegaPartial with fixed value of Pi
      %
      OmegaPartial = InvWishartRndFcn((df+p+1)*BVAROmegaNormal(PiMode,ThetaPartialPsi,Amatrix,qDF,ThetaPi,invOmegaPi,x,X,d,dLag),df);
      %
      % Next we draw vec(Pi) conditional on Psi and Omega
      %
      [MeanThetaBarPi,OmegaBarPi] = BVARPiMeanNormalCond(ThetaBarPsi,ThetaPi,invOmegaPi,x,X,d,dLag);
      thetaBarPi = MultiNormalRndFcn(vec(MeanThetaBarPi),kron(OmegaBarPi,OmegaBar),1);
      ThetaBarPi = reshape(thetaBarPi,p,p*k);
      %
      % now we draw vec(Psi) conditional on Pi and Omega
      %
      if q>0;
         [MeanThetaBarPsi,SigmaBarPsi] = BVARPsiMean(OmegaBar,ThetaBarPi,thetaPsi,invOmegaPsi,x,X,d,dLag);
         thetaBarPsi = MultiNormalRndFcn(vec(MeanThetaBarPsi),SigmaBarPsi,1);
         ThetaBarPsi = reshape(thetaBarPsi,p,q);
         %
         % draw ThetaPartialPsi
         %
         [MeanThetaPartialPsi,SigmaPartialPsi] = BVARPsiMean(OmegaPartial,PiMode,thetaPsi,invOmegaPsi,x,X,d,dLag);
         thetaPartialPsi = MultiNormalRndFcn(vec(MeanThetaPartialPsi),SigmaPartialPsi,1);
         ThetaPartialPsi = reshape(thetaPartialPsi,p,q);
         %
         % compute the height of the log density of PsiMode given OmegaPartial and PiMode
         %
         [V,D] = eig(SigmaPartialPsi);
         LogPsiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PsiMode),vec(MeanThetaPartialPsi),sum(log(diag(D))),(V/D)*V');
%         LogPsiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PsiMode),vec(MeanThetaPartialPsi),log(det(SigmaPartialPsi)),SigmaPartialPsi\eye(p*q));
      end;
      %
      % compute the height of the log density of PiMode given ThetaBarPsi and OmegaBar
      %
      [MeanThetaBarPi,OmegaBarPi] = BVARPiMeanNormalCond(ThetaBarPsi,ThetaPi,invOmegaPi,x,X,d,dLag);
      SigmaBarPi = kron(OmegaBarPi,OmegaBar);
      [V,D] = eig(SigmaBarPi);
      LogPiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PiMode),vec(MeanThetaBarPi),sum(log(diag(D))),(V/D)*V');
%      LogPiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PiMode),vec(MeanThetaBarPi),log(det(SigmaBarPi)),SigmaBarPi\eye(p*p*k));
      %
      % Store draws
      %
      OmegaPostSample(NPostDraws,:) = vech(OmegaBar)';
      OmegaPartialPostSample(NPostDraws,:) = vech(OmegaPartial)';
      PiPostSample(NPostDraws,:) = thetaBarPi';
      if q>0;
         PsiPostSample(NPostDraws,:) = thetaBarPsi';
         PsiPartialPostSample(NPostDraws,:) = thetaPartialPsi';
      end;
      %
      % Check the batch counter
      %
      if NBatchDraws<NumBatchDraws;
         %
         % increase
         %
         NBatchDraws = NBatchDraws+1;
      else;
         %
         % save a batch
         %
         PostSampleFile = [DSGEModel.OutputDirectory '\gibbs\BVAR-Gibbs-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NBatch,'%0.0f') '.' num2str(NumBatches,'%0.0f') '.' num2str(NumPostDraws,'%0.0f') '.mat'];
         BatchOmegaPostSample = OmegaPostSample(BatchIndex,:);
         BatchOmegaPartialPostSample = OmegaPartialPostSample(BatchIndex,:);
         BatchPiPostSample = PiPostSample(BatchIndex,:);
         BatchLogPiDensity = LogPiDensity(BatchIndex,:);
         if q>0;
            BatchPsiPostSample = PsiPostSample(BatchIndex,:);
            BatchPsiPartialPostSample = PsiPartialPostSample(BatchIndex,:);
            BatchLogPsiDensity = LogPsiDensity(BatchIndex,:);
         else;
            BatchPsiPostSample = [];
            BatchPsiPartialPostSample = [];
            BatchLogPsiDensity = [];
         end;
         %
         BatchNState = randn('state');
         BatchUState = rand('state');
         SetLogText(['   Saving: ' GetFilename(PostSampleFile)]);
         save(PostSampleFile,'BatchOmegaPostSample','BatchPiPostSample','BatchPsiPostSample','BatchOmegaPartialPostSample','BatchPsiPartialPostSample','BatchLogPiDensity','BatchLogPsiDensity','BatchNState','BatchUState');
         %
         NBatch = NBatch+1;
         NBatchDraws = 1;
         BatchIndex = BatchIndex+NumBatchDraws;
      end;
      drawnow;
   end;
   GibbsDrawTime = toc;
   %
   % close the current waitbar
   %
   if DSGEModel.ShowProgress==1;
      if ishandle(mywait)==1;
         set(mywait,'UserData','done');
         close(mywait);
         drawnow;
         pause(0.02);
      end;
   end;
else;
   %
   % diffuse prior on Pi
   %
   % Set up initial values.
   % NOTE: Lower case means that we use vectorized values
   %
   ThetaBarPi = PiInit;
   if q>0;
      thetaBarPsi = vec(PsiInit);
      ThetaBarPsi = PsiInit;
      ThetaPartialPsi = PsiPartialInit;
   else;
      thetaBarPsi = [];
      ThetaBarPsi = [];
      ThetaPartialPsi = [];
   end;
   OmegaBar = OmegaInit;
   OmegaPartial = OmegaPartialInit;
   invOmegaPi = 0;
   invOmegaPsi = inv(OmegaPsi);
   df = T+qDF;
   %
   % compute the value of the log likelihood + the log prior at the posterior mode
   %
   LogPosterior = BVARLogPosteriorDiffuse(OmegaMode,PiMode,PsiMode,Amatrix,qDF,thetaPsi,OmegaPsi,x,X,d,dLag);
   %
   % Calculate the log of the full conditional posterior of Omega at the mode
   %
   LogOmegaDensity = logInverseWishartPDF(OmegaMode,(df+p+1)*BVAROmegaMinnesota(PiMode,PsiMode,Amatrix,qDF,x,X,d,dLag),df);
   %
   % check if we should show a progress dialog
   %
   if DSGEModel.ShowProgress==1;
      ProgressStructure.title = ['Progress for ' num2str(NumPostDraws,'%0.0f') ' posterior draws'];
      ProgressStructure.facecolor = CurrINI.progress_facecolor;
      ProgressStructure.startfacecolor = CurrINI.progress_startfacecolor;
      ProgressStructure.edgecolor = CurrINI.progress_edgecolor;
      ProgressStructure.bgcolor = CurrINI.progress_bgcolor;
      ProgressStructure.stop = 0;
      ProgressStructure.clock = DSGEModel.ShowProgressClock;
      ProgressStructure.label = 'Number of draws:';
      %
      ProgressStructure.name = 'Gibbs Sampling of Bayesian VAR';
      ProgressStructure.CurrINI = CurrINI;
      mywait = ProgressDLG(0,ProgressStructure);
      set(mywait,'Color',get(0,'defaultuicontrolbackgroundcolor'));
      drawnow;
   end;
   %
   tic;
   for NPostDraws=StartPostDraws:NumPostDraws;
      %
      % update the progress dialog
      %
      if DSGEModel.ShowProgress==1;
         abort = get(mywait,'UserData');
         if (strcmp(abort,'cancel')==1)|(strcmp(abort,'done')==1);
            break;
         else;
            if NPostDraws<=NumPostDraws;
               ProgressDLG([NPostDraws/NumPostDraws NPostDraws]);
            end;
         end;
      end;
      %
      % we start by drawing OmegaBar conditional on Psi and Pi
      %
      OmegaBar = InvWishartRndFcn((df+p+1)*BVAROmegaMinnesota(ThetaBarPi,ThetaBarPsi,Amatrix,qDF,x,X,d,dLag),df);
      %
      % draw OmegaPartial with fixed value of Pi
      %
      OmegaPartial = InvWishartRndFcn((df+p+1)*BVAROmegaMinnesota(PiMode,ThetaPartialPsi,Amatrix,qDF,x,X,d,dLag),df);
      %
      % Next we draw vec(Pi) conditional on Psi and Omega
      %
      [MeanThetaBarPi,OmegaBarPi] = BVARPiMeanNormalCond(ThetaBarPsi,0,invOmegaPi,x,X,d,dLag);
      thetaBarPi = MultiNormalRndFcn(vec(MeanThetaBarPi),kron(OmegaBarPi,OmegaBar),1);
      ThetaBarPi = reshape(thetaBarPi,p,p*k);
      %
      % now we draw vec(Psi) conditional on Pi and Omega
      %
      if q>0;
         [MeanThetaBarPsi,SigmaBarPsi] = BVARPsiMean(OmegaBar,ThetaBarPi,thetaPsi,invOmegaPsi,x,X,d,dLag);
         thetaBarPsi = MultiNormalRndFcn(vec(MeanThetaBarPsi),SigmaBarPsi,1);
         ThetaBarPsi = reshape(thetaBarPsi,p,q);
         %
         % draw ThetaPartialPsi
         %
         [MeanThetaPartialPsi,SigmaPartialPsi] = BVARPsiMean(OmegaPartial,PiMode,thetaPsi,invOmegaPsi,x,X,d,dLag);
         thetaPartialPsi = MultiNormalRndFcn(vec(MeanThetaPartialPsi),SigmaPartialPsi,1);
         ThetaPartialPsi = reshape(thetaPartialPsi,p,q);
         %
         % compute the height of the log density of PsiMode given OmegaPartial and PiMode
         %
         [V,D] = eig(SigmaPartialPsi);
         LogPsiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PsiMode),vec(MeanThetaPartialPsi),sum(log(diag(D))),(V/D)*V');
%         LogPsiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PsiMode),vec(MeanThetaPartialPsi),log(det(SigmaPartialPsi)),SigmaPartialPsi\eye(p*q));
      end;
      %
      % compute the height of the log density of PiMode given ThetaBarPsi and OmegaBar
      %
      [MeanThetaBarPi,OmegaBarPi] = BVARPiMeanNormalCond(ThetaBarPsi,0,invOmegaPi,x,X,d,dLag);
      SigmaBarPi = kron(OmegaBarPi,OmegaBar);
      [V,D] = eig(SigmaBarPi);
      LogPiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PiMode),vec(MeanThetaBarPi),sum(log(diag(D))),(V/D)*V');
%      LogPiDensity(NPostDraws,1) = logMultiNormalPDF(vec(PiMode),vec(MeanThetaBarPi),log(det(SigmaBarPi)),SigmaBarPi\eye(p*p*k));
      %
      % Store draws
      %
      OmegaPostSample(NPostDraws,:) = vech(OmegaBar)';
      OmegaPartialPostSample(NPostDraws,:) = vech(OmegaPartial)';
      PiPostSample(NPostDraws,:) = thetaBarPi';
      if q>0;
         PsiPostSample(NPostDraws,:) = thetaBarPsi';
         PsiPartialPostSample(NPostDraws,:) = thetaPartialPsi';
      end;
      %
      % Check the batch counter
      %
      if NBatchDraws<NumBatchDraws;
         %
         % increase
         %
         NBatchDraws = NBatchDraws+1;
      else;
         %
         % save a batch
         %
         PostSampleFile = [DSGEModel.OutputDirectory '\gibbs\BVAR-Gibbs-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NBatch,'%0.0f') '.' num2str(NumBatches,'%0.0f') '.' num2str(NumPostDraws,'%0.0f') '.mat'];
         BatchOmegaPostSample = OmegaPostSample(BatchIndex,:);
         BatchOmegaPartialPostSample = OmegaPartialPostSample(BatchIndex,:);
         BatchPiPostSample = PiPostSample(BatchIndex,:);
         BatchLogPiDensity = LogPiDensity(BatchIndex,:);
         if q>0;
            BatchPsiPostSample = PsiPostSample(BatchIndex,:);
            BatchPsiPartialPostSample = PsiPartialPostSample(BatchIndex,:);
            BatchLogPsiDensity = LogPsiDensity(BatchIndex,:);
         else;
            BatchPsiPostSample = [];
            BatchPsiPartialPostSample = [];
            BatchLogPsiDensity = [];
         end;
         %
         BatchNState = randn('state');
         BatchUState = rand('state');
         SetLogText(['   Saving: ' GetFilename(PostSampleFile)]);
         save(PostSampleFile,'BatchOmegaPostSample','BatchPiPostSample','BatchPsiPostSample','BatchOmegaPartialPostSample','BatchPsiPartialPostSample','BatchLogPiDensity','BatchLogPsiDensity','BatchNState','BatchUState');
         %
         NBatch = NBatch+1;
         NBatchDraws = 1;
         BatchIndex = BatchIndex+NumBatchDraws;
      end;
      drawnow;
   end;
   GibbsDrawTime = toc;
   %
   % close the current waitbar
   %
   if DSGEModel.ShowProgress==1;
      if ishandle(mywait)==1;
         set(mywait,'UserData','done');
         close(mywait);
         drawnow;
         pause(0.02);
      end;
   end;
end;
%
% check if we cancelled or not
%
if strcmp(lower(abort),'cancel')==1;
   SetLogText(strvcat(' ','Gibbs sampling of Bayesian VAR was aborted.'));
   return;
end;
SetLogText(strvcat(' ','Gibbs sampling of Bayesian VAR completed!',[' * Time:                              ' StringTrim(datestr(now,14)) ', ' GetCurrentDate],[' * Posterior sampling time:           ' SecondsToHoursStr(GibbsDrawTime) ' (hrs:min:sec)']));
clear('PsiPartialPostSample','OmegaPartialPostSample');
%
% Calculate the marginal likelihood
%
SetLogText(strvcat(' ','Computing the marginal likelihood...',' * Chib''s marginal likelihood identity based estimator',[' * Time:                              ' StringTrim(datestr(now,14))]));
tic;
LogMargs = MargLikeChib(LogPiDensity,LogPsiDensity,LogPosterior,LogOmegaDensity,q,NumBurnin,DSGEModel.SequentialML,DSGEModel.SequentialStartIterationValue,DSGEModel.SequentialStepLengthValue);
ChibTime = toc;
if isempty(LogMargs)==1;
   SetLogText(strvcat('Estimation error...'));
   return;
end;
SetLogText(strvcat('Estimation completed!',[' * Estimation time:                   ' SecondsToHoursStr(ChibTime) ' (hrs:min:sec)']));
clear('LogPiDensity','LogPsiDensity');
%
% check if we should show a plot of the log marginal likelihood
%
GraphDirectory = [DSGEModel.OutputDirectory '\graphs'];
if DSGEModel.SequentialML==1;
   SetLogText(strvcat(' ','Creating plot of the sequential Chib estimator...'));
   FigHandle =figure('Visible','off', ...
          'Position',[(CurrINI.scrsz(3)-700)/2 (CurrINI.scrsz(4)-500)/2 700 500], ...
          'Tag','MargLikeChib', ...
          'Name','Sequential Marginal Likelihood - Chib');
   %
   SetFigureProperties(FigHandle,CurrINI);
   %
   Xvalues = LogMargs(:,1)+NumBurnin;
   Yvalues = LogMargs(:,2);
   %
   plot(Xvalues,Yvalues,'LineStyle','-','Color',[0.75 0 0],'LineWidth',2,'Marker','.');
   hold('off');
   SetAxesFonts(CurrINI,'Title','Chib estimator','XLabel',['Gibbs draw number (with ' num2str(NumBurnin,'%0.0f') ' burn in draws)'],'YLabel','Log marginal likelihood');
   axis('tight');
   %
   % check the YLim
   %
   ylim = get(gca,'YLim');
   add_on = 0.1*(ylim(2)-ylim(1));
   ylim = [ylim(1)-add_on ylim(2)+add_on];
   set(gca,'YLim',ylim);
   %
   % save the grasph
   %
   if DirExist(GraphDirectory);
      file = [GraphDirectory '\LogMarginalLikelihood-BVAR-Chib-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NumPostDraws,'%0.0f') '.fig'];
      SetLogText(strvcat('   Saving plot to:',['   ' GraphDirectory '...']));
      saveas(FigHandle,file);
   end;
   %
   % display the graph
   %
   SetToolbarImages(FigHandle,CurrINI.images);
   set(FigHandle,'Visible','on','CloseRequestFcn','delete(gcf); drawnow; pause(0.02);');
   drawnow;
   drawnow;
end;
%
% compute summary statistics about the parameters
%
SetLogText(strvcat(' ','Computing summary information about the parameters...',[' * Time:                              ' StringTrim(datestr(now,14))]));
tic;
if q>0;
   SetLogText(' * Psi - The steady state parameters...');
   PsiParData = ParameterSummary(PsiPostSample(NumBurnin+1:NumPostDraws,:));
   SetLogText('   DONE');
   clear('PsiPostSample');
else;
   PsiParData = [];
end;
SetLogText(' * Pi - The parameters on lagged endogenous...');
PiParData = ParameterSummary(PiPostSample(NumBurnin+1:NumPostDraws,:));
SetLogText(strvcat('   DONE',' * Pi - The covariance parameters...'));
clear('PiPostSample');
OmegaParData = ParameterSummary(OmegaPostSample(NumBurnin+1:NumPostDraws,:));
ParTime = toc;
SetLogText(strvcat('Computation completed!',[' * Estimation time:                   ' SecondsToHoursStr(ParTime) ' (hrs:min:sec)']));
clear('OmegaPostSample');
%
% Ready the print the results
%
SetLogText(strvcat(' ','Writing posterior sampling summary results to text file:',['   Directory: ' DSGEModel.OutputDirectory],['   File name: BVAR-GibbsSummary-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NumPostDraws,'%0.0f') '.txt']));
PrintBVARPosteriorSampling(DSGEModel,CurrINI,LogMargs,PsiParData,PiParData,OmegaParData,PsiMode,PiMode,OmegaMode,NumPostDraws,NumBurnin,OverallTightness,CrossEqTightness,HarmonicLagDecay,Amatrix,qDF,NewStartYear,NewStartPeriod,GibbsDrawTime);
%
% Display the results file
%
TextFile = [DSGEModel.OutputDirectory '\BVAR-GibbsSummary-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '-' num2str(NumPostDraws,'%0.0f') '.txt'];
TextGUI(TextFile,['Bayesian VAR Posterior Sampling - ' DSGEModel.NameOfModel],[(CurrINI.scrsz(3)-1000)/2 32 1000 CurrINI.scrsz(4)-100],'Results',50,CurrINI,1,1);

%
% end of BVARPosteriorSampling.m
%
