function BVARPosteriorModeEstimation(DSGEModel,CurrINI,controls)
% BVARPosteriorModeEstimation: Estimates the posterior mode of the joint posterior density
%                              of the parameters of the BVAR model.
%
% USAGE:
%
%       BVARPosteriorModeEstimation(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-2009 European Central Bank.
%
%                       First version: January 12, 2007.
%                        This version: March 9, 2009.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:
%
% * 15-01-2007: Completed the Minnesota-style posterior mode estimation part.
%
% * 17-01-2006: Added the normal conditional on the covariance matrix prior and the diffuse
%               prior posterior mode estimation part.
%
% * 18-01-2007: Corrected the calls to the Omega estimator.
%
% * 19-01-2007: Made sure that the log posterior value is also printed.
%
% * 25-01-2007: Updated the documentation.
%
% * 30-01-2007: Took the function name change from "BVARPsiMeanMinnesota" to "BVARPsiMean" into account. Added some
%               more information about the estimation to the log window.
%
% * 31-01-2007: The files with the posterior mode estimation results and the posterior mode estimates now have
%               new names based on the prior used.
%
% * 05-02-2007: Added 10 new variables that are saved to the output file.
%
% * 08-02-2007: Fixed OmegaInit so that it properly takes the assumed prior into account.
%
% * 09-02-2007: Added additional variables that are stored in the posterior mode stimation output file.
%               These outputs are used by the posterior sampling routine.
%
% * 13-11-2007: Updated the documentation.
%
% * 02-01-2008: Removed code using "DSGEModel.InitializeKalman".
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 23-05-2008: Updated the documentation.
%
% * 28-05-2008: Added a wait dialog when the posterior mode estimation routine is running without
%               displaying a progress dialog.
%
% * 19-12-2008: Updated the documentation.
%
% * 09-03-2009: Added quotation marks around paths for the dos copy command.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

SetLogText(strvcat(' ','Initializing Bayesian VAR posterior mode estimation...',' * Setting up the data matrices'));
%
% 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);
%
% Determine the first observation to use when computing log-likelihood with the Kalman filter
%
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;
%
% Setup all the basic data matrices
%
x = DSGEModel.Y(DSGEModel.BVARY,FirstPeriod+KalmanFirstObservation-1:LastPeriod);
T = size(x,2);
X = zeros(p*k,T);
for i=1:k;
   X(1+(p*(i-1)):p*i,:) = DSGEModel.Y(DSGEModel.BVARY,FirstPeriod+KalmanFirstObservation-1-i:LastPeriod-i);
end;
if isempty(DSGEModel.BVARX)==0;
   d = DSGEModel.X(DSGEModel.BVARX,FirstPeriod+KalmanFirstObservation-1:LastPeriod);
   dLag = zeros(q*k,T);
   for i=1:k;
      dLag(1+(q*(i-1)):q*i,:) = -DSGEModel.X(DSGEModel.BVARX,FirstPeriod+KalmanFirstObservation-1-i:LastPeriod-i);
   end;
else;
   d = [];
   dLag = [];
end;
%
% Compute the Maximum likelihood estimate of the residual covariance matrix of the VAR
%
Z = [d;X];
BML = (x*Z')/(Z*Z');
Epsilon = x-(BML*Z);
OmegaML = (1/T)*Epsilon*Epsilon';
PiInit = BML(:,q+1:size(BML,2));
if q>0;
   PsiInit = (eye(p)-(PiInit*SumMatrix(p,k)))*BML(:,1:q);
else;
   PsiInit = [];
end;
%
% actual sample that we will use
%
[NewStartYear,NewStartPeriod] = AdjustSampleStart(DSGEModel.SubBeginYear,DSGEModel.SubBeginPeriod,DSGEModel.DataFrequency,KalmanFirstObservation-1);
%
SetLogText(strvcat('   STATUS: OK',' * Collecting steady state prior data'));
%
% Collect all the data for the prior
%
[stat,msg] = dos(['copy /Y "' DSGEModel.SteadyStatePriorFile '" "' pwd '\tmp"']);
%
% collect the steady state prior data
%
eval(['PriorMatrix = ' GetName(DSGEModel.SteadyStatePriorFile) '(DSGEModel.BVARY,DSGEModel.BVARX);']);
%
thetaPsi = PriorMatrix(:,1);
OmegaPsi = diag(PriorMatrix(:,2).*PriorMatrix(:,2));
if det(OmegaPsi)<1.0e-10;
   SetLogText(strvcat('   STATUS: Determinant of steady state prior covariance','           matrix is too small',' ','Posterior mode estimation of Bayesian VAR was aborted.'));
   txt = ['The determinant of the prior covariance matrix is: ' num2str(det(OmegaPsi)) '. This is too small. YADA has therefore aborted posterior mode estimation of the Bayesian VAR.'];
   About(txt,'information','YADA - BVAR Model Problem',200,500,CurrINI);
   return;
end;
%
% Lagged parameters
%
SetLogText(strvcat('   STATUS: OK',' * Collecting prior data for parameters on lagged variables'));
%
% Collect the hyperparameters
%
PriorLevelMeanStr = get(controls.bvar.priorlevelmean,'String');
PriorLevelMean = str2num(StringTrim(PriorLevelMeanStr(DSGEModel.PriorLevelMeanValue,:)));
%
PriorDiffMeanStr = get(controls.bvar.priordiffmean,'String');
PriorDiffMean = str2num(StringTrim(PriorDiffMeanStr(DSGEModel.PriorDiffMeanValue,:)));
%
thetaPi = zeros(p*p*k,1);
DiffVar = 1-DSGEModel.levels;
DiffVar = DiffVar(DSGEModel.BVARY);
LevVar = DSGEModel.levels(DSGEModel.BVARY);
Pi1 = diag((PriorLevelMean*LevVar)+(PriorDiffMean*DiffVar));
thetaPi(1:p*p,1) = vec(Pi1);
%
OverallTightnessStr = get(controls.bvar.overalltightness,'String');
OverallTightness = str2num(StringTrim(OverallTightnessStr(DSGEModel.OverallTightnessValue,:)));
%
CrossEqTightnessStr = get(controls.bvar.crosseqtightness,'String');
CrossEqTightness = str2num(StringTrim(CrossEqTightnessStr(DSGEModel.CrossEqTightnessValue,:)));
%
HarmonicLagDecayStr = get(controls.bvar.harmoniclagdecay,'String');
HarmonicLagDecay = str2num(StringTrim(HarmonicLagDecayStr(DSGEModel.HarmonicLagDecayValue,:)));
%
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)]));
   OmegaPi = MinnesotaPrior(OverallTightness,CrossEqTightness,HarmonicLagDecay,diag(OmegaML),k);
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)]));
   OmegaPi = NormalConditionalPrior(OverallTightness,HarmonicLagDecay,p,k);
else;
   %
   % disffuse prior on Pi
   %
   thetaPi = [];
   OmegaPi = [];
   SetLogText('   PRIOR: Diffuse');
end;
SetLogText(strvcat('   STATUS: OK',' * Collecting prior data for covariance parameters'));
%
% Covariance matrix parameters
%
VarianceTightnessStr = get(controls.bvar.variancetightness,'String');
VarianceTightness = str2num(StringTrim(VarianceTightnessStr(DSGEModel.VarianceTightnessValue,:)));
%
if DSGEModel.OmegaPriorType==1;
   %
   % diffuse prior
   %
   Amatrix = 0;
   qDF = 0;
   SetLogText('   PRIOR: Diffuse');
else;
   %
   % inverted Wishart
   %
   SetLogText('   PRIOR: Inverse Wishart');
   qDF = p+get(controls.bvar.wishartdf,'Value')+1;
   if DSGEModel.WishartType==1;
      %
      % We use the ML estimate of Omega for the A matrix
      %
      Amatrix = OmegaML;
      SetLogText(strvcat('          - ML estimate of Omega for A',['          - Degrees of freedom:       ' num2str(qDF,'%0.0f')]));
   else;
      %
      % collect the variance tightness hyperparameter
      %
      Amatrix = VarianceTightness*eye(p);
      SetLogText(strvcat(['          - variance tightness:       ' num2str(VarianceTightness)],['          - Degrees of freedom:       ' num2str(qDF,'%0.0f')]));
   end;
end;
%
% Start the optimization to locate the posterior mode
%
MaxIterStr = get(controls.optimize.maximumiterations,'String');
MaxIterStr = StringTrim(strrep(MaxIterStr(DSGEModel.OptMaxIterationsValue,:),',',''));
MaxIter = str2num(MaxIterStr);
%
ToleranceStr = get(controls.optimize.tolerance,'String');
ToleranceStr = StringTrim(ToleranceStr(DSGEModel.OptToleranceValue,:));
ToleranceValue = str2num(ToleranceStr);
%
% We're OK. Start up the posterior mode estimation
%
SetLogText(strvcat('   STATUS: OK',' ','Starting posterior mode estimation.',[' * Time:                              ' StringTrim(datestr(now,14))],[' * Tolerance level:                   ' ToleranceStr],[' * Maximum number of iterations:      ' MaxIterStr],[' * 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 lags:                    ' num2str(k,'%0.0f')],[' * Number of exogenous variables:     ' num2str(q,'%0.0f')]));
%
abort = '';
it = 1;
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;
   else;
      thetaBarPsi = [];
      ThetaBarPsi = [];
   end;
   OmegaInit = (1/(T+qDF+p+1))*((Epsilon*Epsilon')+Amatrix);
   OmegaBar = OmegaInit;
   %
   % Evaluate the log posterior at these values
   %
   logPost = BVARLogPosteriorMinnesota(OmegaBar,ThetaBarPi,ThetaBarPsi,Amatrix,qDF,thetaPsi,OmegaPsi,thetaPi,OmegaPi,x,X,d,dLag);
   ConvCrit = 1;
   invOmegaPi = inv(OmegaPi);
   invOmegaPsi = inv(OmegaPsi);
   %
   % check if we should show a progress dialog
   %
   if DSGEModel.ShowProgress==1;
      ProgressStructure.title = ['Progress for ' MaxIterStr ' 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 = 'Convergence Value:';
      %
      ProgressStructure.name = 'Posterior Mode Estimation of Bayesian VAR';
      ProgressStructure.CurrINI = CurrINI;
      mywait = ProgressDLG(0,ProgressStructure);
      set(mywait,'Color',get(0,'defaultuicontrolbackgroundcolor'));
      drawnow;
   else;
      %
      % show a wait dialog
      %
      txt = ['Please wait while YADA attempts to estimate the posterior mode. A maximum of ' MaxIterStr ' iterations are used.'];
      WaitHandle = WaitDLG(txt,'information','Posterior Mode Estimation of Bayesian VAR',500,200,CurrINI,0);
      WaitControls = get(WaitHandle,'UserData');
   end;
   %
   tic;
   while (it<=MaxIter)&(ConvCrit>ToleranceValue);
      oldlogPost = logPost;
      %
      % check if we should update the progress dialog
      %
      if DSGEModel.ShowProgress==1;
         abort = get(mywait,'UserData');
         if (strcmp(abort,'cancel')==1)|(strcmp(abort,'done')==1);
            break;
         else;
            ProgressDLG([it/MaxIter ConvCrit]);
         end;
      end;
      %
      % compute new values for Omega, Pi, and Psi
      %
      % Step 1: New (conditional mean) value for Pi
      %
      ThetaBarPi = BVARPiMeanMinnesota(OmegaBar,ThetaBarPsi,thetaPi,invOmegaPi,x,X,d,dLag);
      %
      % Step 2: New (conditional mean) for Psi.
      %
      if q>0;
         ThetaBarPsi = BVARPsiMean(OmegaBar,ThetaBarPi,thetaPsi,invOmegaPsi,x,X,d,dLag);
      end;
      %
      % Step 3: New value for OmegaBar
      %
      OmegaBar = BVAROmegaMinnesota(ThetaBarPi,ThetaBarPsi,Amatrix,qDF,x,X,d,dLag);
      %
      % Evaluate the log posterior at these values
      %
      logPost = BVARLogPosteriorMinnesota(OmegaBar,ThetaBarPi,ThetaBarPsi,Amatrix,qDF,thetaPsi,OmegaPsi,thetaPi,OmegaPi,x,X,d,dLag);
      ConvCrit = abs(logPost-oldlogPost);
      it = it+1;
   end;
   MaxTime = toc;
   %
   % done with the optimization
   %
   NumIter = it-1;
   %
   % 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 prior 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;
   else;
      thetaBarPsi = [];
      ThetaBarPsi = [];
   end;
   invOmegaPi = inv(OmegaPi);
   invOmegaPsi = inv(OmegaPsi);
   OmegaInit = (1/(T+qDF+(p*k)+p+1))*((Epsilon*Epsilon')+Amatrix+((PiInit-ThetaPi)*invOmegaPi*(PiInit-ThetaPi)'));
   OmegaBar = OmegaInit;
   %
   % Evaluate the log posterior at these values
   %
   logPost = BVARLogPosteriorNormalCond(OmegaBar,ThetaBarPi,ThetaBarPsi,Amatrix,qDF,thetaPsi,OmegaPsi,ThetaPi,OmegaPi,x,X,d,dLag);
   ConvCrit = 1;
   %
   % check if we should show a progress dialog
   %
   if DSGEModel.ShowProgress==1;
      ProgressStructure.title = ['Progress for ' MaxIterStr ' 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 = 'Convergence Value:';
      %
      ProgressStructure.name = 'Posterior Mode Estimation of Bayesian VAR';
      ProgressStructure.CurrINI = CurrINI;
      mywait = ProgressDLG(0,ProgressStructure);
      set(mywait,'Color',get(0,'defaultuicontrolbackgroundcolor'));
      drawnow;
   else;
      %
      % show a wait dialog
      %
      txt = ['Please wait while YADA attempts to estimate the posterior mode. A maximum of ' MaxIterStr ' iterations are used.'];
      WaitHandle = WaitDLG(txt,'information','Posterior Mode Estimation of Bayesian VAR',500,200,CurrINI,0);
      WaitControls = get(WaitHandle,'UserData');
   end;
   %
   tic;
   while (it<=MaxIter)&(ConvCrit>ToleranceValue);
      oldlogPost = logPost;
      %
      % check if we should update the progress dialog
      %
      if DSGEModel.ShowProgress==1;
         abort = get(mywait,'UserData');
         if (strcmp(abort,'cancel')==1)|(strcmp(abort,'done')==1);
            break;
         else;
            ProgressDLG([it/MaxIter ConvCrit]);
         end;
      end;
      %
      % compute new values for Omega, Pi, and Psi
      %
      % Step 1: New (conditional mean) value for Pi
      %
      ThetaBarPi = BVARPiMeanNormalCond(ThetaBarPsi,ThetaPi,invOmegaPi,x,X,d,dLag);
      %
      % Step 2: New (conditional mean) for Psi.
      %
      if q>0;
         ThetaBarPsi = BVARPsiMean(OmegaBar,ThetaBarPi,thetaPsi,invOmegaPsi,x,X,d,dLag);
      end;
      %
      % Step 3: New value for OmegaBar
      %
      OmegaBar = BVAROmegaNormal(ThetaBarPi,ThetaBarPsi,Amatrix,qDF,ThetaPi,invOmegaPi,x,X,d,dLag);
      %
      % Evaluate the log posterior at these values
      %
      logPost = BVARLogPosteriorNormalCond(OmegaBar,ThetaBarPi,ThetaBarPsi,Amatrix,qDF,thetaPsi,OmegaPsi,ThetaPi,OmegaPi,x,X,d,dLag);
      ConvCrit = abs(logPost-oldlogPost);
      it = it+1;
   end;
   MaxTime = toc;
   %
   % done with the optimization
   %
   NumIter = it-1;
   %
   % 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;
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;
   else;
      thetaBarPsi = [];
      ThetaBarPsi = [];
   end;
   OmegaInit = (1/(T+qDF+p+1))*((Epsilon*Epsilon')+Amatrix);
   OmegaBar = OmegaInit;
   %
   % Evaluate the log posterior at these values
   %
   logPost = BVARLogPosteriorDiffuse(OmegaBar,ThetaBarPi,ThetaBarPsi,Amatrix,qDF,thetaPsi,OmegaPsi,x,X,d,dLag);
   ConvCrit = 1;
   invOmegaPi = 0;
   invOmegaPsi = inv(OmegaPsi);
   %
   % check if we should show a progress dialog
   %
   if DSGEModel.ShowProgress==1;
      ProgressStructure.title = ['Progress for ' MaxIterStr ' 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 = 'Convergence Value:';
      %
      ProgressStructure.name = 'Posterior Mode Estimation of Bayesian VAR';
      ProgressStructure.CurrINI = CurrINI;
      mywait = ProgressDLG(0,ProgressStructure);
      set(mywait,'Color',get(0,'defaultuicontrolbackgroundcolor'));
      drawnow;
   else;
      %
      % show a wait dialog
      %
      txt = ['Please wait while YADA attempts to estimate the posterior mode. A maximum of ' MaxIterStr ' iterations are used.'];
      WaitHandle = WaitDLG(txt,'information','Posterior Mode Estimation of Bayesian VAR',500,200,CurrINI,0);
      WaitControls = get(WaitHandle,'UserData');
   end;
   %
   tic;
   while (it<=MaxIter)&(ConvCrit>ToleranceValue);
      oldlogPost = logPost;
      %
      % check if we should update the progress dialog
      %
      if DSGEModel.ShowProgress==1;
         abort = get(mywait,'UserData');
         if (strcmp(abort,'cancel')==1)|(strcmp(abort,'done')==1);
            break;
         else;
            ProgressDLG([it/MaxIter ConvCrit]);
         end;
      end;
      %
      % compute new values for Omega, Pi, and Psi
      %
      % Step 1: New (conditional mean) value for Pi
      %
      ThetaBarPi = BVARPiMeanNormalCond(ThetaBarPsi,0,invOmegaPi,x,X,d,dLag);
      %
      % Step 2: New (conditional mean) for Psi.
      %
      if q>0;
         ThetaBarPsi = BVARPsiMean(OmegaBar,ThetaBarPi,thetaPsi,invOmegaPsi,x,X,d,dLag);
      end;
      %
      % Step 3: New value for OmegaBar
      %
      OmegaBar = BVAROmegaMinnesota(ThetaBarPi,ThetaBarPsi,Amatrix,qDF,x,X,d,dLag);
      %
      % Evaluate the log posterior at these values
      %
      logPost = BVARLogPosteriorDiffuse(OmegaBar,ThetaBarPi,ThetaBarPsi,Amatrix,qDF,thetaPsi,OmegaPsi,x,X,d,dLag);
      ConvCrit = abs(logPost-oldlogPost);
      it = it+1;
   end;
   MaxTime = toc;
   %
   % done with the optimization
   %
   NumIter = it-1;
   %
   % 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;
end;
if strcmp(lower(abort),'cancel')==1;
   SetLogText(strvcat(' ','Posterior mode estimation of Bayesian VAR was aborted.'));
   return;
end;
if it>MaxIter;
   SetLogText(strvcat(' ','Posterior mode estimation of Bayesian VAR didn''t converge.'));
   return;
end;
%
% the posterior mode estimation routine for the Bayesian VAR converged
%
SetLogText(strvcat(' ','Posterior mode estimation completed.',[' * Time:                              ' StringTrim(datestr(now,14))],[' * Optimization time:                 ' SecondsToHoursStr(MaxTime) ' (hrs:min:sec)'],[' * Number of iterations:              ' num2str(NumIter)],[' * Convergence value:                 ' num2str(ConvCrit)],[' * Value of the log posterior:        ' num2str(logPost)]));
%
% save the results to file
%
if q>0;
   Psi = ThetaBarPsi;
else;
   Psi = [];
end;
Pi = ThetaBarPi;
Omega = OmegaBar;
VariableNamesBVAR = DSGEModel.VariableNames(DSGEModel.BVARY,:);
if isempty(DSGEModel.BVARX)==0;
   XVariableNamesBVAR = DSGEModel.XVariableNames(DSGEModel.BVARX,:);
else;
   XVariableNamesBVAR = '';
end;
%
ModeFileDir = [DSGEModel.OutputDirectory '\mode'];
stat = MakeDir(DSGEModel.OutputDirectory,'mode');
if stat==1;
   SetLogText(strvcat(' ','Saving posterior mode data to disk using directory:',['   "' ModeFileDir '"...']));
   ModeFile = [ModeFileDir '\PosteriorMode-BVAR-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '.mat'];
   save(ModeFile,'Psi','Pi','Omega','VariableNamesBVAR','XVariableNamesBVAR','Amatrix','qDF','thetaPsi','OmegaPsi','thetaPi','OmegaPi','x','X','d','dLag','DSGEModel','OverallTightness','CrossEqTightness','HarmonicLagDecay','VarianceTightness');
end;
%
% Write output to file
%
PrintBVARPosteriorMode(DSGEModel,CurrINI,Psi,Pi,Omega,VariableNamesBVAR,XVariableNamesBVAR,thetaPsi,OmegaPsi,thetaPi,OmegaPi,OverallTightness,CrossEqTightness,HarmonicLagDecay,Amatrix,qDF,NewStartYear,NewStartPeriod,MaxTime,ToleranceValue,MaxIter,NumIter,ConvCrit,logPost);
%
% Display the results
%
ModeTextFile = [DSGEModel.OutputDirectory '\PosteriorModeBVAR-' DSGEModel.NameOfModel '-' num2str(DSGEModel.PriorType,'%0.0f') '.' num2str(DSGEModel.OmegaPriorType,'%0.0f') '.txt'];
TextGUI(ModeTextFile,['Bayesian VAR Posterior Mode Estimation - ' DSGEModel.NameOfModel],[(CurrINI.scrsz(3)-1000)/2 32 1000 CurrINI.scrsz(4)-100],'Results',50,CurrINI,1,1);

%
% end of BVARPosteriorModeEstimation.m
%
