function PosteriorModeEstimation(DSGEModel,CurrINI,controls)
% PosteriorModeEstimation: Runs the posterior mode estimation procedures.
%
% USAGE:
%
%       PosteriorModeEstimation(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: August 10, 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:
%
% * 11-08-2006: Added the check optimum features.
%
% * 14-08-2006: Added some initial features regarding additional and calibrated parameters.
%
% * 15-08-2006: Added reading the sub-sample setting and the "CurrINI.GridWidth" option (when checking the optimum)
%
% * 16-08-2006: Removed an unncessary file exist check before copying the measurement equation file to the tmp directory.
%
% * 17-08-2006: Added the two new inputs for the measurement equation file with string matrices containing the names
%               of the observed variables and of the exogenous variables.
%
% * 18-08-2006: Changed the StrTrim function to StringStrim.
%
% * 21-08-2006: Added automatic saving of graphs to "GraphDirectory".
%
% * 22-08-2006: Added calculations of the largest eigenvalue of F, the state transition matrix, at the mode and at the
%               initial values.
%
% * 23-08-2006: Added collection of results from the optimization routines through the iterations. Forecast errors are
%               now also computed the observed variables.
%
% * 24-08-2006: Removed the output "InitialValues" from the call to "VerifyPriorData". Furthermore, an unnecessary call
%               to the measurement equation function was removed.
%
% * 25-08-2006: Added changes required to handle a uniform prior.
%
% * 28-08-2006: Added the input CloseSelf=1 to the TextGUI call. This allows the latter function to avoid using the waitfor
%               command.
%
% * 14-09-2006: Added the new inputs for the measurement equation call.
%
% * 15-09-2006: Checking the ModelParameters structure. No changes made!
%
% * 21-09-2006: Checked the PriorDist.normal matrix. No changes made! Bug fixed in VerifyPriorData function.
%
% * 22-09-2006: The AiMSolver call is now run twice to figure out the solution time.
%
% * 16-10-2006: Collected information about the initialized parameters names in order to separate them from the
%               additional parameters that are updated.
%
% * 17-10-2006: Added a try, catch call when attempting to load mat-files. This is needed if matlab can't
%               read the mat-file, e.g., because it was create with a later version of matlab.
%
% * 20-10-2006: Added the call to the summary output from the estimation.
%
% * 25-10-2006: The field KalmanFirstObservation to the DSGEModel structure now uses the round function instead of the
%               ceil function.
%
% * 26-10-2006: Added the DSGEModel fields "FirstPeriod" and "LastPeriod" to the data stored on disk from the posterior
%               mode estimation.
%
% * 09-11-2006: Added information sent to output window about the estimation sample. The function also checks the gradient
%               of the log posterior at the initial values to capture potential bad gradients.
%
% * 10-11-2006: The value of the log prior and the log Jacobian at the initial value are now also sent to the output window.
%
% * 20-11-2006: The initial value for the state vector is now user defined through the variable DSGEModel.InitialStateVector.
%
% * 22-11-2006: Additional output from failed optimizations are now written to file.
%
% * 24-11-2006: Took the name changes "logPosteriorPhiDSGEOptimize" to "logPosteriorPhiDSGE", and "logPosteriorThetaDSGEOptimize"
%               to "logPosteriorThetaDSGE" into account. The "phiEstimates" matrix is now obtained from the optimizer and stored
%               in the output file.
%
% * 06-12-2006: Made sure that the "CalculateDSGEStateVariables" function receives the new input variable.
%
% * 07-12-2006: Took the changes is fields for CurrINI and DSGEModel into account, where CurrINI.MaximumIterationsValue
%               is now given by DSGEModel.OptMaxIterationsValue, CurrINI.ToleranceValue is now DSGEModel.OptToleranceValue.
%               Moreover, CurrINI.MaximizeAlgorithmValue is now DSGEModel.MaximizeAlgorithmValue; CurrINI.StepLengthHessian is
%               DSGEModel.StepLengthHessian; CurrINI.NumberOfGridPoints is DSGEModel.NumberOfGridPoints;
%               CurrINI.FiniteDifferenceHessian is DSGEModel.FiniteDifferenceHessian; CurrINI.CheckOptimum is
%               DSGEModel.CheckOptimum; CurrINI.CheckTransformedOptimum is DSGEModel.CheckTransformedOptimum; CurrINI.GridWidth
%               is DSGEModel.GridWidth.
%
% * 20-02-2007: The matrix with partial derivatives "ThetaPhiPartial", i.e., d Theta/ d phi', is now saved to the posterior
%               mode results file.
%
% * 09-03-2007: Made sure that "CurrINI" is sent to the "CalculateDSGEMaximumEigenvalue" function.
%
% * 20-04-2007: Added the input "phiFittedStd" for the call to the "PrintPosteriorModeSummary" function.
%
% * 12-11-2007: Updated the documentation.
%
% * 11-12-2007: Introduced the variable "NumParam" that replaced calls to the length function for theta, phi,
%               thetaMode and phiMode.
%
% * 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.
%
% * 16-05-2008: The function now attempts to collect data from AiM output from the tmp directory for error
%               message reporting. Took the "DSGEModel.AIMTolerance" input for AiMSolver into account.
%
% * 19-05-2008: Made sure that the "AiMErrorData.mat" file is cleaned from the tmp directory after estimation
%               has finished.
%
% * 20-05-2008: The gradient vector is now sent to the "PrintOptimizationErrorReport" function.
%
% * 23-05-2008: Updated the documentation.
%
% * 27-05-2008: The mcode value 7 is now supported.
%
% * 03-07-2008: Made sure that the measurement equation receives only the active state variable names.
%               Added try-catch code for AiMSolver.
%
% * 22-07-2008: Took the lower bound for the gamma prior into account.
%
% * 23-07-2008: Took the change from 0-1 lower and upper bounds to general bounds for the beta distribution
%               into account.
%
% * 23-10-2008: The method for computing the total number of parameters have been changed.
%
% * 21-01-2009: Updated the documentation.
%
% * 30-01-2009: Added code for dealing with user choice of method to initialize the inverse Hessian
%               when using the csminwel optimization routine.
%
% * 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: The mcode values 7 and 8 are now properly mentioned when encountered.
%
% * 09-03-2009: Added quotation marks around paths for the dos del and copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

EstimationStatus = 0;
if FileExist([pwd '\tmp\AiMErrorData.mat'])==1;
   dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
end;
%
% Check the file with prior information
%
SetLogText(strvcat(' ','Initializing posterior mode estimation...',[' * Checking "' GetFilename(DSGEModel.PriorFile) '" with data on prior...']));
%
[status,theta,thetaIndex,thetaDist,LowerBound,ModelParameters,thetaPositions,PriorDist,ParameterNames,UniformBounds] = VerifyPriorData(DSGEModel,CurrINI);
if (status==1)&(isempty(theta)==0);
   OriginalModelParameters = ModelParameters;
   %
   % Transform PriorDist.gamma and PriorDist.beta to the relevant parameters that are used by the logGammaPDF and
   % the logBetaPDF functions
   %
   if isempty(PriorDist.beta)==0;
      [a,b] = MomentToParamStdbetaPDF((PriorDist.beta(:,1)-PriorDist.beta(:,3))./(PriorDist.beta(:,4)-PriorDist.beta(:,3)),PriorDist.beta(:,2)./(PriorDist.beta(:,4)-PriorDist.beta(:,3)));
      PriorDist.beta_ab = [a b];
   else;
      PriorDist.beta_ab = [];
   end;
   %
   if isempty(PriorDist.gamma)==0;
      [a,b] = MomentToParamGammaPDF(PriorDist.gamma(:,1)-PriorDist.gamma(:,3),PriorDist.gamma(:,2));
      PriorDist.gamma_ab = [a b];
   else;
      PriorDist.gamma_ab = [];
   end;
   %
   % 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;
   if FileExist(DSGEModel.InitializeParameterFile)==1;
      [stat,msg] = dos(['copy /Y "' DSGEModel.InitializeParameterFile '" "' pwd '\tmp"']);
   end;
   ErrorStr = '';
   ErrorFile = '';
   try;
      AIMData = load(DSGEModel.AIMDataFile);
   catch;
      ErrorStr = ['Unable to load the file "' DSGEModel.AIMDataFile '". Message caught is: ' lasterr];
      ErrorFile = DSGEModel.AIMDataFile;
   end;
   %
   % The structure AIMData contains fields such as neq, nlead and nlag
   %
   firstUpdatePos = 0;
   lastUpdatePos = 0;
   firstInitializePos = 0;
   lastInitializePos = 0;
   if isempty(ErrorStr)==1;
      if DSGEModel.RunInitializeFirst==1;
         %
         % check if the file with additional parameters that only needs to be updated is OK
         %
         SetLogText(strvcat('   STATUS: OK',[' * Checking ' GetName(DSGEModel.InitializeParameterFile) '...']));
         %
         % we need to run the file DSGEModel.InitializeParameterFile first
         %
         try;
            firstInitializePos = length(fieldnames(ModelParameters))+1;
            eval(['ModelParameters = ' GetName(DSGEModel.InitializeParameterFile) '(ModelParameters);']);
            lastInitializePos = length(fieldnames(ModelParameters));
         catch;
            ErrorStr = ['Error found in "' GetName(DSGEModel.InitializeParameterFile) '.m". Message caught is: ' lasterr];
            ErrorFile = DSGEModel.InitializeParameterFile;
         end;
         if (FileExist(DSGEModel.UpdateParameterFile)==1)&(isempty(ErrorStr)==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;
      else;
         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;
         if (FileExist(DSGEModel.InitializeParameterFile)==1)&(isempty(ErrorStr)==1);
            SetLogText(strvcat('   STATUS: OK',[' * Checking ' GetName(DSGEModel.InitializeParameterFile) '...']));
            try;
               firstInitializePos = length(fieldnames(ModelParameters))+1;
               eval(['ModelParameters = ' GetName(DSGEModel.InitializeParameterFile) '(ModelParameters);']);
               lastInitializePos = length(fieldnames(ModelParameters));
            catch;
               ErrorStr = ['Error found in "' GetName(DSGEModel.InitializeParameterFile) '.m". Message caught is: ' lasterr];
               ErrorFile = DSGEModel.InitializeParameterFile;
            end;
         end;
      end;
   end;
   %
   if isempty(ErrorStr)==1;
      %
      % we can go on with checking if the initial values permit a unique solution of the DSGE model
      %
      SetLogText(strvcat('   STATUS: OK',' * Checking if the DSGE model has a unique convergent','   solution for initial parameter values...'));
      %
      % we can now run the AiM solver:
      %
      % To measure the solution time more accurately we call it twice, with the first call only being
      % performed for checking if the model has a unique convergent solution at the initial values.
      %
      try;
         [mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
      catch;
         ErrorStr = ['There was an error running AiM. YADA received the following message: ' lasterr];
         About(ErrorStr,'error','AiM Error',300,500,CurrINI);
         drawnow;
         return;
      end;
      if mcode==1;
         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)==0);
            %
            % check the measurement equation file for internal validity
            %
            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;
               SetLogText(strvcat(['   STATUS: Bad input file ' GetName(ErrorFile)],' '));
               About(ErrorStr,'error',['Error - ' GetFilename(ErrorFile)],200,500,CurrINI);
               drawnow;
               return;
            end;
            %
            % set the InitialStateVector value
            %
            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;
            %
            % The model has a unique solution for the initial values
            %
            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);
            %
            StepLengthStr = get(controls.optimize.steplength,'String');
            StepLength = str2num(StringTrim(StepLengthStr(DSGEModel.StepLengthHessian,:)));
            %
            NumGridPointsStr = get(controls.optimize.numbergridpoints,'String');
            NumGridPoints = str2num(StringTrim(NumGridPointsStr(DSGEModel.NumberOfGridPoints,:)));
            %
            SetLogText(strvcat('   STATUS: OK',' * Checking if the log posterior can be calculated...'));
            %
            % Collect additional parameters in ParameterNames.additional
            %
            ParameterNames = ReadAdditionalParameterNames(ParameterNames,ModelParameters,OriginalModelParameters,[firstUpdatePos lastUpdatePos firstInitializePos lastInitializePos]);
            %
            % set up vectors with calibrated parameters and additional parameters. These are only used for printing to the output file
            %
            [thetaCalibrated,thetaAdditionalInitial] = ReadAdditionalParameterValues(ModelParameters,ParameterNames);
            %
            % Set up the sub-sample
            %
            [DSGEModel.FirstPeriod,DSGEModel.LastPeriod] = CreateSubSample(DSGEModel);
            %
            % Determine the first observation to use when computing log-likelihood with the Kalman filter
            %
            [NewStartYear,NewStartPeriod] = AdjustSampleStart(DSGEModel.SubBeginYear,DSGEModel.SubBeginPeriod,DSGEModel.DataFrequency,DSGEModel.KalmanFirstObservation-1);
            %
            NumParam = length(theta);
            %
            ReturnCodeText = '';
            %
            % Starting the optimization
            %
            if DSGEModel.MaximizeAlgorithmValue==1;
               %
               % We use Sims' csminwel routine with transformed parameters.
               %
               OptAlgorithm = 'csminwel';
               phi = ThetaToPhi(theta,thetaIndex,UniformBounds,LowerBound);
               ErrorStr = '';
               try;
                  tic;
                  logPost = logPosteriorPhiDSGE(phi,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                  LogPostTime = toc;
                  LogPosteriorPhiIV = -logPost;
               catch;
                  ErrorStr = ['Unable to calculate the log posterior. Message caught is: ' lasterr];
               end;
               %
               if ((exist('logPost')==0)|(logPost==1000000))&(isempty(ErrorStr)==1);
                  SetLogText(strvcat('   STATUS: Unable to compute log posterior at initial values',' '));
                  txt = 'Unable to compute the log posterior at the initial values. You probably need to change your initial parameter values.';
                  About(txt,'information','YADA - DSGE Model Problem',200,500,CurrINI);
                  drawnow;
                  return;
               else;
                  %
                  % Compute a numerical gradient
                  %
                  [g, badg] = YADAnumgrad('logPosteriorPhiDSGE',phi,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                  if badg~=0;
                     SetLogText(strvcat('   STATUS: Bad gradient of log posterior at initial values',' '));
                     [badpar,maxpos] = min(abs(g));
                     %
                     if badpar>0;
                        ValueStr = num2str(badpar);
                     else;
                        ValueStr = 'NaN';
                     end;
                     %
                     ParamStr = ['"' StringTrim(ParameterNames.estimated(maxpos(1),:)) '"'];
                     AddParamStr = ' ';
                     if length(maxpos)>=2;
                        AddParamStr = 's ';
                        for i=2:length(maxpos);
                           ParamStr = [ParamStr ', "' StringTrim(ParameterNames.estimated(maxpos(i),:)) '"'];
                        end;
                     end;
                     %
                     txt = ['Unable to compute a meaningful numerical gradient of the log posterior at the initial values. The largest absolute value or an improper value (NaN or Inf) of the gradient is found for the parameter' AddParamStr ParamStr ', with value of ' ValueStr '. You probably need to change your initial parameter values.'];
                     About(txt,'information','YADA - DSGE Model Problem',200,500,CurrINI);
                     drawnow;
                     return;
                  end;
               end;
               %
               if isempty(ErrorStr)==1;
                  LogJacobianIV = logJacobian(phi,thetaIndex,UniformBounds);
                  LogPriorIV = logPriorDSGE(PhiToTheta(phi,thetaIndex,UniformBounds,LowerBound),thetaDist,PriorDist,LowerBound);
                  LogPosteriorThetaIV = LogPosteriorPhiIV-LogJacobianIV;
                  LogLikelihoodIV = LogPosteriorThetaIV-LogPriorIV;
                  %
                  % compute the initial inverse Hessian
                  %
                  if DSGEModel.InitializeHessian==1;
                     InitialInvHessMat = 10*eye(NumParam);
                  elseif DSGEModel.InitializeHessian==2;
                     InitialInvHessMat = 1*eye(NumParam);
                  elseif DSGEModel.InitializeHessian==3;
                     InitialInvHessMat = 10^(-2)*eye(NumParam);
                  elseif DSGEModel.InitializeHessian==4;
                     InitialInvHessMat = 10^(-3)*eye(NumParam);
                  elseif DSGEModel.InitializeHessian==5;
                     InitialInvHessMat = 10^(-4)*eye(NumParam);
                  elseif DSGEModel.InitializeHessian==6;
                     InitialInvHessMat = 10^(-5)*eye(NumParam);
                  elseif DSGEModel.InitializeHessian==7;
                     InitialInvHessMat = 10^(-6)*eye(NumParam);
                  else;
                     %
                     % sample from the prior
                     %
                     InitialInvHessMat = PriorHessianMatrix(thetaDist,thetaIndex,PriorDist,UniformBounds,LowerBound);
                  end;
                  %
                  % Run the measurement equation function just so that we can determine the number of measurement errors
                  %
                  if FileExist([pwd '\tmp\AiMErrorData.mat'])==1;
                     dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
                  end;
                  SetLogText(strvcat('   STATUS: OK',['   Computation time:       ' num2str(LogPostTime) ' seconds'],['   Log posterior value:    ' num2str(LogPosteriorPhiIV,'%0.4f')],['   Log prior value:        ' num2str(LogPriorIV,'%0.4f')],['   Log Jacobian value:     ' num2str(LogJacobianIV,'%0.4f')],['   Log likelihood value:   ' num2str(LogLikelihoodIV,'%0.4f')],' ','Optimization routine:','   Sims'' csminwel for transformed parameters...',[' * Time:                             ' StringTrim(datestr(now,14))],[' * Number of parameters to estimate: ' num2str(NumParam)],[' * Total number of parameters:       ' num2str(size(ParameterNames.all,1),'%0.0f')],[' * Tolerance level:                  ' ToleranceStr],[' * Maximum number of iterations:     ' MaxIterStr],[' * 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')]));
                  tic;
                  [fh,phiMode,gh,InverseHessian,NumIter,fcount,FuncDiff,AbortEstimation,OptResults,phiEstimates,ReturnCode] = YADAcsminwel('logPosteriorPhiDSGE',phi,InitialInvHessMat,[],ToleranceValue,MaxIter,DSGEModel.ShowProgress,DSGEModel.ShowProgressClock,'Posterior Mode Estimation - csminwel',CurrINI, ...
                     thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                  MaxTime = toc;
                  GradientVector = abs(gh);
                  %
                  % gh is the numerical gradient, fh the function value at minimum
                  %
                  % ReturnCode: 0 = minimum found
                  %             1 = zero gradient
                  %        2 or 4 = back and forth on step length never finished
                  %             3 = smallest step still improving too slow
                  %             5 = largest step still improving too fast
                  %             6 = smallest step still improving too slow, reversed gradient
                  %             7 = warning: possible inaccuracy in H matrix
                  %
                  if (ReturnCode==0)&(strcmp(lower(AbortEstimation),'cancel')~=1)&(NumIter<=MaxIter);
                     LogPostDensity = -fh;
                     %
                     % we're OK
                     %
                     % Calculate the log posterior at theta
                     %
                     thetaMode = PhiToTheta(phiMode,thetaIndex,UniformBounds,LowerBound);
                     %
                     LogJacobianMode = logJacobian(phiMode,thetaIndex,UniformBounds);
                     LogPriorMode = logPriorDSGE(thetaMode,thetaDist,PriorDist,LowerBound);
                     LogLikelihoodMode = LogPostDensity-(LogPriorMode+LogJacobianMode);
                     %
                     SetLogText(strvcat(' ','Optimization completed!',[' * Time:                             ' StringTrim(datestr(now,14))],[' * Optimization time:                ' SecondsToHoursStr(MaxTime) ' (hrs:min:sec)'],[' * Number of iterations:             ' num2str(NumIter)],[' * Convergence value:                ' num2str(FuncDiff)],[' * Log posterior of phi at mode:     ' num2str(LogPostDensity,'%0.4f')],[' * Log posterior of theta at mode:   ' num2str(LogPostDensity-LogJacobianMode,'%0.4f')],[' * Log prior at mode:                ' num2str(LogPriorMode,'%0.4f')],[' * Log Jacobian at mode:             ' num2str(LogJacobianMode,'%0.4f')],[' * Log likelihood at mode:           ' num2str(LogLikelihoodMode,'%0.4f')]));
                     %
                     EstimationStatus = 1;
                     %
                  elseif strcmp(lower(AbortEstimation),'cancel')==1;
                     %
                     % posteriod mode estimation was aborted
                     %
                     SetLogText(strvcat(' ','Posterior mode estimation was aborted.'));
                  elseif (ReturnCode==0)&(strcmp(lower(AbortEstimation),'cancel')~=1)&(NumIter>MaxIter);
                     %
                     % maximum number of iterations was reached
                     %
                     txt = 'The maximum number of iterations was reached without csminwel converging. YADA has aborted from posterior mode estimation.';
                     SetLogText(strvcat(' ','Maximum number of iteration reached.'));
                     About(txt,'information','Optimization Problem',200,500,CurrINI);
                     drawnow;
                  elseif ReturnCode==1;
                     txt = ['Csminwel return code: zero gradient. YADA has aborted from posterior mode estimation. Writing optimization information to file: "OptErrorReport-' DSGEModel.NameOfModel '.txt" in "' DSGEModel.OutputDirectory '".'];
                     SetLogText(strvcat(' ','Csminwel return code: zero gradient.'));
                     About(txt,'information','Optimization Problem',200,500,CurrINI);
                     drawnow;
                     ReturnCodeText = 'Csminwel return code: zeros gradient.';
                  elseif (ReturnCode==2)|(ReturnCode==4);
                     txt = ['Csminwel return code: back and forth on step length never finished. YADA has aborted from posterior mode estimation. Writing optimization information to file: "OptErrorReport-' DSGEModel.NameOfModel '.txt" in "' DSGEModel.OutputDirectory '".'];
                     SetLogText(strvcat(' ','Csminwel return code: back and forth on step length never finished.'));
                     About(txt,'information','Optimization Problem',200,500,CurrINI);
                     drawnow;
                     ReturnCodeText = 'Csminwel return code: back and forth on step length never finished.';
                  elseif ReturnCode==3;
                     txt = ['Csminwel return code: smallest step still improving too slow. YADA has aborted from posterior mode estimation. Writing optimization information to file: "OptErrorReport-' DSGEModel.NameOfModel '.txt" in "' DSGEModel.OutputDirectory '".'];
                     SetLogText(strvcat(' ','Csminwel return code: smallest step still improving too slow.'));
                     About(txt,'information','Optimization Problem',200,500,CurrINI);
                     drawnow;
                     ReturnCodeText = 'Csminwel return code: smallest step still improving too slow.';
                  elseif ReturnCode==5;
                     txt = ['Csminwel return code: largest step still improving too fast. YADA has aborted from posterior mode estimation. Writing optimization information to file: "OptErrorReport-' DSGEModel.NameOfModel '.txt" in "' DSGEModel.OutputDirectory '".'];
                     SetLogText(strvcat(' ','Csminwel return code: largest step still improving too fast.'));
                     About(txt,'information','Optimization Problem',200,500,CurrINI);
                     drawnow;
                     ReturnCodeText = 'Csminwel return code: largest step still improving too fast.';
                  elseif ReturnCode==6;
                     txt = ['Csminwel return code: smallest step still improving too slow, reversed gradient. YADA has aborted from posterior mode estimation. Writing optimization information to file: "OptErrorReport-' DSGEModel.NameOfModel '.txt" in "' DSGEModel.OutputDirectory '".'];
                     SetLogText(strvcat(' ','Csminwel return code: smallest step still improving too slow, reversed gradient.'));
                     About(txt,'information','Optimization Problem',200,500,CurrINI);
                     drawnow;
                     ReturnCodeText = 'Csminwel return code: smallest step still improving too slow, reversed gradient.';
                  elseif ReturnCode==7;
                     txt = ['Csminwel return code: possible inaccuracy in inverse Hessian. YADA has aborted from posterior mode estimation. Writing optimization information to file: "OptErrorReport-' DSGEModel.NameOfModel '.txt" in "' DSGEModel.OutputDirectory '".'];
                     SetLogText(strvcat(' ','Csminwel return code: possible inaccuracy in inverse Hessian.'));
                     About(txt,'information','Optimization Problem',200,500,CurrINI);
                     drawnow;
                     ReturnCodeText = 'Csminwel return code: possible inaccuracy in inverse Hessian.';
                  end;
               else;
                  SetLogText(strvcat('   STATUS: Uanble to calculate log posterior.',' '));
                  About(ErrorStr,'error','Error - Log Posterior',200,500,CurrINI);
                  drawnow;
               end;
            elseif DSGEModel.MaximizeAlgorithmValue==2;
               phi = ThetaToPhi(theta,thetaIndex,UniformBounds,LowerBound);
               %
               % use matlab's fminunc
               %
               OptAlgorithm = 'fminunc';
               ErrorStr = '';
               try;
                  tic;
                  logPost = logPosteriorPhiDSGE(phi,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                  LogPosteriorPhiIV = -logPost;
                  LogPostTime = toc;
               catch;
                  ErrorStr = ['Unable to calculate the log posterior. Message caught is: ' lasterr];
               end;
               %
               if ((exist('logPost')==0)|(logPost==1000000))&(isempty(ErrorStr)==1);
                  SetLogText(strvcat('   STATUS: Unable to compute log posterior at initial values',' '));
                  txt = 'Unable to compute the log posterior at the initial values. You need to change your initial parameter values.';
                  About(txt,'information','YADA - DSGE Model Problem',200,500,CurrINI);
                  drawnow;
                  return;
               else;
                  %
                  % Compute a numerical gradient
                  %
                  [g, badg] = YADAnumgrad('logPosteriorPhiDSGE',phi,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                  if badg~=0;
                     SetLogText(strvcat('   STATUS: Bad gradient of log posterior at initial values',' '));
                     [badpar,maxpos] = max(abs(g));
                     %
                     txt = ['Unable to compute a meaningful numerical gradient of the log posterior at the initial values. The largest absolute value of the gradient is found for the parameter "' StringTrim(ParameterNames.estimated(maxpos,:)) '", with an gradient value of ' num2str(badpar) '. You probably need to change your initial parameter values.'];
                     About(txt,'information','YADA - DSGE Model Problem',200,500,CurrINI);
                     drawnow;
                     return;
                  end;
               end;
               %
               if isempty(ErrorStr)==1;
                  LogJacobianIV = logJacobian(phi,thetaIndex,UniformBounds);
                  LogPriorIV = logPriorDSGE(PhiToTheta(phi,thetaIndex,UniformBounds,LowerBound),thetaDist,PriorDist,LowerBound);
                  LogPosteriorThetaIV = LogPosteriorPhiIV-LogJacobianIV;
                  LogLikelihoodIV = LogPosteriorThetaIV-LogPriorIV;
                  %
                  % set the options
                  %
                  Options = optimset('MaxIter',MaxIter,'MaxFunEvals',100*MaxIter,'Diagnostics','off','LargeScale','off','TolFun',ToleranceValue,'TolX',ToleranceValue);
                  %
                  % call the function
                  %
                  if FileExist([pwd '\tmp\AiMErrorData.mat'])==1;
                     dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
                  end;
                  SetLogText(strvcat('   STATUS: OK',['   Computation time:       ' num2str(LogPostTime) ' seconds'],['   Log posterior value:    ' num2str(LogPosteriorPhiIV,'%0.4f')],['   Log prior value:        ' num2str(LogPriorIV,'%0.4f')],['   Log Jacobian value:     ' num2str(LogJacobianIV,'%0.4f')],['   Log likelihood value:   ' num2str(LogLikelihoodIV,'%0.4f')],' ','Optimization routine:','   Matlab''s fminunc for transformed parameters...',[' * Time:                             ' StringTrim(datestr(now,14))],[' * Number of parameters to estimate: ' num2str(NumParam)],[' * Total number of parameters:       ' num2str(size(ParameterNames.all,1),'%0.0f')],[' * Tolerance level:                  ' ToleranceStr],[' * Maximum number of iterations:     ' MaxIterStr],[' * 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')]));
                  if MatlabNumber>=7;
                     tic;
                     [phiMode,fh,ExitFlag,OUTPUT,Grad,HessianMatrix] = YADAfminunc7('logPosteriorPhiDSGE',phi,Options,DSGEModel.ShowProgress,DSGEModel.ShowProgressClock,'Posterior Mode Estimation - fminunc',CurrINI, ...
                        thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                     MaxTime = toc;
                     GradientVector = abs(Grad);
                     AbortEstimation = OUTPUT.abort;
                     NumIter = OUTPUT.iterations;
                     OptResults = OUTPUT.OptResults;
                     if (ExitFlag>0)&(strcmp(lower(AbortEstimation),'cancel')~=1);
                        FuncDiff = OUTPUT.convergence;
                        LogPostDensity = -fh;
                        InverseHessian = inv(HessianMatrix);
                        phiEstimates = OUTPUT.phiEstimates;
                        %
                        % we're OK
                        %
                        % Calculate the log posterior at theta
                        %
                        thetaMode = PhiToTheta(phiMode,thetaIndex,UniformBounds,LowerBound);
                        %
                        LogJacobianMode = logJacobian(phiMode,thetaIndex,UniformBounds);
                        LogPriorMode = logPriorDSGE(thetaMode,thetaDist,PriorDist,LowerBound);
                        LogLikelihoodMode = LogPostDensity-(LogPriorMode+LogJacobianMode);
                        %
                        SetLogText(strvcat(' ','Optimization completed!',[' * Time:                             ' StringTrim(datestr(now,14))],[' * Optimization time:                ' SecondsToHoursStr(MaxTime) ' (hrs:min:sec)'],[' * Number of iterations:             ' num2str(NumIter)],[' * Convergence value:                ' num2str(FuncDiff)],[' * Log posterior of phi at mode:     ' num2str(LogPostDensity,'%0.4f')],[' * Log posterior of theta at mode:   ' num2str(LogPostDensity-LogJacobianMode,'%0.4f')],[' * Log prior at mode:                ' num2str(LogPriorMode,'%0.4f')],[' * Log Jacobian at mode:             ' num2str(LogJacobianMode,'%0.4f')],[' * Log likelihood at mode:           ' num2str(LogLikelihoodMode,'%0.4f')]));
                        %
                        EstimationStatus = 1;
                        %
                     elseif strcmp(lower(AbortEstimation),'cancel')==1;
                        %
                        % posteriod mode estimation was aborted
                        %
                        SetLogText(strvcat(' ','Posterior mode estimation was aborted.'));
                     elseif ExitFlag==0;
                        %
                        % Reached maximum number of iteration
                        %
                        txt = 'The maximum number of function evaluations or iterations was reached without fminunc converging (ExitFlag is 0). YADA has aborted from posterior mode estimation.';
                        SetLogText(strvcat(' ','Maximum number of function evaluations or iteration reached.',[' * Time:                             ' StringTrim(datestr(now,14))],[' * Optimization time:                ' SecondsToHoursStr(MaxTime) ' (hrs:min:sec)']));
                        About(txt,'information','Optimization Problem',200,500,CurrINI);
                        drawnow;
                        ReturnCodeText = 'Fminunc return code: the maximum number of function evaluations or iteration was reached without fminunc coverging (ExitFlag is 0).';
                     elseif ExitFlag<0;
                        %
                        % no convergence
                        %
                        txt = ['Fminunc didn''t converge (ExitFlag is ' int2str(ExitFlag) '). YADA has aborted from posterior mode estimation. Writing optimization information to file: "OptErrorReport-' DSGEModel.NameOfModel '.txt" in "' DSGEModel.OutputDirectory '".'];
                        SetLogText(strvcat(' ',['Fminunc didn''t converge (ExitFlag is ' int2str(ExitFlag) ').'],[' * Time:                             ' StringTrim(datestr(now,14))],[' * Optimization time:                ' SecondsToHoursStr(MaxTime) ' (hrs:min:sec)']));
                        About(txt,'information','Optimization Problem',200,500,CurrINI);
                        drawnow;
                        ReturnCodeText = ['Fminunc return code: No covergence (ExitFlag is ' int2str(ExitFlag) ').'];
                     end;
                  else;
                     tic;
                     [phiMode,fh,ExitFlag,OUTPUT,Grad,HessianMatrix] = YADAfminunc5('logPosteriorPhiDSGE',phi,Options,DSGEModel.ShowProgress,DSGEModel.ShowProgressClock,'Posterior Mode Estimation - fminunc',CurrINI, ...
                        thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                     MaxTime = toc;
                     AbortEstimation = OUTPUT.abort;
                     NumIter = OUTPUT.iterations;
                     OptResults = OUTPUT.OptResults;
                     GradientVector = abs(Grad);
                     if (ExitFlag>0)&(strcmp(lower(AbortEstimation),'cancel')~=1);
                        FuncDiff = OUTPUT.convergence;
                        LogPostDensity = -fh;
                        InverseHessian = inv(HessianMatrix);
                        phiEstimates = OUTPUT.phiEstimates;
                        %
                        % we're OK
                        %
                        % Calculate the log posterior at theta
                        %
                        thetaMode = PhiToTheta(phiMode,thetaIndex,UniformBounds,LowerBound);
                        %
                        LogJacobianMode = logJacobian(phiMode,thetaIndex,UniformBounds);
                        LogPriorMode = logPriorDSGE(thetaMode,thetaDist,PriorDist,LowerBound);
                        LogLikelihoodMode = LogPostDensity-(LogPriorMode+LogJacobianMode);
                        %
                        SetLogText(strvcat(' ','Optimization completed!',[' * Time:                             ' StringTrim(datestr(now,14))],[' * Optimization time:                ' SecondsToHoursStr(MaxTime) ' (hrs:min:sec)'],[' * Number of iterations:             ' num2str(NumIter)],[' * Convergence value:                ' num2str(FuncDiff)],[' * Log posterior of phi at mode:     ' num2str(LogPostDensity,'%0.4f')],[' * Log posterior of theta at mode:   ' num2str(LogPostDensity-LogJacobianMode,'%0.4f')],[' * Log prior at mode:                ' num2str(LogPriorMode,'%0.4f')],[' * Log Jacobian at mode:             ' num2str(LogJacobianMode,'%0.4f')],[' * Log likelihood at mode:           ' num2str(LogLikelihoodMode,'%0.4f')]));
                        %
                        EstimationStatus = 1;
                        %
                     elseif strcmp(lower(AbortEstimation),'cancel')==1;
                        %
                        % posteriod mode estimation was aborted
                        %
                        SetLogText(strvcat(' ','Posterior mode estimation was aborted.'));
                     elseif ExitFlag==0;
                        %
                        % Reached maximum number of iteration
                        %
                        txt = 'The maximum number of function evaluations or iterations was reached without fminunc converging (ExitFlag is 0). YADA has aborted from posterior mode estimation.';
                        SetLogText(strvcat(' ','Maximum number of function evaluations or iteration reached.',[' * Time:                             ' StringTrim(datestr(now,14))],[' * Optimization time:                ' SecondsToHoursStr(MaxTime) ' (hrs:min:sec)']));
                        About(txt,'information','Optimization Problem',200,500,CurrINI);
                        drawnow;
                        ReturnCodeText = 'Fminunc return code: the maximum number of function evaluations or iteration was reached without fminunc coverging (ExitFlag is 0).';
                     elseif ExitFlag<0;
                        %
                        % no convergence
                        %
                        txt = ['Fminunc didn''t converge (ExitFlag is ' int2str(ExitFlag) '). YADA has aborted from posterior mode estimation.'];
                        SetLogText(strvcat(' ',['Fminunc didn''t converge (ExitFlag is ' int2str(ExitFlag) ').'],[' * Time:                             ' StringTrim(datestr(now,14))],[' * Optimization time:                ' SecondsToHoursStr(MaxTime) ' (hrs:min:sec)']));
                        About(txt,'information','Optimization Problem',200,500,CurrINI);
                        drawnow;
                        ReturnCodeText = ['Fminunc return code: No covergence (ExitFlag is ' int2str(ExitFlag) ').'];
                     end;
                  end;
               else;
                  SetLogText(strvcat('   STATUS: Uanble to calculate log posterior.',' '));
                  About(ErrorStr,'error','Error - Log Posterior',200,500,CurrINI);
                  drawnow;
               end;
            end;
            %
            % we can continue from here given that estimation went well
            %
            if EstimationStatus==1;
               if FileExist([pwd '\tmp\AiMErrorData.mat'])==1;
                  dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
               end;
               %
               % Record the first observation used by the Kalman filter when calculating the log-likelihood
               %
               KalmanFirstObservation = DSGEModel.KalmanFirstObservation;
               FirstPeriod = DSGEModel.FirstPeriod;
               LastPeriod = DSGEModel.LastPeriod;
               %
               % Compute the largest eigenvalue of F in the state equation at the mode and at the initial values
               %
               MaxEigenvalueIV = CalculateDSGEMaximumEigenvalue(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI);
               MaxEigenvalueMode = CalculateDSGEMaximumEigenvalue(thetaMode,thetaPositions,ModelParameters,DSGEModel,CurrINI);
               %
               % Compute the partial derivatives of the original parameters wrt the transformed parameters, evaluated at the
               % posterior mode
               %
               ThetaPhiPartial = PartialThetaPartialPhi(phiMode,thetaIndex,UniformBounds);
               %
               % Check if we should compute the inverse Hessian with finite difference methods
               %
               if DSGEModel.FiniteDifferenceHessian==1;
                  SetLogText(strvcat(' ','Computing posterior standard deviations from inverse Hessian','using finite differences...'));
                  %
                  % for the transformed parameters
                  %
                  VechHessian = YADAHessianMatrix('logPosteriorPhiDSGE',phiMode,StepLength,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                  HessianMat = GeneralizedCholeski(reshape(VechHessian,NumParam,NumParam));
                  phiInverseHessianFinDiff = inv(HessianMat);
                  phiPostModeStdFinDiff = sqrt(diag(phiInverseHessianFinDiff));
                  %
                  % for the orginal parameters
                  %
                  VechHessian = YADAHessianMatrix('logPosteriorThetaDSGE',thetaMode,StepLength,thetaIndex,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                  HessianMat = GeneralizedCholeski(reshape(VechHessian,NumParam,NumParam));
                  thetaInverseHessianFinDiff = inv(HessianMat);
                  thetaPostModeStdFinDiff = sqrt(diag(thetaInverseHessianFinDiff));
               else;
                  phiInverseHessianFinDiff = NaN;
                  phiPostModeStdFinDiff = NaN;
                  thetaInverseHessianFinDiff = NaN;
                  thetaPostModeStdFinDiff = NaN;
               end;
               %
               % Check if we should check the optimum
               %
               if DSGEModel.CheckOptimum==1;
                  SetLogText(strvcat(' ','Checking the optimum...',' * Transformed parameters'));
                  %
                  % Check if GraphDirectory exists
                  %
                  GraphDirectory = [DSGEModel.OutputDirectory '\graphs'];
                  MinMaxGrid = zeros(NumParam,2);
                  stat = MakeDir(DSGEModel.OutputDirectory,'graphs');
                  %
                  % create lower and upper bounds from InverseHessian and phiMode
                  %
                  MinMaxGrid(:,1) = phiMode-(DSGEModel.GridWidth*sqrt(diag(InverseHessian)));
                  MinMaxGrid(:,2) = phiMode+(DSGEModel.GridWidth*sqrt(diag(InverseHessian)));
                  phiFittedStd = VillaniCheckOptimum('logPosteriorPhiDSGE',phiMode,MinMaxGrid,InverseHessian,NumGridPoints,0,[],ParameterNames.estimated,'Checking Curvature at Posterior Mode - Transformed Parameters',CurrINI,GraphDirectory,'phi', ...
                     thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                  if DSGEModel.CheckTransformedOptimum==0;
                     SetLogText(' * Original parameters');
                     MinMaxGrid(:,1) = PhiToTheta(MinMaxGrid(:,1),thetaIndex,UniformBounds,LowerBound);
                     MinMaxGrid(:,2) = PhiToTheta(MinMaxGrid(:,2),thetaIndex,UniformBounds,LowerBound);
                     thetaFittedStd = VillaniCheckOptimum('logPosteriorThetaDSGE',thetaMode,MinMaxGrid,ThetaPhiPartial*InverseHessian*ThetaPhiPartial,NumGridPoints,0,[],ParameterNames.estimated,'Checking Curvature at Posterior Mode - Original Parameters',CurrINI,GraphDirectory,'theta', ...
                        thetaIndex,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
                  else;
                     thetaFittedStd = NaN;
                  end;
               else;
                  phiFittedStd = NaN;
                  thetaFittedStd = NaN;
               end;
               if FileExist([pwd '\tmp\AiMErrorData.mat'])==1;
                  dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
               end;
               %
               % check the parameters with most extreme correlations based on the inverse Hessian
               %
               ExtremeCorrBound = 0.7;
               CorrHessian = diag(1./sqrt(diag(InverseHessian)))*InverseHessian*diag(1./sqrt(diag(InverseHessian)));
               CorrHessian = tril(CorrHessian,-1);
               [RowIndex,ColIndex] = find(abs(CorrHessian)>=ExtremeCorrBound);
               if isempty(RowIndex)==0;
                  ExtremeCorrParameters = cell(length(RowIndex),3);
                  for i=1:length(RowIndex)
                     ExtremeCorrParameters{i,1} = ParameterNames.estimated(RowIndex(i),:);
                     ExtremeCorrParameters{i,2} = ParameterNames.estimated(ColIndex(i),:);
                     ExtremeCorrParameters{i,3} = CorrHessian(RowIndex(i),ColIndex(i));
                  end;
               else;
                  ExtremeCorrParameters = [];
               end;
               %
               % calculate the log-likelihood at theta and the log prior
               %
               ModelParameters = ThetaToModelParameters(ModelParameters,thetaMode,thetaPositions);
               if FileExist(DSGEModel.UpdateParameterFile)==1;
                  eval(['ModelParameters = ' GetName(DSGEModel.UpdateParameterFile) '(ModelParameters);']);
               end;
               [LogLikeDensity,mcode,status] = logLikelihoodDSGE(ModelParameters,DSGEModel,AIMData);
               LogPriorDensity = logPriorDSGE(thetaMode,thetaDist,PriorDist,LowerBound);
               LogJacobianValue = logJacobian(phiMode,thetaIndex,UniformBounds);
               LogPriorDensityIV = logPriorDSGE(theta,thetaDist,PriorDist,LowerBound);
               %
               % Compute forecasts of the observed variables at mode and at initial values
               %
               StateVarStructureMode = CalculateDSGEStateVariables(thetaMode,thetaPositions,ModelParameters,DSGEModel,0);
               StateVarStructureIV = CalculateDSGEStateVariables(theta,thetaPositions,ModelParameters,DSGEModel,0);
               T = length(StateVarStructureMode.Y);
               MSEYMode = (1/T)*(StateVarStructureMode.Y-StateVarStructureMode.Yhat)*(StateVarStructureMode.Y-StateVarStructureMode.Yhat)';
               MSEYIV = (1/T)*(StateVarStructureIV.Y-StateVarStructureIV.Yhat)*(StateVarStructureIV.Y-StateVarStructureIV.Yhat)';
               MFEYMode = (1/T)*sum(StateVarStructureMode.Y-StateVarStructureMode.Yhat,2);
               MFEYIV = (1/T)*sum(StateVarStructureIV.Y-StateVarStructureIV.Yhat,2);
               %
               % get the mode values of any additional parameters
               %
               [thetaCalibratedMode,thetaAdditionalMode] = ReadAdditionalParameterValues(ModelParameters,ParameterNames);
               %
               % Save to disk
               %
               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-' DSGEModel.NameOfModel '.mat'];
                  save(ModeFile,'phiMode','thetaMode','phi','theta','phiEstimates','thetaCalibrated','thetaAdditionalInitial','thetaAdditionalMode','ParameterNames','ExtremeCorrParameters','LogPostDensity','LogLikeDensity',...
                       'LogPriorDensity','LogJacobianValue','InverseHessian','ThetaPhiPartial','thetaIndex','LowerBound','thetaPositions','thetaDist','PriorDist','ModelParameters','UniformBounds',...
                       'DSGEModel','AIMData','phiInverseHessianFinDiff','phiPostModeStdFinDiff','thetaInverseHessianFinDiff','thetaPostModeStdFinDiff',...
                       'phiFittedStd','thetaFittedStd','KalmanFirstObservation','FirstPeriod','LastPeriod');
               end;
               %
               % write output to file
               %
               SetLogText(strvcat(' ','Writing posterior mode estimation results to text file:',['   Directory: ' DSGEModel.OutputDirectory],['   File name: PosteriorModeResults-' DSGEModel.NameOfModel '.txt']));
               PrintPosteriorModeResults(DSGEModel,CurrINI,ParameterNames,phiMode,thetaMode,phi,theta,thetaCalibrated,thetaAdditionalInitial,thetaAdditionalMode,thetaDist,PriorDist,ThetaPhiPartial,InverseHessian,phiFittedStd,thetaFittedStd,phiPostModeStdFinDiff,thetaPostModeStdFinDiff,LogPostDensity,LogLikeDensity,LogPriorDensity,LogJacobianValue,ExtremeCorrParameters,OptAlgorithm,MaxTime,ToleranceValue,MaxIter,NumIter,FuncDiff,MaxEigenvalueIV,MaxEigenvalueMode,OptResults,LogPosteriorPhiIV,LogPosteriorThetaIV,LogPriorDensityIV,rank(R),MFEYMode,MSEYMode,MFEYIV,MSEYIV);
               %
               % write summary information to file
               %
               SetLogText(strvcat(' ','Writing posterior mode estimation summary results to text file:',['   Directory: ' DSGEModel.OutputDirectory],['   File name: PosteriorModeSummary-' DSGEModel.NameOfModel '.txt']));
               PrintPosteriorModeSummary(DSGEModel,CurrINI,ParameterNames,thetaMode,theta,thetaDist,PriorDist,ThetaPhiPartial,InverseHessian,phiFittedStd,LogPostDensity,LogLikeDensity,LogPriorDensity,LogJacobianValue,ExtremeCorrParameters,OptAlgorithm,MaxTime,ToleranceValue,MaxIter,NumIter,FuncDiff,MaxEigenvalueIV,MaxEigenvalueMode,LogPosteriorPhiIV,LogPosteriorThetaIV,LogPriorDensityIV,rank(R),MFEYMode,MSEYMode,MFEYIV,MSEYIV);
               SetLogText('All done!');
               %
               % open one of the results files
               %
               txt = 'Would you like to view only a summary of the posterior mode estimation output? [No means all printed results]';
               answer = Query(txt,'question',120,'Output Selection',500,'yes',CurrINI);
               drawnow;
               if strcmp(lower(answer),'yes')==1;
                  ModeTextFile = [DSGEModel.OutputDirectory '\PosteriorModeSummary-' DSGEModel.NameOfModel '.txt'];
               else;
                  ModeTextFile = [DSGEModel.OutputDirectory '\PosteriorModeResults-' DSGEModel.NameOfModel '.txt'];
               end;
               TextGUI(ModeTextFile,['Posterior Mode Estimation - ' DSGEModel.NameOfModel],[(CurrINI.scrsz(3)-1000)/2 32 1000 CurrINI.scrsz(4)-100],'Results',50,CurrINI,1,1);
            elseif (EstimationStatus==0)&(strcmp(lower(AbortEstimation),'cancel')==0)&(NumIter<=MaxIter);
               %
               % write error information to file
               %
               SetLogText(strvcat('Writing optimization error summary to disk.',['   Directory: ' DSGEModel.OutputDirectory],['   File name: OptErrorReport-' DSGEModel.NameOfModel '.txt']));
               %
               % transform phiMode to thetaMode taking into account that certain entries may be Inf or Nan's
               %
               NotFinite = 1-isfinite(phiMode);
               if sum(NotFinite)>=1;
                  HelpIndex = (1:NumParam)';
                  SwapParams = (NotFinite .* HelpIndex);
                  SwapParams = SwapParams(SwapParams>0);
                  phiModeOld = phiMode;
                  phiMode(SwapParams) = 1;
               end;
               thetaMode = PhiToTheta(phiMode,thetaIndex,UniformBounds,LowerBound);
               %
               % compute the partial of theta
               %
               phiGradientVector = abs(inv(PartialThetaPartialPhi(phiMode,thetaIndex,UniformBounds))*GradientVector);
               %
               % Now we have to reset all those elelemts that were either NaN or Inf
               %
               if sum(NotFinite)>=1;
                  thetaMode(SwapParams) = phiModeOld(SwapParams);
                  phiGradientVector(SwapParams) = abs(phiModeOld(SwapParams));
               end;
               PrintOptimizationErrorReport(DSGEModel,CurrINI,AIMData,ParameterNames,PriorDist,thetaDist,thetaMode,theta,phiGradientVector,GradientVector,rank(R),OptResults,ReturnCodeText);
               %
               % display the file
               %
               TxtFileName = [DSGEModel.OutputDirectory '\OptErrorReport-' DSGEModel.NameOfModel '.txt'];
               TextGUI(TxtFileName,['Optimization Error Summary - ' DSGEModel.NameOfModel],[(CurrINI.scrsz(3)-1000)/2 32 1000 CurrINI.scrsz(4)-100],'Results',50,CurrINI,1,1);
            else;
               if FileExist([pwd '\tmp\AiMErrorData.mat'])==1;
                  dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
               end;
            end;
         else;
            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. You need to change your initial parameter values.'];
            About(txt,'information','YADA - DSGE Model Problem',200,500,CurrINI);
            drawnow;
         end;
      else;
         %
         % 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 '" You have to change your initial parameter values.'];
         About(txt,'information','YADA - DSGE Model Problem',200,500,CurrINI);
         drawnow;
      end;
   else;
      SetLogText(strvcat(['   STATUS: Bad input file ' GetName(ErrorFile)],' '));
      About(ErrorStr,'error',['Error - ' GetFilename(ErrorFile)],200,500,CurrINI);
      drawnow;
   end;
elseif (status==1)&(isempty(theta)==1);
   SetLogText(strvcat('   STATUS: No parameters to estimate!',' '));
   txt = 'All parameters of your model have status calibrated. Hence, YADA doesn''t have any parameters to estimate.';
   About(txt,'information','Calibrated DSGE Model',200,500,CurrINI);
   drawnow;
else;
   SetLogText(strvcat('   STATUS: Bad prior file or inadequate data.',' '));
end;

%
% end of PosteriorModeEstimation.m
%
