function PrintOptimizationErrorReport(DSGEModel,CurrINI,AIMData,ParameterNames,PriorDist,thetaDist,thetaMode,theta,phiGradient,thetaGradient,NumMeasureErrors,OptResults,ReturnCodeText)
% PrintOptimizationErrorReport: Prints selective information about failed optimizations to file.
%
% USAGE:
%
%       PrintOptimizationErrorReport(DSGEModel,CurrINI,AIMData,ParameterNames,PriorDist,thetaDist,thetaMode,theta,phiGradient,thetaGradient,NumMeasureErrors,OptResults,ReturnCodeText)
%
% 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.
%
%                 AIMData (structure) with fields given by neq, nlag, nlead etc that is produced when running "compute_aim_data".
%                           The neq, nlag, and nlead fields are needed by AiMSolver.
%
%                 ParameterNames (structure) with fields "calibrated", "beta", "gamma", "normal", "invgamma", "truncnormal",
%                           "uniform", "all" and "estimated". Each field returns a string matrix with the parameter names.
%
%                 PriorDist (structure) with fields "beta", "beta_ab", "gamma", "gamma_ab", "normal", "invgamma",
%                           "truncnormal", and "uniform" where the values are matrices with the parameter values for the
%                           distributions. In the case of the normal, beta, and gamma these parameter values are
%                           the mean and the standard deviation. For the invgamma it is the s and q parameters,
%                           where q is a positive integer (degrees of freedom). For the left truncated normal
%                           there's a third parameter, the lower bound. Finally, for the uniform the parameters
%                           are the lower and the upper bound.
%
%                 thetaDist (vector) with integers 1 if gamma, 2 if beta, 3 if normal, 4 if invgamma, 5 if
%                           truncnormal, and 6 if uniform.
%
%                 thetaMode (vector) with the last reported parameter values.
%
%                 theta (vector) with the initial parameter values.
%
%                 phiGradient (vector) with the absoulte value of the gradient for the log posterior with respect
%                            to the transformed (phi) parameters.
%
%                 thetaGradient (vector) with the absolute value of the gradient for the log posterior with respect
%                            to the original (theta) parameters.
%
%                 NumMeasureErrors (integer), giving the number of measurement errors in the model.
%
%                 OptResults (matrix), of dimension (NumIter x 4) containing iteration no, value of log posterior, value
%                            of convergence measure, and the largest absolute value of the gradient.
%
%                 ReturnCodeText (string vector) with a message from the optimization routine.
%
%
%                       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 22, 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:
%
% * 13-11-2007: Updated the documentation.
%
% * 2-1-2008:   Removed code using "DSGEModel.InitializeKalman".
%
% * 16-5-2008:  The function now also tries to check if an error report was written to disk after running AiM.
%
% * 19-5-2008:  Fixed a spelling bug for the "AIMData" input variable.
%
% * 20-5-2008:  Added the two input variables "phiGradient" and "thetaGradient".
%
% * 23-5-2008:  Updated the documentation.
%
% * 26-5-2008:  Added a newline command before displaying the eigenvalues.
%
% * 27-5-2008:  The mcode value 7 is now taken into account.
%
% * 3-7-2008:   Added the input variable "ReturnCodeText".
%
% * 22-7-2008:  Took into account that the gamma and inverse gamma distributions now use
%               the lower bound parameter.
%
% * 23-7-2008:  Took into account that the beta distribution supports general lower and upper bounds.
%
% * 21-01-2009: Updated the documentation.
%
% * 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 command.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

TxtFileName = [DSGEModel.OutputDirectory '\OptErrorReport-' DSGEModel.NameOfModel '.txt'];
fid = fopen(TxtFileName,'wt');
%
fprintf(fid,'*********************************************************************************************************\n');
fprintf(fid,'*                                                                                                       *\n');
fprintf(fid,'*                  P O S T E R I O R   M O D E   E S T I M A T I O N   W I T H   Y A D A                *\n');
fprintf(fid,'*                                         E R R O R   R E P O R T                                       *\n');
fprintf(fid,'*                                                                                                       *\n');
fprintf(fid,'*********************************************************************************************************\n\n');
%
% Model info
%
fprintf(fid,'Name of model:                    %s\n',DSGEModel.NameOfModel);
fprintf(fid,'Data construction file:           %s\n',GetFilename(DSGEModel.DataConstructionFile));
fprintf(fid,'Measurement equation file:        %s\n',GetFilename(DSGEModel.MeasurementEquationFile));
fprintf(fid,'Prior distribution data file:     %s\n',GetFilename(DSGEModel.PriorFile));
fprintf(fid,'Sheet name:                       %s\n',DSGEModel.PriorFileSheet);
if FileExist(DSGEModel.UpdateParameterFile)==1;
   fprintf(fid,'Update parameters file:           %s\n',GetFilename(DSGEModel.UpdateParameterFile));
end;
if FileExist(DSGEModel.InitializeParameterFile)==1;
   fprintf(fid,'Initialize parameters file:       %s\n',GetFilename(DSGEModel.InitializeParameterFile));
end;
fprintf(fid,'AiM model file:                   %s\n\n',GetFilename(DSGEModel.AIMFile));
fprintf(fid,'*********************************************************************************************************\n');
fprintf(fid,'*                                                                                                       *\n');
fprintf(fid,'*                           S A M P L E   &   D A T A   I N F O R M A T I O N                           *\n');
fprintf(fid,'*                                                                                                       *\n');
fprintf(fid,'*********************************************************************************************************\n\n');
%
% Sample information
%
[NewStartYear,NewStartPeriod] = AdjustSampleStart(DSGEModel.SubBeginYear,DSGEModel.SubBeginPeriod,DSGEModel.DataFrequency,DSGEModel.KalmanFirstObservation-1);
%
fprintf(fid,'Selected sample:                  %s:%s -- %s:%s\n',DSGEModel.SubBeginYear,DSGEModel.SubBeginPeriod,DSGEModel.SubEndYear,DSGEModel.SubEndPeriod);
fprintf(fid,'Number of observations:           %0.0f\n',DSGEModel.LastPeriod-DSGEModel.FirstPeriod+1);
fprintf(fid,'Number of observations used to\n');
fprintf(fid,'initialize the Kalman filter:     %0.0f\n',DSGEModel.KalmanFirstObservation-1);
fprintf(fid,'Effective number of observations: %0.0f\n',DSGEModel.LastPeriod-DSGEModel.FirstPeriod+1-(DSGEModel.KalmanFirstObservation-1));
fprintf(fid,'Estimation sample:                %s:%s -- %s:%s\n',NewStartYear,NewStartPeriod,DSGEModel.SubEndYear,DSGEModel.SubEndPeriod);
%
fprintf(fid,'Number of observed variables:     %0.0f\n',min(size(DSGEModel.Y)));
fprintf(fid,'Number of exogenous variables:    %0.0f\n',min(size(DSGEModel.X)));
fprintf(fid,'Number of state variables:        %0.0f\n',size(DSGEModel.StateVariableNames,1));
fprintf(fid,'Number of state shocks:           %0.0f\n',size(DSGEModel.StateShockNames,1));
fprintf(fid,'Number of measurement errors:     %0.0f\n\n',NumMeasureErrors);
%
% observed variable names, state variable names and state shocks
%
VarNames = StrWrap(vecstr(DSGEModel.VariableNames,' '),60);
fprintf(fid,'Names of observed variables:      %s\n',VarNames(1,:));
for i=2:size(VarNames,1);
   if strcmp(StringTrim(VarNames(i,:)),'')==0;
      fprintf(fid,'                                  %s\n',VarNames(i,:));
   end;
end;
if min(size(DSGEModel.X))>0;
   XVarNames = StrWrap(vecstr(DSGEModel.XVariableNames,' '),60);
   fprintf(fid,'Names of exogenous variables:     %s\n',XVarNames(1,:));
   for i=2:size(XVarNames,1);
      if strcmp(StringTrim(XVarNames(i,:)),'')==0;
         fprintf(fid,'                                  %s\n',XVarNames(i,:));
      end;
   end;
end;
StateVarNames = StrWrap(vecstr(DSGEModel.StateVariableNames,' '),60);
fprintf(fid,'Names of state variables:         %s\n',StateVarNames(1,:));
for i=2:size(StateVarNames,1);
   if strcmp(StringTrim(StateVarNames(i,:)),'')==0;
      fprintf(fid,'                                  %s\n',StateVarNames(i,:));
   end;
end;
StateShockNames = StrWrap(vecstr(DSGEModel.StateShockNames,' '),60);
fprintf(fid,'Names of the state shocks:        %s\n',StateShockNames(1,:));
for i=2:size(StateShockNames,1);
   if strcmp(StringTrim(StateShockNames(i,:)),'')==0;
      fprintf(fid,'                                  %s\n',StateShockNames(i,:));
   end;
end;
%
%
prt_val = ['%' num2str(6+CurrINI.decimals,'%0.0f') '.' num2str(CurrINI.decimals,'%0.0f') 'f'];
prt_str = ['|    %10.0f          ' prt_val '           ' prt_val '         ' prt_val '     |\n'];
fprintf(fid,'\n---------------------------------------------------------------------------------------\n');
fprintf(fid,'| Iteration Number   log posterior of phi    convergence measure   max(abs(gradient)) |\n');
fprintf(fid,'|-------------------------------------------------------------------------------------|\n');
fprintf(fid,prt_str,OptResults');
fprintf(fid,'---------------------------------------------------------------------------------------\n\n');
%
% check if we have information from the optimizer
%
if isempty(ReturnCodeText)==0;
   fprintf(fid,'********************************************************************************\n');
   fprintf(fid,'*                                                                              *\n');
   fprintf(fid,'*       O P T I M I Z A T I O N   R O U T I N E   R E T U R N   C O D E        *\n');
   fprintf(fid,'*                                                                              *\n');
   fprintf(fid,'********************************************************************************\n\n');
   fprintf(fid,'%s\n\n',ReturnCodeText);
end;
%
% Parameter information
%
fprintf(fid,'********************************************************************************\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'*                             P A R A M E T E R S                              *\n');
fprintf(fid,'*                                                                              *\n');
fprintf(fid,'********************************************************************************\n\n');
%
%
% check the length of parameter names string
%
ns = size(ParameterNames.estimated,2);
if ns<20;
   add_str = SpaceStr(20-ns);
else;
   add_str = '';
end;
prt_str_1 = ['%s%s    ' prt_val '    ' prt_val '    ' prt_val '    ' prt_val];
prt_str_2 = ['        %s(%s    ' prt_val '    ' prt_val];
%
% Calculate an approximation of the inverse Hessian for theta, the original parameters, using the
% inverse Hessian at the mode for phi, the transformed parameters, and the partial derivatives of
% theta wrt phi
%
%
fprintf(fid,'======================================================================================\n');
fprintf(fid,'                                 Original Parameters                                  \n');
fprintf(fid,'      parameter          theta last        theta init   phi gradient  theta gradient\n');
if (isempty(PriorDist.truncnormal)==1)&(isempty(PriorDist.invgamma)==1)&(isempty(PriorDist.gamma)==1)&(isempty(PriorDist.beta)==1);
   fprintf(fid,'     (prior dist        parameter 1     parameter 2)\n');
else;
   if isempty(PriorDist.beta)==1;
      fprintf(fid,'     (prior dist        parameter 1     parameter 2     parameter 3)\n');
   else;
      fprintf(fid,'     (prior dist        parameter 1     parameter 2     parameter 3\n');
      fprintf(fid,'                                                        parameter 4)\n');
   end;
end;
%
n_gamma = 0;
n_beta = 0;
n_normal = 0;
n_invgamma = 0;
n_truncnormal = 0;
n_uniform = 0;
for i=1:length(thetaMode);
   if thetaDist(i)==1;
      dist_str = 'gamma';
      dist_add = SpaceStr(6);
      n_gamma = n_gamma+1;
      prior_val = PriorDist.gamma(n_gamma,:);
   elseif thetaDist(i)==2;
      dist_str = 'beta';
      dist_add = SpaceStr(7);
      n_beta = n_beta+1;
      prior_val = PriorDist.beta(n_beta,:);
   elseif thetaDist(i)==3;
      dist_str = 'normal';
      dist_add = SpaceStr(5);
      n_normal = n_normal+1;
      prior_val = PriorDist.normal(n_normal,:);
   elseif thetaDist(i)==4;
      dist_str = 'invgamma';
      dist_add = SpaceStr(3);
      n_invgamma = n_invgamma+1;
      prior_val = PriorDist.invgamma(n_invgamma,:);
   elseif thetaDist(i)==5;
      dist_str = 'truncnormal';
      dist_add = '';
      n_truncnormal = n_truncnormal+1;
      prior_val = PriorDist.truncnormal(n_truncnormal,:);
   elseif thetaDist(i)==6;
      dist_str = 'uniform';
      dist_add = SpaceStr(4);
      n_uniform = n_uniform+1;
      prior_val = PriorDist.uniform(n_uniform,:);
   end;
   fprintf(fid,'--------------------------------------------------------------------------------------\n');
   fprintf(fid,[prt_str_1 '\n'],add_str,ParameterNames.estimated(i,:),thetaMode(i),theta(i),phiGradient(i),thetaGradient(i));
   if (strcmp(dist_str,'truncnormal')==0)&(strcmp(dist_str,'invgamma')==0)&(strcmp(dist_str,'gamma')==0)&(strcmp(dist_str,'beta')==0);
      fprintf(fid,[prt_str_2 ')\n'],dist_add,dist_str,prior_val(1,1),prior_val(1,2));
   else;
      if strcmp(dist_str,'beta')==0;
         fprintf(fid,[prt_str_2 '    ' prt_val ')\n'],dist_add,dist_str,prior_val(1,1),prior_val(1,2),prior_val(1,3));
      else;
         fprintf(fid,[prt_str_2 '    ' prt_val '\n'],dist_add,dist_str,prior_val(1,1),prior_val(1,2),prior_val(1,3));
         fprintf(fid,['                                                        ' prt_val ')\n'],prior_val(1,4));
      end;
   end;
end;
fprintf(fid,'======================================================================================\n\n');
%
% check if error data from AiM was stored in the tmp directory
%
AiMErrorFile = [pwd '\tmp\AiMErrorData.mat'];
if FileExist(AiMErrorFile)==1;
   AiMErrorData = load(AiMErrorFile);
   %
   if isfield(AiMErrorData,'mcode')==1;
      mcode = AiMErrorData.mcode;
      ModelEigenvalues = AiMErrorData.ModelEigenvalues;
      %
      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;
      %
      % write 
      %
      fprintf(fid,'********************************************************************************\n');
      fprintf(fid,'*                                                                              *\n');
      fprintf(fid,'*                  A I M   E R R O R   I N F O R M A T I O N                   *\n');
      fprintf(fid,'*                                                                              *\n');
      fprintf(fid,'********************************************************************************\n\n');
      fprintf(fid,'The AiM solver provided the return code: %s\n\n',mcodeStr);
      fprintf(fid,'                                         at the last parameter vector.\n');
      %
      % write ModelEigenvalues data to file
      %
      fprintf(fid,'AiM Model File:                            %s\n',GetFilename(DSGEModel.AIMFile));
      fprintf(fid,'Number of equations:                       %0.0f\n',AIMData.neq);
      fprintf(fid,'Number of lags:                            %0.0f\n',AIMData.nlag);
      fprintf(fid,'Number of leads:                           %0.0f\n',AIMData.nlead);
      fprintf(fid,'Number of exact shiftrights:               %0.0f\n',ModelEigenvalues.NumExact);
      fprintf(fid,'Number of numeric shiftrights:             %0.0f\n',ModelEigenvalues.NumNumeric);
      fprintf(fid,'Number of large roots:                     %0.0f\n',ModelEigenvalues.LargeRoots);
      fprintf(fid,'Number of stability conditions:            %0.0f\n',ModelEigenvalues.NumExact+ModelEigenvalues.NumNumeric+ModelEigenvalues.LargeRoots);
      fprintf(fid,'Number of required stability conditions:   %0.0f\n',AIMData.neq*AIMData.nlead);
      fprintf(fid,'Difference:                                %0.0f\n',ModelEigenvalues.DimStability);
      fprintf(fid,'Dimension of DSGE state transition matrix: %0.0f\n',ModelEigenvalues.DimCompanion);
      %
      if mcode~=7;
         fprintf(fid,'\n                           Eigenvalues\n');
         fprintf(fid,'Number       Real          Imaginary      Amplitude         Period\n');
         prt_str = [' %3.0f     ' prt_val '   ' prt_val '   ' prt_val '   ' prt_val ' \n'];
         for i=1:length(ModelEigenvalues.Roots);
            if i<=length(ModelEigenvalues.Amplitude);
               fprintf(fid,prt_str,i,real(ModelEigenvalues.Roots(i)),imag(ModelEigenvalues.Roots(i)),ModelEigenvalues.Amplitude(i),ModelEigenvalues.Period(i));
            else;
               fprintf(fid,prt_str,i,real(ModelEigenvalues.Roots(i)),imag(ModelEigenvalues.Roots(i)),0,NaN);
            end;
         end;
      end;
      %
      fprintf(fid,'\n********************************************************************************\n\n');
   end;
   %
   % delete the error data file.
   %
   dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
end;
%
% Some information about the filename and its location
%
fprintf(fid,'Directory for file: %s\n',GetPath(TxtFileName));
fprintf(fid,'Name of file:       %s\n',GetFilename(TxtFileName));
fprintf(fid,'Output created on:  %s\n\n',datestr(now,0));
fclose(fid);

%
% end of PrintOptimizationErrorReport.m
%
