function DSGEPosteriorSampling(DSGEModel,CurrINI,controls)
% DSGEPosteriorSampling: Runs the random walk Metropolis algorithm to sample from the
%                        posterior distribution of the parameters of the DSGE model.
%
% USAGE:
%
%       DSGEPosteriorSampling(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: November 27, 2006.
%                        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:
%
% * 28-11-2006: Continued work on the function.
%
% * 29-11-2006: Finished the random walk Metropolis posterior draws part. Added the modified harmonic mean
%               estimator of the marginal likelihood. Started work on the Chib & Jeliazkov estimator.
%
% * 30-11-2006: Finished work on the Chib & jeliazkov estimator.
%
% * 01-12-2006: Removed some unnecessary stuff.
%
% * 07-12-2006: Added the call to the print function and wrote the graph code.
%
% * 08-12-2006: Finished the call to the print function.
%
% * 11-12-2006: Added the possibility of loading data run using fewer draws.
%
% * 13-12-2006: Fixed the y-axis properties when plotting the log marginal likelihoods.
%
% * 18-12-2006: Added code for taking the value of DSGEModel.RandomNumberValue into account.
%
% * 21-12-2006: Made sure that the variable "PostSampleTime" is initialized.
%
% * 08-01-2007: Added a call to the clear function.
%
% * 27-02-2007: Added the call "SetAxesFonts(CurrINI,legHandle)" after the legend call.
%
% * 28-02-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.
%
% * 29-03-2007: Updated the documentation.
%
% * 10-04-2007: Added two new inverse Hessian estimators.
%
% * 11-04-2007: The function can now calculate results for multiple chains.
%
% * 12-04-2007: Fixed a bug with phiInit initialization when draws are recovered from a run with
%               fewer total posterior draws. The code had also been adjusted such that loading of
%               posterior draws for a lower number of parallel chains is possible. This case requires
%               that the total number of draws is equal across the number of parallel chains that are
%               checked.
%
% * 17-04-2007: The function now takes the selected maximum correlation for the inverse Hessian estimator
%               into account.
%
% * 18-04-2007: Added the average line to the sequential estimates of the log marginal likelihood using the
%               modified harmonic mean estimator. Also, the line widths have been made smaller and the
%               markers have been removed.
%
% * 13-11-2007: Updated the documentation.
%
% * 07-12-2007: Fixed a bug with initialization of the number of chains. When "ExitLoop" is 1, then the
%               "NChains" loop should also be stopped.
%
% * 11-12-2007: Added the variable "NumParam" instead of calling the length function for the number of
%               parameters. This has had a significant effect (about 5 percent!) on the computational
%               speed in simple models.
%
% * 21-12-2007: Specified unit-root state processes can now be handled.
%
% * 02-01-2008: Removed code using "DSGEModel.InitializeKalman".
%
% * 07-01-2007: Added drawnow calls after Query and About calls.
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 18-02-2008: The function now takes the option "ModifiedHessian" into account. When this variable
%               is 1 and the standard deviations from the modified Hessian are used, these standard
%               deviation are scaled up using the correlation structure to ensure that they correspond
%               to "marginal" rather than "conditional" standard deviations.
%
% * 05-03-2008: Made sure that matrices that are created when the variable "DSGEModel.RandomNumberValue"
%               is equal to unity are properly cleared. This variable is unity when the option "Use fixed
%               state to initialize random number generators" on the Miscellaneous tab is check marked.
%
% * 10-03-2008: The acceptance ratio is now also written to the estimation log window when saving a batch
%               of posterior draws.
%
% * 14-03-2008: Made sure that large matrices are cleared from memory once they've served their purpose.
%
% * 25-03-2008: Fixed a bug where "NChains" (full number of chains) instead of "NChain" (current chain)
%               was claimed to have been completed. Old chains can now be reused also when "OverwriteDraws"
%               is zero, given that a lower number of total chains is considered for copying. Fixed a bug
%               to ensure that phiInit is initialized when "OverwriteDraws" is zero.
%
% * 04-04-2008: The function now estimates the posterior mode of the marginal distributions via kernel
%               density estimation.
%
% * 13-05-2008: The check of running "DSGEModel.InitializeParameterFile" is no longer conducted. The values
%               from the posterior mode output are used instead.
%
% * 16-05-2008: Took the "DSGEModel.AIMTolerance" input for AiMSolver into account.
%
% * 19-06-2008: Made sure that the "AiMErrorData.mat" file is cleaned from the tmp directory after estimation
%               has finished.
%
% * 23-05-2008: Updated the documentation.
%
% * 27-05-2008: The mcode value 7 from AiM is now taken on board.
%
% * 28-05-2008: Added a wait dialog when the posterior sampler is run and no progress dialog is displayed.
%
% * 29-05-2008: Added the calls to GetCurrentDate.
%
% * 03-07-2008: Made sure that the measurement equation receives only the active state variable names.
%
% * 23-10-2008: The method for calculating the total number of parameters have been changed.
%
% * 29-10-2008: Allowed for loading of posterior draws from fewer draws with 1 chain also when the user has selected not
%               to overwrite old draws.
%
% * 21-01-2009: Updated the documentation.
%
% * 09-02-2009: Changed the input "length(DSGEModel.StateVariablePositions)" to "size(F,1)" when calling the
%               measurement equation function. This required that the call to the measurement equation function
%               was moved to after trying the solve the DSGE model.
%
% * 13-02-2009: Added code to allow for a user defined covariance matrix for the proposal density. The mcode
%               values 7 and 8 are now properly mentioned when encountered.
%
% * 19-02-2009: Fixed a bug when importing data from fewer number of posterior draws and using one batch
%               only. In this case, the "NBatchDraws" variable was not set correctly resulting in not
%               saving the draws to mat-file when they had been completed. Bug discovered by Peter Welz.
%
% * 26-02-2009: Added code for preparing output of the parameters that are defined via the update function,
%               i.e., the parameters that are a function of the estimated parameters.
%
% * 02-03-2009: Made sure that the call to the function "GraphPosteriorDensity" includes the 3 empty matrices
%               for the prior densities of the update parameters.
%
% * 09-03-2009: Added quotation marks around paths for the dos del and copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if FileExist([pwd '\tmp\AiMErrorData.mat'])==1;
   dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
end;
SetLogText(strvcat(' ','Initializing posterior sampling using random walk Metropolis...',[' * Attempting to load posterior mode data file...']));
%
% First we load the posterior mode estimates
%
ModeFile = [DSGEModel.OutputDirectory '\mode\PosteriorMode-' DSGEModel.NameOfModel '.mat'];
ErrorStr = '';
try;
   ModeData = load(ModeFile);
catch;
   ErrorStr = ['Unable to load the file "' ModeFile '". Message caught is: ' lasterr];
end;
if isempty(ErrorStr)==1;
   SetLogText(strvcat('   STATUS: OK',' * Attempting to load DSGE model data file...'));
else;
   SetLogText(strvcat(' ',['Unable to load the file "' GetFilename(ModeFile) '".']))
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   drawnow;
   return;
end;
%
phiMode = ModeData.phiMode;
NumParam = length(phiMode);
if DSGEModel.InverseHessianEstimator==1;
   SigmaMode = ModeData.InverseHessian;
   InvHessianTxt = 'Optimization routine';
elseif DSGEModel.InverseHessianEstimator==2;
   if sum(isnan(ModeData.phiFittedStd))==0;
      SigmaStd = sqrt(abs(diag(ModeData.InverseHessian)));
      SigmaCorr = ModeData.InverseHessian ./ (SigmaStd * SigmaStd');
      phiFittedStd = ModeData.phiFittedStd;
      oldphiFittedStd = phiFittedStd;
      if DSGEModel.ModifiedHessian==1;
         FullSet = (1:length(phiFittedStd));
         for i=1:length(phiFittedStd);
            iExcluded = setdiff(FullSet,i);
            phiFittedStd(i) = phiFittedStd(i)/sqrt(1-(SigmaCorr(i,iExcluded)*inv(SigmaCorr(iExcluded,iExcluded))*SigmaCorr(iExcluded,i)));
         end;
      end;
      SigmaMode = SigmaCorr .* (phiFittedStd * phiFittedStd');
   else;
      SetLogText(strvcat(' ','Quadratic fit to log posterior unavailable.'));
      txt = 'YADA has not fitted a quadratic to the evaluated log posterior. Please run the "Check posterior mode" function on the View menu first. Alternatively, you can use the estimator of the inverse Hessian from the optimization routine.';
      About(txt,'information','YADA - Modified Inverse Hessian',120,500,CurrINI);
      drawnow;
      return;
   end;
   InvHessianTxt = 'Quadratic modification';
elseif DSGEModel.InverseHessianEstimator==3;
   if sum(sum(isnan(ModeData.phiInverseHessianFinDiff)))==0;
      SigmaMode = ModeData.phiInverseHessianFinDiff;
   else;
      txt = 'YADA has not estimated the inverse Hessian using finite difference during the posterior mode estimation. Would you like to begin the finite difference estimation now?';
      answer = Query(txt,'information',140,'YADA - Finite Difference Inverse Hessian',500,'no',CurrINI);
      drawnow;
      if strcmp(lower(answer),'yes')==1;
         SigmaMode = [];
      else;
         SetLogText(strvcat(' ','Finite difference inverse Hessian unavailable.'));
         return;
      end;
   end;
   InvHessianTxt = 'Finite difference';
elseif DSGEModel.InverseHessianEstimator==4;
   %
   % try to load the file
   %
   CovMatData = load(DSGEModel.ParameterCovMatrix);
   SigmaMode = CovMatData.ParameterCovarianceMatrix;
   if size(SigmaMode,1)==NumParam;
      if IsPositiveDefinite(SigmaMode)==1;
         InvHessianTxt = 'Own matrix';
      else;
         SetLogText(strvcat(' ','The parameter covariance matrix is not positive definite.'));
         txt = ['The parameter covariance matrix you have selected from the file "' GetName(DSGEModel.ParameterCovMatrix) '" isn''t positive definite. YADA is therefore unable to continue with posterior sampling.'];
         About(txt,'information','YADA - Own Parameter Covariance Matrix',200,500,CurrINI);
         return;
      end;
   else;
      SetLogText(strvcat(' ',['The dimension of your parameter covariance matrix (' int2str(size(SigmaMode,1)) ') is not'],['compatible with the number of parameters to estimate (' int2str(NumParam) ').']));
      txt = ['The dimension of the parameter covariance matrix you have selected from the file "' GetName(DSGEModel.ParameterCovMatrix) '" is ' int2str(size(SigmaMode,1)) ', while the number of parameters to estimate is ' int2str(NumParam) '. YADA is therefore unable to continue with the posterior sampling.'];
      About(txt,'information','YADA - Own Parameter Covariance Matrix',200,500,CurrINI);
      return;
   end;
end;
%
% check if we should manipulate the correlations for the proposal density
%
if DSGEModel.MaxCorrelationValue>1;
   MaxCorrStr = get(controls.posterior.maxcorrelation,'String');
   MaxCorr = str2num(StringTrim(MaxCorrStr(DSGEModel.MaxCorrelationValue,:)));
   %
   % create the correlation matrix
   %
   StdVec = sqrt(diag(SigmaMode));
   CorrMat = SigmaMode ./ (StdVec*StdVec');
   for i=1:length(StdVec);
      CorrMat(i,i) = 0;
   end;
   ScaleFactor = max(max(abs(CorrMat)));
   ScaleMatrix = min(1,MaxCorr/ScaleFactor)*ones(length(StdVec),length(StdVec));
   for i=1:length(StdVec);
      ScaleMatrix(i,i) = 1;
   end;
   SigmaMode = SigmaMode .* ScaleMatrix;
else;
   MaxCorr = 1;
end;
%
% load additional information about the model.
%
thetaIndex = ModeData.thetaIndex;
UniformBounds = ModeData.UniformBounds;
LowerBound = ModeData.LowerBound;
thetaPositions = ModeData.thetaPositions;
thetaDist = ModeData.thetaDist;
PriorDist = ModeData.PriorDist;
ModelParameters = ModeData.ModelParameters;
ParameterNames = ModeData.ParameterNames;
thetaMode = ModeData.thetaMode;
%
% sampling data
%
DSGEModel.FirstPeriod = ModeData.FirstPeriod;
DSGEModel.LastPeriod = ModeData.LastPeriod;
DSGEModel.KalmanFirstObservation = ModeData.DSGEModel.KalmanFirstObservation;
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;
[NewStartYear,NewStartPeriod] = AdjustSampleStart(DSGEModel.SubBeginYear,DSGEModel.SubBeginPeriod,DSGEModel.DataFrequency,DSGEModel.KalmanFirstObservation-1);
%
% state vector initialization data
%
DSGEModel.UseOwnInitialState = ModeData.DSGEModel.UseOwnInitialState;
DSGEModel.InitialStateValues = ModeData.DSGEModel.InitialStateValues;
%
% Coverage probabilities for the Modified Harmonic Mean estimator of
% the marginal likelihood
%
CovStartStr = get(controls.dsge.coveragestart,'String');
DSGEModel.CovStartValue = str2num(StringTrim(CovStartStr(DSGEModel.CoverageStart,:)));
CovIncStr = get(controls.dsge.coverageincrement,'String');
DSGEModel.CovIncValue = str2num(StringTrim(CovIncStr(DSGEModel.CoverageIncrement,:)));
CovEndStr = get(controls.dsge.coverageend,'String');
DSGEModel.CovEndValue = str2num(StringTrim(CovEndStr(DSGEModel.CoverageEnd,:)));
%
% 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,:)),',',''));
%
% Now we need to copy files to the tmp directory and load data from DSGEModel.AIMDataFile
%
[stat,msg] = dos(['copy /Y "' GetPath(DSGEModel.AIMDataFile) 'compute_aim_matrices.m" "' pwd '\tmp"']);
[stat,msg] = dos(['copy /Y "' DSGEModel.MeasurementEquationFile '" "' pwd '\tmp"']);
if FileExist(DSGEModel.UpdateParameterFile)==1;
   [stat,msg] = dos(['copy /Y "' DSGEModel.UpdateParameterFile '" "' pwd '\tmp"']);
end;
%
ErrorFile = '';
try;
   AIMData = load(DSGEModel.AIMDataFile);
catch;
   ErrorStr = ['Unable to load the file "' DSGEModel.AIMDataFile '". Message caught is: ' lasterr];
   ErrorFile = DSGEModel.AIMDataFile;
end;
%
% Test if there are any problems
%
firstUpdatePos = 0;
lastUpdatePos = 0;
firstInitializePos = 0;
lastInitializePos = 0;
if isempty(ErrorStr)==1;
   if DSGEModel.UseOwnInitialState==1;
      if length(DSGEModel.InitialStateValues)==size(AIMData.endog,1);
         DSGEModel.InitialStateVector = DSGEModel.InitialStateValues(DSGEModel.StateVariablePositions);
      else;
         DSGEModel.InitialStateVector = zeros(length(DSGEModel.StateVariablePositions),1);
      end;
   else;
      DSGEModel.InitialStateVector = zeros(length(DSGEModel.StateVariablePositions),1);
   end;
   %
   if FileExist(DSGEModel.UpdateParameterFile)==1;
      SetLogText(strvcat('   STATUS: OK',[' * Checking ' GetName(DSGEModel.UpdateParameterFile) '...']));
      try;
         firstUpdatePos = length(fieldnames(ModelParameters))+1;
         eval(['ModelParameters = ' GetName(DSGEModel.UpdateParameterFile) '(ModelParameters);']);
         lastUpdatePos = length(fieldnames(ModelParameters));
      catch;
         ErrorStr = ['Error found in "' GetName(DSGEModel.UpdateParameterFile) '.m". Message caught is: ' lasterr];
         ErrorFile = DSGEModel.UpdateParameterFile;
      end;
   end;
end;
%
% Check if we should exit
%
if isempty(ErrorStr)==0;
   SetLogText(strvcat(['   STATUS: Bad input file ' GetName(ErrorFile)],' '));
   About(ErrorStr,'error',['Error - ' GetFilename(ErrorFile)],120,500,CurrINI);
   drawnow;
   return;
end;
%
% Try to solve the DSGE model at the posterior mode
%
SetLogText(strvcat('   STATUS: OK',' * Checking if the DSGE model has a unique convergent','   solution for posterior mode parameter values...'));
[mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
if mcode~=1;
   %
   % There doesn't exist a unique convergent solution for the initial parameter values.
   %
   if mcode==2;
      mcodeStr = 'Roots are not correctly computed by real_schur.'
   elseif mcode==3;
      mcodeStr = 'Too many big roots.';
   elseif mcode==35;
      mcodeStr = 'Too many big roots, and q(:,right) is singular.';
   elseif mcode==4;
      mcodeStr = 'Too few big roots.';
   elseif mcode==45;
      mcodeStr = 'Too few big roots, and q(:,right) is singular.';
   elseif mcode==5;
      mcodeStr = 'q(:,right) is singular.';
   elseif mcode==61;
      mcodeStr = 'Too many exact shiftrights.';
   elseif mcode==62;
      mcodeStr = 'Too many numeric shiftrights.';
   elseif mcode==7;
      mcodeStr = 'Infinite or NaN values detected.';
   elseif mcode==8;
      mcodeStr = 'The function "compute_aim_matrices" returns complex numbers.';
   else;
      mcodeStr = 'Return code not properly specified.';
   end;
   SetLogText(strvcat(['   STATUS: ' mcodeStr],' '));
   txt = ['The AiM solver provided the return code: ' num2str(mcode) ', i.e., "' mcodeStr '" The posterior mode values are no longer valid.'];
   About(txt,'information','YADA - DSGE Model Problem',120,500,CurrINI);
   drawnow;
   return;
end;
%
% do it again to measure the time
%
tic;
[mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
SolveTime = toc;
%
% check the state equations
%
[F,B0] = AiMtoStateSpace(BMatrix,SMatrix,DSGEModel.StateVariablePositions,DSGEModel.StateShockPositions);
if isempty(DSGEModel.UnitRootStates)==1;
   F1 = F;
else;
   KeepPos = GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates);
   F1 = F(KeepPos,KeepPos);
end;
if (max(abs(eig(F1)))>=1)|(IsAlmostEqual(max(abs(eig(F1))),1)==1);
   SetLogText(strvcat('   STATUS: Non-stationary state equation',' '));
   txt = ['The largest eigenvalue of the state equation is ' num2str(max(abs(eig(F1)))) '. The state equation is therefore non-stationary. The posterior mode values are no longer valid.'];
   About(txt,'information','YADA - DSGE Model Problem',120,500,CurrINI);
   drawnow;
   return;
end;
%
% check the measurement equation file for internal validity
%
if isempty(ErrorStr)==1;
   SetLogText(strvcat('   STATUS: OK',['   Maximum absolute error: ' num2str(MaxAbsError)],['   Solution time:          ' num2str(SolveTime) ' seconds'],[' * Checking ' GetName(DSGEModel.MeasurementEquationFile) '...']));
   try;
      eval(['[A,H,R] = ' GetName(DSGEModel.MeasurementEquationFile) '(ModelParameters,AIMData.endog(DSGEModel.StateVariablePositions,:),DSGEModel.VariableNames,DSGEModel.XVariableNames,DSGEModel.n,size(F,1),DSGEModel.k);']);
   catch;
      ErrorStr = ['Error found in "' GetName(DSGEModel.MeasurementEquationFile) '.m". Message caught is: ' lasterr];
      ErrorFile = DSGEModel.MeasurementEquationFile;
   end;
end;
%
% Check if we should exit
%
if isempty(ErrorStr)==0;
   SetLogText(strvcat(['   STATUS: Bad input file ' GetName(ErrorFile)],' '));
   About(ErrorStr,'error',['Error - ' GetFilename(ErrorFile)],120,500,CurrINI);
   drawnow;
   return;
end;
%
% Try to calculate the the log posterior at the posterior mode
%
SetLogText(strvcat('   STATUS: OK',' * Checking if the log posterior can be calculated...'));
%
ErrorStr = '';
try;
   logPostMode = logPosteriorPhiDSGE(phiMode,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
catch;
   ErrorStr = ['Unable to calculate the log posterior. Message caught is: ' lasterr];
end;
if FileExist([pwd '\tmp\AiMErrorData.mat'])==1;
   dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
end;
%
if ((exist('logPostMode')==0)|(logPostMode==1000000))&(isempty(ErrorStr)==1);
   SetLogText(strvcat('   STATUS: Unable to compute log posterior at posterior mode',' '));
   txt = 'Unable to compute the log posterior at the posterior mode values. You probably need to re-estimate the posterior mode.';
   About(txt,'information','YADA - DSGE Model Problem',120,500,CurrINI);
   drawnow;
   return;
end;
%
% there doesn't seem to be any problems with the model
%
tic;
logPostMode = -logPosteriorPhiDSGE(phiMode,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
LogPostTime = toc;
SetLogText(strvcat('   STATUS: OK',['   Computation time:       ' num2str(LogPostTime) ' seconds']));
%
% check if we should calculate a finite difference Hessian
%
StepLengthStr = get(controls.optimize.steplength,'String');
StepLength = str2num(StringTrim(StepLengthStr(DSGEModel.StepLengthHessian,:)));
if isempty(SigmaMode)==1;
   SetLogText(' * Estimating inverse Hessian with finite differences...');
   tic;
   VechHessian = YADAHessianMatrix('logPosteriorPhiDSGE',phiMode,StepLength,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
   FinDiffTime = toc;
   if sum(isnan(VechHessian))==0;
      SetLogText(strvcat('   STATUS: OK',['   Computation time:       ' num2str(FinDiffTime) ' seconds']));
      HessianMat = GeneralizedCholeski(reshape(VechHessian,NumParam,NumParam));
      phiInverseHessianFinDiff = inv(HessianMat);
      %
      % save to post mode file
      %
      save(ModeFile,'phiInverseHessianFinDiff','-APPEND');
      SigmaMode = phiInverseHessianFinDiff;
   else;
      %
      % bad Hessian encountered
      %
      SetLogText(strvcat('   STATUS: Bad estimator of inverse Hessian using','           finite differences',' '));
      txt = ['The inverse Hessian estimator using finite differences with a step length of ' num2str(StepLength) ' has at least one NaN entry. Please change the step length or select a different estimator of the inverse Hessian.'];
      About(txt,'information','YADA - Inverse Hessian Estimation Problem',150,500,CurrINI);
      drawnow;
      return;
   end;
end;
%
% Get the base values number of draws, 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 parallel chains
%
ChainsStr = get(controls.posterior.chains,'String');
NumChains = str2num(StringTrim(ChainsStr(DSGEModel.ParallelChainsValue,:)));
%
% 4. 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);
   drawnow;
   return;
end;
%
% 5. Scale factor for initializing the random walk Metropolis algorithm, i.e., c0 in "YADA Manual - Computational Details".
%
InitScaleStr = get(controls.posterior.initialscale,'String');
c0 = str2num(StringTrim(InitScaleStr(DSGEModel.MHInitialScaleFactor,:)));
%
% 6. Scale factor for the random walk Metropolis algorithm, i.e., c in "YADA Manual - Computational Details".
%
ScaleStr = get(controls.posterior.scalefactor,'String');
c = str2num(StringTrim(ScaleStr(DSGEModel.MHScaleFactor,:)));
%
% 7. Randomization weight for choosing initial values (multiple chains only)
%    If 0, then posterior mode is used, if 1 then random draws are used.
%
RWStr = get(controls.posterior.randomweight,'String');
RWvalue = str2num(StringTrim(RWStr(DSGEModel.RandomWeightValue,:)));
RWStr = '';
if NumChains>1;
   RWStr = [' * Weight on randomization:            ' num2str(RWvalue)];
end;
%
% make sure the output directory is exists
%
stat = MakeDir(DSGEModel.OutputDirectory,'rwm');
if stat~=1;
   SetLogText(strvcat(' ','Couldn''t create the directory:',['  ' DSGEModel.OutputDirectory '\rwm'],'Posterior sampling was aborted.'));
   txt = ['YADA was for some reason unable to create the directory "' DSGEModel.OutputDirectory '\rwm". Posterior sampling has therefore been aborted.'];
   About(txt,'information','YADA - Directory Creation Problem',180,500,CurrINI);
   drawnow;
   return;
end;
%
SetLogText(strvcat(' ',[' * Time:                               ' StringTrim(datestr(now,14)) ', ' GetCurrentDate],[' * Draws per chain:                    ' num2str(NumPostDraws,'%0.0f')],[' * Sample batches per chain:           ' num2str(NumBatches,'%0.0f')],[' * Number of parallel chains:          ' num2str(NumChains,'%0.0f')],[' * Burn-in draws per chain:            ' num2str(NumBurnin,'%0.0f')],[' * Inverse Hessian estimator:          ' InvHessianTxt],[' * Maximum absolute correlation:       ' num2str(MaxCorr)],[' * Scale factor for initial phi value: ' num2str(c0)],[' * Scale factor in proposal density:   ' num2str(c)],RWStr,[' * Number of parameters to estimate:   ' num2str(NumParam)],[' * Total number of parameters:         ' num2str(size(ParameterNames.all,1),'%0.0f')],[' * Selected sample:                    ' DSGEModel.SubBeginYear ':' DSGEModel.SubBeginPeriod ' - ' DSGEModel.SubEndYear ':' DSGEModel.SubEndPeriod],[' * Initializing Kalman filter:         ' num2str(DSGEModel.KalmanFirstObservation-1,'%0.0f') ' observations'],[' * Effective sample size:              ' num2str(DSGEModel.LastPeriod-DSGEModel.FirstPeriod+1-(DSGEModel.KalmanFirstObservation-1),'%0.0f')],[' * Estimation sample:                  ' NewStartYear ':' NewStartPeriod ' - ' DSGEModel.SubEndYear ':' DSGEModel.SubEndPeriod],[' * Number of observed variables:       ' num2str(DSGEModel.n,'%0.0f')],[' * Number of exogenous variables:      ' num2str(DSGEModel.k,'%0.0f')],[' * Number of state variables:          ' num2str(length(DSGEModel.StateVariablePositions),'%0.0f')],[' * Number of state shocks:             ' num2str(size(DSGEModel.StateShockNames,1),'%0.0f')],[' * Number of measurement errors:       ' num2str(rank(R),'%0.0f')],' '));
%
% Check if we should keep old draws
%
StartChain = 1;
NBatchDraws = 1;
BatchIndex = (1:NumBatchDraws);
StartBatch = 1;
StartPostDraws = 1;
PostSample = zeros(NumPostDraws,NumParam);
thetaPostSample = PostSample;
LogPost = zeros(NumPostDraws,1);
AcceptedDraws = zeros(NumPostDraws,1);
LoadedOldData = 0;
if DSGEModel.OverwriteDraws==0;
   %
   % check which files may already have been stored.
   %
   ExitLoop = 0;
   for NChains=1:NumChains;
      BatchIndex = (1:NumBatchDraws);
      %
      % Initialize matrices to store parameters etc in.
      %
      PostSample = zeros(NumPostDraws,NumParam);
      thetaPostSample = PostSample;
      LogPost = zeros(NumPostDraws,1);
      AcceptedDraws = zeros(NumPostDraws,1);
      for NBatch=1:NumBatches;
         PostSampleFile = [DSGEModel.OutputDirectory '\rwm\RWM-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-' num2str(NChains,'%0.0f') '.' num2str(NumChains,'%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
            %
            RWMData = load(PostSampleFile);
            PostSample(BatchIndex,:) = RWMData.BatchPostSample;
            thetaPostSample(BatchIndex,:) = RWMData.BatchthetaPostSample;
            LogPost(BatchIndex,1) = RWMData.BatchLogPost;
            AcceptedDraws(BatchIndex,1) = RWMData.BatchAcceptedDraws;
            phi = RWMData.BatchPostSample(size(RWMData.BatchPostSample,1),:)';
            theta = RWMData.BatchthetaPostSample(size(RWMData.BatchthetaPostSample,1),:)';
            if NBatch==1;
               if isfield(RWMData,'phiInit')==1;
                  RWM(NChains).phiInit = RWMData.phiInit;
               else;
                  RWM(NChains).phiInit = RWMData.BatchPostSample(1,:)';
               end;
               phiInit = RWM(NChains).phiInit;
            end;
            %
            BatchIndex = BatchIndex+NumBatchDraws;
            LoadedOldData = 1;
            clear('RWMData');
         else;
            ExitLoop = 1;
            break;
         end;
      end;
      if (ExitLoop==1)|(NBatch==NumBatches);
         StartChain = NChains;
         StartBatch = NBatch;
         StartPostDraws = max(BatchIndex)-NumBatchDraws+1;
         %
         % if a full chain has been found at the end add 1 to StartChain
         %
         if isequal(sum(abs(PostSample(NumPostDraws,:))),0)==0;
            StartChain = NChains+1;
            StartBatch = 1;
            StartPostDraws = 1;
            BatchIndex = (1:NumBatchDraws);
            if StartChain<=NumChains;
               PostSample = zeros(NumPostDraws,NumParam);
               thetaPostSample = PostSample;
               LogPost = zeros(NumPostDraws,1);
               AcceptedDraws = zeros(NumPostDraws,1);
            end;
         end;
         if ExitLoop==1;
            break;
         end;
      end;
   end;
end;
%
if ((DSGEModel.OverwriteDraws==1)|(LoadedOldData==0))&(NumChains==1)&(DSGEModel.PosteriorDrawsValue>1);
   %
   % Check if we should try to find data to load
   % 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 '\rwm\RWM-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-1.1-' num2str(NBatch,'%0.0f') '.' num2str(NBatches,'%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);
      drawnow;
      if strcmp(lower(answer),'yes')==1;
         %
         % first we load the data
         %
         NTempBatchDraws = NPDraws/NBatches;
         %
         BatchInd = (1:NTempBatchDraws);
         %
         for i=1:NBatch;
            PostSampleFile = [DSGEModel.OutputDirectory '\rwm\RWM-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-1.1-' num2str(i,'%0.0f') '.' num2str(NBatches,'%0.0f') '.' num2str(NPDraws,'%0.0f') '.mat'];
            %
            RWMData = load(PostSampleFile);
            PostSample(BatchInd,:) = RWMData.BatchPostSample;
            thetaPostSample(BatchInd,:) = RWMData.BatchthetaPostSample;
            LogPost(BatchInd,1) = RWMData.BatchLogPost;
            AcceptedDraws(BatchInd,1) = RWMData.BatchAcceptedDraws;
            BatchInd = BatchInd+NTempBatchDraws;
         end;
         %
         StartPostDraws = ceil(NPDraws*(NBatch/NBatches))+1;
         phi = PostSample(StartPostDraws-1,:)';
         theta = thetaPostSample(StartPostDraws-1,:)';
         if isfield(RWMData,'phiInit')==1;
            RWM(NumChains).phiInit = RWMData.phiInit;
         else;
            RWM(NumChains).phiInit = PostSample(1,:)';
         end;
         phiInit = RWM(NumChains).phiInit;
         clear('RWMData');
         %
         % 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 = LogPost(BatchIndex,1);
               if sum(BatchLogPost~=0)==NumBatchDraws;
                  %
                  % we have a full sample
                  %
                  PostSampleFile = [DSGEModel.OutputDirectory '\rwm\RWM-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-1.1-' num2str(i,'%0.0f') '.' num2str(NumBatches,'%0.0f') '.' num2str(NumPostDraws,'%0.0f') '.mat'];
                  BatchPostSample = PostSample(BatchIndex,:);
                  BatchthetaPostSample = thetaPostSample(BatchIndex,:);
                  BatchLogPost = LogPost(BatchIndex,1);
                  BatchAcceptedDraws = AcceptedDraws(BatchIndex,1);
                  %
                  % save to file
                  %
                  SetLogText(strvcat(['   Saving:           ' GetFilename(PostSampleFile)],['   Acceptance ratio: ' num2str(AcceptedDraws(BatchIndex(length(BatchIndex)),1)/BatchIndex(length(BatchIndex)))]));
                  save(PostSampleFile,'BatchPostSample','BatchthetaPostSample','BatchLogPost','BatchAcceptedDraws','phiInit');
                  drawnow;
                  pause(0.02);
                  %
                  % clear batched matrices from memeory
                  %
                  clear('BatchPostSample','BatchthetaPostSample','BatchLogPost','BatchAcceptedDraws');
                  %
                  BatchIndex = BatchIndex+NumBatchDraws;
               else;
                  NBatchDraws = sum(BatchLogPost~=0)+1;
                  StartBatch = i;
                  break;
               end;
            end;
         else;
            NBatchDraws = StartPostDraws;
         end;
      end;
   end;
elseif (LoadedOldData==0)&(NumChains>1);
   %
   % check if a lower number of chains can be located and draws obtained from them
   % The same number of batches is required
   %
   NPDrawsValue = DSGEModel.PosteriorDrawsValue;
   PostSampleFile = '';
   %
   NPDraws = str2num(StringTrim(strrep(DrawsStr(NPDrawsValue,:),',','')));
   NBatches = NumBatches;
   if NBatches>NPDraws;
      NBatches = NPDraws;
   end;
   %
   % check if output exists
   %
   BreakChain = 0;
   UseChain = 0;
   TotalChains = 0;
   for NChains=1:NumChains-1;
      for ThisChain=1:NChains;
         for NBatch=1:NBatches;
            PostSampleFileI = [DSGEModel.OutputDirectory '\rwm\RWM-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-' num2str(ThisChain,'%0.0f') '.' num2str(NChains,'%0.0f') '-' num2str(NBatch,'%0.0f') '.' num2str(NBatches,'%0.0f') '.' num2str(NPDraws,'%0.0f') '.mat'];
            if FileExist(PostSampleFileI)==1;
               PostSampleFile = PostSampleFileI;
               UseChain = ThisChain;
               TotalChains = NChains;
            else;
               NBatch = NBatch-1;
               BreakChain = 1;
               break;
            end;
         end;
      end;
   end;
   %
   % UseChain TotalChains NPDrawsValue
   %
   if isempty(PostSampleFile)==0;
      txt = ['YADA has located ' num2str(ceil(UseChain*NPDraws*(NBatch/NBatches)),'%0.0f') ' previous posterior draws when running ' num2str(UseChain,'%0.0f') ' out of ' num2str(TotalChains,'%0.0f') ' random walk Metropolis chains. Would you like to include these in your ' num2str(NumChains,'%0.0f') ' chains with ' num2str(NumPostDraws,'%0.0f') ' draws?'];
      answer = Query(txt,'question',150,'Posterior Sampling',500,'no',CurrINI);
      drawnow;
      if strcmp(lower(answer),'yes')==1;
         for NChain=1:UseChain;
            if (NChain<UseChain)|(NBatch==NumBatches);
               %
               % we can just copy the files
               %
               for NB=1:NBatches;
                  OldFile = [DSGEModel.OutputDirectory '\rwm\RWM-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-' num2str(NChain,'%0.0f') '.' num2str(TotalChains,'%0.0f') '-' num2str(NB,'%0.0f') '.' num2str(NBatches,'%0.0f') '.' num2str(NPDraws,'%0.0f') '.mat'];
                  NewFile = [DSGEModel.OutputDirectory '\rwm\RWM-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-' num2str(NChain,'%0.0f') '.' num2str(NumChains,'%0.0f') '-' num2str(NB,'%0.0f') '.' num2str(NBatches,'%0.0f') '.' num2str(NPDraws,'%0.0f') '.mat'];
                  [status,message] = CopyFile(OldFile,NewFile);
                  RWMData = load(OldFile);
                  AcceptedDraws = RWMData.BatchAcceptedDraws;
                  if NB==1;
                     RWMData = load(OldFile);
                     if isfield(RWMData,'phiInit')==1;
                        RWM(NChain).phiInit = RWMData.phiInit;
                     else;
                        RWM(NChain).phiInit = RWMData.BatchPostSample(1,:)';
                     end;
                  end;
                  SetLogText(strvcat(['   Saving:           ' GetFilename(NewFile)],['   Acceptance ratio: ' num2str(AcceptedDraws(length(AcceptedDraws))/(length(AcceptedDraws)*NB))]));
                  drawnow;
                  pause(0.02);
                  clear('RWMData');
               end;
               %
               % Write some final results to the estimation log window
               %
               SetLogText(strvcat(' ',['Random Walk Metropolis sampling chain no. ' num2str(NChain,'%0.0f') ' completed.'],[' * Time:                                  ' StringTrim(datestr(now,14)) ', ' GetCurrentDate],' * Posterior sampling time:               00:00:00.00 (hrs:min:sec)',[' * Acceptance ratio:                      ' num2str(AcceptedDraws(length(AcceptedDraws),1)/NumPostDraws)],' '));
               StartChain = StartChain+1;
            else;
               NTempBatchDraws = NPDraws/NBatches;
               %
               BatchInd = (1:NTempBatchDraws);
               %
               % we copy until NBatch, but we shall also load the data
               %
               for NB=1:NBatch;
                  OldFile = [DSGEModel.OutputDirectory '\rwm\RWM-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-' num2str(NChain,'%0.0f') '.' num2str(TotalChains,'%0.0f') '-' num2str(NB,'%0.0f') '.' num2str(NBatches,'%0.0f') '.' num2str(NPDraws,'%0.0f') '.mat'];
                  NewFile = [DSGEModel.OutputDirectory '\rwm\RWM-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-' num2str(NChain,'%0.0f') '.' num2str(NumChains,'%0.0f') '-' num2str(NB,'%0.0f') '.' num2str(NBatches,'%0.0f') '.' num2str(NPDraws,'%0.0f') '.mat'];
                  [status,message] = CopyFile(OldFile,NewFile);
                  %
                  % now we load the data
                  %
                  RWMData = load(OldFile);
                  PostSample(BatchInd,:) = RWMData.BatchPostSample;
                  thetaPostSample(BatchInd,:) = RWMData.BatchthetaPostSample;
                  LogPost(BatchInd,1) = RWMData.BatchLogPost;
                  AcceptedDraws(BatchInd,1) = RWMData.BatchAcceptedDraws;
                  BatchInd = BatchInd+NTempBatchDraws;
                  SetLogText(strvcat(['   Saving:           ' GetFilename(NewFile)],['   Acceptance ratio: ' num2str(AcceptedDraws(BatchInd(length(BatchInd)))/BatchInd(length(BatchInd)))]));
                  drawnow;
                  pause(0.02);
               end;
               %
               StartPostDraws = ceil(NPDraws*(NBatch/NBatches))+1;
               phi = PostSample(StartPostDraws-1,:)';
               theta = thetaPostSample(StartPostDraws-1,:)';
               if isfield(RWMData,'phiInit')==1;
                  RWM(NChain).phiInit = RWMData.phiInit;
               else;
                  RWM(NChain).phiInit = PostSample(1,:)';
               end;
               phiInit = RWM(NChain).phiInit;
               %
               BatchIndex = BatchInd;
               StartBatch = NBatch+1;
               StartChain = UseChain;
               clear('RWMData');
            end;
         end;
      end;
   end;
end;
%
% Now we are ready to start the chains
%
abort = '';
CholeskiSigma = chol(SigmaMode)';
if StartChain>NumChains;
   %
   % Everything has already been gone through
   %
   SetLogText('Posterior sampling with random walk Metropolis already completed!');
   txt = 'The random walk Metropolis 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);
   drawnow;
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 (StartPostDraws-1+(NumPostDraws*(StartChain-1)))>0;
      for i=1:StartChain;
         if c0>0;
            randn(NumParam,1);
         end;
         if (NumChains>1)&(RWvalue>0);
            randn(NumParam,1);
         end;
         if i<StartChain;
            ATemp = rand(1,NumPostDraws);
            BTemp = randn(NumParam,NumPostDraws);
         else;
            if StartPostDraws>1;
               ATemp = rand(1,StartPostDraws-1);
               BTemp = randn(NumParam,StartPostDraws-1);
            end;
         end;
      end;
      clear('ATemp','BTemp');
   end;
else;
   %
   % we set the state based on the current time
   %
   rand('state',sum(100*clock));
   randn('state',sum(100*clock));
end;
%
% initialize the PostSampleTime variable
%
tic;
PostSampleTime = toc;
%
for NChains=StartChain:NumChains;
   if NChains>StartChain;
      %
      % reset overything that needs resetting
      %
      BatchIndex = (1:NumBatchDraws);
      NBatchDraws = 1;
      NBatch = 1;
      %
      % Initialize matrices to store parameters etc in.
      %
      PostSample = zeros(NumPostDraws,NumParam);
      thetaPostSample = PostSample;
      LogPost = zeros(NumPostDraws,1);
      AcceptedDraws = zeros(NumPostDraws,1);
      StartPostDraws = 1;
   else;
      NBatch = StartBatch;
   end;
   %
   % Initialize phi if it doesn't already exist
   %
   if StartPostDraws==1;
      if c0>0;
         phi = phiMode+(c0*CholeskiSigma*randn(NumParam,1));
      else;
         phi = phiMode;
      end;
      if NumChains>1;
         if RWvalue>0;
            phi = phiMode+(4*RWvalue*CholeskiSigma*randn(NumParam,1));
         else;
            phi = phiMode;
         end;
      end;
      RWM(NChains).phiInit = phi;
      phiInit = phi;
      theta = PhiToTheta(phi,thetaIndex,UniformBounds,LowerBound);
   end;
   %
   % Now evaluate the log posterior at the starting value
   %
   logPostInit = -logPosteriorPhiDSGE(phi,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
   logPostOld = logPostInit;
   %
   % 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 = 'Acceptance Ratio:';
      %
      ProgressStructure.name = ['Random Walk Metropolis - Chain ' num2str(NChains,'%0.0f') ' of ' num2str(NumChains,'%0.0f')];
      ProgressStructure.CurrINI = CurrINI;
      mywait = ProgressDLG(0,ProgressStructure);
      set(mywait,'Color',get(0,'defaultuicontrolbackgroundcolor'));
      drawnow;
   else;
      %
      % show a wait dialog
      %
      txt = ['Please wait while YADA runs the random walk Metropolis sampler for ' int2str(NChains) ' of ' int2str(NumChains) ' Markov chains. The total number of posterior draws per chain is ' int2str(NumPostDraws) ', while the number of burn-in draws per chain is ' int2str(NumBurnin) '. The current chain was started at: ' StringTrim(datestr(now,14)) ' on ' GetCurrentDate '.'];
      WaitHandle = WaitDLG(txt,'information',['Random Walk Metropolis - Chain ' num2str(NChains,'%0.0f') ' of ' num2str(NumChains,'%0.0f')],500,200,CurrINI,0);
      WaitControls = get(WaitHandle,'UserData');
   end;
   %
   tic;
   for NPostDraws=StartPostDraws:NumPostDraws;
      %
      % draw a proposal phi
      %
      phiProp = phi+(c*CholeskiSigma*randn(NumParam,1));
      %
      % Evaluate the log posterior at phiProp
      %
      logPostProp = -logPosteriorPhiDSGE(phiProp,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
      AcceptProb = min(1,exp(logPostProp-logPostOld));
      if rand<AcceptProb;
         if NPostDraws>1;
            AcceptedDraws(NPostDraws,1) = AcceptedDraws(NPostDraws-1,1)+1;
         else;
            AcceptedDraws(1,1) = 1;
         end;
         phi = phiProp;
         logPostOld = logPostProp;
         theta = PhiToTheta(phi,thetaIndex,UniformBounds,LowerBound);
      else;
         if NPostDraws>1;
            AcceptedDraws(NPostDraws,1) = AcceptedDraws(NPostDraws-1,1);
         else;
            AcceptedDraws(1,1) = 0;
         end;
      end;
      %
      % 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 AcceptedDraws(NPostDraws,1)/NPostDraws]);
            end;
         end;
      end;
      %
      % store draws
      %
      PostSample(NPostDraws,:) = phi';
      thetaPostSample(NPostDraws,:) = theta';
      LogPost(NPostDraws,1) = logPostOld;
      %
      if NBatchDraws<NumBatchDraws;
         %
         % increase
         %
         NBatchDraws = NBatchDraws+1;
      else;
         %
         % save a batch
         %
         PostSampleFile = [DSGEModel.OutputDirectory '\rwm\RWM-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-' num2str(NChains,'%0.0f') '.' num2str(NumChains,'%0.0f') '-' num2str(NBatch,'%0.0f') '.' num2str(NumBatches,'%0.0f') '.' num2str(NumPostDraws,'%0.0f') '.mat'];
         BatchPostSample = PostSample(BatchIndex,:);
         BatchthetaPostSample = thetaPostSample(BatchIndex,:);
         BatchLogPost = LogPost(BatchIndex,1);
         BatchAcceptedDraws = AcceptedDraws(BatchIndex,1);
         %
         SetLogText(strvcat(['   Saving:           ' GetFilename(PostSampleFile)],['   Acceptance ratio: ' num2str(AcceptedDraws(NPostDraws,1)/NPostDraws)]));
         save(PostSampleFile,'BatchPostSample','BatchthetaPostSample','BatchLogPost','BatchAcceptedDraws','phiInit');
         drawnow;
         pause(0.02);
         %
         % clear batched matrices from memory
         %
         clear('BatchPostSample','BatchthetaPostSample','BatchLogPost','BatchAcceptedDraws');
         %
         NBatch = NBatch+1;
         NBatchDraws = 1;
         BatchIndex = BatchIndex+NumBatchDraws;
      end;
      drawnow;
   end;
   PostSampleTime = 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;
   if strcmp(abort,'cancel')==0;
      if NumChains==1;
         SetLogText(strvcat(' ','Posterior sampling with the random walk Metropolis completed!',[' * Time:                               ' StringTrim(datestr(now,14)) ', ' GetCurrentDate],[' * Posterior sampling time:            ' SecondsToHoursStr(PostSampleTime) ' (hrs:min:sec)'],[' * Acceptance ratio:                   ' num2str(AcceptedDraws(NumPostDraws,1)/NumPostDraws)]));
      elseif NChains<NumChains;
         SetLogText(strvcat(' ',['Random Walk Metropolis sampling chain no. ' num2str(NChains,'%0.0f') ' completed.'],[' * Time:                                  ' StringTrim(datestr(now,14)) ', ' GetCurrentDate],[' * Posterior sampling time:               ' SecondsToHoursStr(PostSampleTime) ' (hrs:min:sec)'],[' * Acceptance ratio:                      ' num2str(AcceptedDraws(NumPostDraws,1)/NumPostDraws)],' '));
      else;
         SetLogText(strvcat(' ',['Random Walk Metropolis sampling chain no. ' num2str(NChains,'%0.0f') ' completed.'],[' * Time:                                  ' StringTrim(datestr(now,14)) ', ' GetCurrentDate],[' * Posterior sampling time:               ' SecondsToHoursStr(PostSampleTime) ' (hrs:min:sec)'],[' * Acceptance ratio:                      ' num2str(AcceptedDraws(NumPostDraws,1)/NumPostDraws)]));
      end;
      drawnow;
      pause(0.02);
   else;
      SetLogText('Posterior sampling aborted.');
      break;
   end;
end;
%
% clean up the tmp directory for AiM related error data.
%
if FileExist([pwd '\tmp\AiMErrorData.mat'])==1;
   dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
end;
%
if strcmp(abort,'cancel')==1;
   return;
end;
%
% check if we have just run 1 chain. if so display certain information.
%
for CurrChain=1:NumChains;
   %
   % make sure that we load the correct results
   %
   if NumChains>1;
      %
      % load the data
      %
      SetLogText(strvcat(' ','Reloading posterior draws data...'));
      [PostSample,thetaPostSample,LogPost,AcceptedDraws,NumBurnin] = LoadDSGEPosteriorDraws(DSGEModel,CurrINI,CurrChain);
      %
      % temporarily remove sequential estimation
      %
      DSGEModel.SequentialML = 0;
      SetLogText(strvcat('Computing summary information about the parameters...',[' * Chain number:                       ' num2str(CurrChain,'%0.0f')],[' * Time:                               ' StringTrim(datestr(now,14))]));
   else;
      SetLogText(strvcat(' ','Computing summary information about the parameters...',[' * Time:                               ' StringTrim(datestr(now,14))]));
   end;
   %
   NumAcceptedDraws = AcceptedDraws(NumPostDraws,1)/NumPostDraws;
   clear('AcceptedDraws');
   %
   % compute summary information about the parameter draws for the original parameters
   %
   tic;
   ParData = ParameterSummary(thetaPostSample(NumBurnin+1:NumPostDraws,:));
   ParTime = toc;
   SetLogText(strvcat('Computation completed!',[' * Estimation time:                    ' SecondsToHoursStr(ParTime) ' (hrs:min:sec)']));
   %
   % check if we have parameters to update
   %
   if FileExist(DSGEModel.UpdateParameterFile)==1;
      tic;
      SetLogText(strvcat(' ','Preparing posterior values for parameters to update...',[' * Time:                               ' StringTrim(datestr(now,14))]));
      [UpdatePostSample,UpdateMode,UpdateNames] = GetUpdatePosteriorDraws(thetaPostSample,NumBurnin,CurrChain,NumChains,0,DSGEModel,CurrINI);
      SetLogText(strvcat('Computing summary information...',[' * Time:                               ' StringTrim(datestr(now,14))]));
      UpdateParData = ParameterSummary(UpdatePostSample);
      UpParTime = toc;
      SetLogText(strvcat('Computation completed!',[' * Total preparation time:             ' SecondsToHoursStr(UpParTime) ' (hrs:min:sec)']));
   else;
      UpdatePostSample = [];
      UpdateParData = [];
      UpdateMode = [];
      UpdateNames = '';
   end;
   %
   % try to estimate the posterior mode of the marginal distribution.
   %
   if DSGEModel.PosteriorDensityValue==1;
      ModeAlgorithm = ' * Algorithm:                          Gaussian';
   elseif DSGEModel.PosteriorDensityValue==2;
      ModeAlgorithm = strvcat(' * Algorithm:                          Silverman-type and','                                       Sköld-Roberts correction');
   elseif DSGEModel.PosteriorDensityValue==3;
      ModeAlgorithm = ' * Algorithm:                          Sheather-Jones bandwidth';
   elseif DSGEModel.PosteriorDensityValue==4;
      ModeAlgorithm = ' * Algorithm:                          Bump killing bandwidth';
   end;
   SetLogText(strvcat(' ','Computing the posterior mode of the marginal distribution...',ModeAlgorithm,[' * Time:                               ' StringTrim(datestr(now,14))]));
   tic;
   [thetaMargMode,UpdateMargMode] = GraphPosteriorDensity(thetaPostSample,NumBurnin,CurrChain,UpdatePostSample,UpdateMode,UpdateNames,[],[],[],DSGEModel,CurrINI);
   ModeTime = toc;
   if isempty(thetaMargMode)==1;
      thetaMargMode = thetaMode;
      UpdateMargMode = UpdateMode;
      SetLogText(strvcat('Estimation aborted.',' * Using the posterior mode of the joint distribution.'));
   else;
      SetLogText(strvcat('Estimation completed!',[' * Estimation time:                    ' SecondsToHoursStr(ModeTime) ' (hrs:min:sec)']));
   end;
   %
   % clear original parameter draws from memory
   %
   clear('thetaPostSample','UpdatePostSample');
   %
   if DSGEModel.MarginalLikelihoodValue~=4;
      %
      % compute the marginal likelihood
      %
      SetLogText(strvcat(' ','Computing the marginal likelihood...'));
      if DSGEModel.MarginalLikelihoodValue~=3;
         %
         % calculate the marginal likelihood with the modified harmonic mean estimator
         %
         SetLogText(strvcat(' * Geweke''s modified harmonic mean estimator',[' * Time:                               ' StringTrim(datestr(now,14))]));
         tic;
         [LogMargsMH,CovProb] = MargLikeModifiedHarmonic(PostSample(NumBurnin+1:NumPostDraws,:),LogPost(NumBurnin+1:NumPostDraws,:),DSGEModel.SequentialML,[],DSGEModel);
         ModHarmTime = toc;
         SetLogText(strvcat('Estimation completed!',[' * Estimation time:                    ' SecondsToHoursStr(ModHarmTime) ' (hrs:min:sec)']));
      else;
         LogMargsMH = [];
         CovProb = [];
      end;
      if DSGEModel.MarginalLikelihoodValue~=2;
         %
         % calculate the marginal likelihood with the Chib and Jeliazkov estimator
         %
         SetLogText(strvcat(' * Chib and Jeliazkov''s estimator',[' * Time:                               ' StringTrim(datestr(now,14))]));
         tic;
         LogMargsCJ = MargLikeChibJeliazkov(PostSample,LogPost,phiMode,SigmaMode,c,logPostMode,NumBurnin,DSGEModel.SequentialML,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
         CJTime = toc;
         SetLogText(strvcat('Estimation completed!',[' * Estimation time:                    ' SecondsToHoursStr(CJTime) ' (hrs:min:sec)']));
      else;
         LogMargsCJ = [];
      end;
      %
      % clear draws of transformed parameters from memory
      %
      clear('PostSample','LogPost');
      %
      % prepare plots
      %
      GraphDirectory = [DSGEModel.OutputDirectory '\graphs'];
      if (DSGEModel.MarginalLikelihoodValue~=3)&(DSGEModel.SequentialML==1);
         SetLogText(strvcat(' ','Creating plot of sequential modified harmonic mean estimate...'));
         FigHandle =figure('Visible','off', ...
                'Position',[(CurrINI.scrsz(3)-700)/2 (CurrINI.scrsz(4)-500)/2 700 500], ...
                'Tag','MargLikeMHM', ...
                'Name','Sequential Marginal Likelihood - Modified Harmonic Mean');
         %
         SetFigureProperties(FigHandle,CurrINI);
         %
         Xvalues = LogMargsMH(:,1)+NumBurnin;
         Yvalues = LogMargsMH(:,2:size(LogMargsMH,2));
         CurrMap = colormap;
         NumColor = 1;
         NumStyle = 1;
         ColorOrder = [(1:8:57) (2:8:58) (3:8:59) (4:8:60) (5:8:61) (6:8:62) (7:8:63) (8:8:64)];
         LegendStr = '';
         for i=1:size(Yvalues,2);
            if NumColor<=length(CurrMap);
               CurrColor = CurrMap(ColorOrder(i),:);
            else;
               CurrColor = CurrMap(ColorOrder(1),:);
               NumColor = 1;
            end;
            if NumStyle==1;
               CurrStyle = '-';
            elseif NumStyle==2;
               CurrStyle = '--';
            elseif NumStyle==3;
               CurrStyle = ':';
            elseif NumStyle==4;
               NumStyle = 0;
               CurrStyle = '-.';
            end;
            plot(Xvalues,Yvalues(:,i),'LineStyle',CurrStyle,'Color',CurrColor,'LineWidth',1);
            hold('on');
            NumColor = NumColor+1;
            NumStyle = NumStyle+1;
            LegendStr = strvcat(LegendStr,[num2str(CovProb(i)*100,'%0.0f') ' percent coverage prob']);
         end;
         %
         % plot the average
         %
         plot(Xvalues,(1/size(Yvalues,2))*sum(Yvalues,2),'LineStyle','-','Color',[0 0 0],'LineWidth',1);
         LegendStr = strvcat(LegendStr,'Overall mean');
         hold('off');
         %
         SetAxesFonts(CurrINI,'Title','Modified Harmonic Mean Estimates','XLabel',['RWM 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);
         %
         legHandle = legend(LegendStr,0);
         SetAxesFonts(CurrINI,legHandle);
         %
         % save the grasph
         %
         if DirExist(GraphDirectory);
            file = [GraphDirectory '\LogMarginalLikelihood-ModifiedHarmonic-' num2str(NumPostDraws,'%0.0f') '-' int2str(CurrChain) '.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;
         pause(0.02);
      end;
      if (DSGEModel.MarginalLikelihoodValue~=2)&(DSGEModel.SequentialML==1);
         SetLogText(strvcat(' ','Creating plot of sequential Chib and Jeliazkov estimate...'));
         FigHandle =figure('Visible','off', ...
                'Position',[(CurrINI.scrsz(3)-700)/2 (CurrINI.scrsz(4)-500)/2 700 500], ...
                'Tag','MargLikeCJ', ...
                'Name','Sequential Marginal Likelihood - Chib and Jeliazkov');
         %
         SetFigureProperties(FigHandle,CurrINI);
         %
         Xvalues = LogMargsCJ(:,1)+NumBurnin;
         Yvalues = LogMargsCJ(:,2);
         %
         plot(Xvalues,Yvalues,'LineStyle','-','Color',[0.75 0 0],'LineWidth',2,'Marker','.');
         hold('off');
         SetAxesFonts(CurrINI,'Title','Chib and Jeliazkov Estimates','XLabel',['RWM 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-ChibJeliazkov-' num2str(NumPostDraws,'%0.0f') '-' int2str(CurrChain) '.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;
         pause(0.02);
      end;
   else;
      LogMargsMH = [];
      CovProb = [];
      LogMargsCJ = [];
   end;
   %
   % write sample statistics to file
   %
   SetLogText(strvcat(' ','Writing posterior sampling summary results to text file:',['   Directory: ' DSGEModel.OutputDirectory],['   File name: PosteriorSamplingSummary-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-' num2str(CurrChain,'%0.0f') '.' num2str(NumChains,'%0.0f') '-' num2str(NumPostDraws,'%0.0f') '.txt']));
   PrintPosteriorSamplingSummary(DSGEModel,CurrINI,ParameterNames,PriorDist,thetaDist,ParData,thetaMargMode,UpdateParData,UpdateMargMode,UpdateNames,NumPostDraws,NumBurnin,rank(R),c0,c,PostSampleTime,NumAcceptedDraws,LogMargsMH,CovProb,LogMargsCJ,StepLength,CurrChain,NumChains,PhiToTheta(RWM(CurrChain).phiInit,thetaIndex,UniformBounds,LowerBound),RWvalue,MaxCorr);
   PostSampleTextFile = [DSGEModel.OutputDirectory '\PosteriorSamplingSummary-' DSGEModel.NameOfModel '-' num2str(DSGEModel.InverseHessianEstimator,'%0.0f') '-' num2str(CurrChain,'%0.0f') '.' num2str(NumChains,'%0.0f') '-' num2str(NumPostDraws,'%0.0f') '.txt'];
   clear('LogMargsMH','LogMargsCJ');
   %
   % display the results
   %
   if CurrChain==NumChains;
      TextGUI(PostSampleTextFile,['Posterior Sampling - ' DSGEModel.NameOfModel],[(CurrINI.scrsz(3)-1000)/2 32 1000 CurrINI.scrsz(4)-100],'Results',50,CurrINI,1,1);
   end;
end;

%
% end of DSGEPosteriorSampling.m
%
