function PlotMargLikeChibJeliazkov(PostSample,LogPost,NumBurnin,CurrChain,DSGEModel,CurrINI)
% PlotMargLikeChibJeliazkov: Calculates sequentially the Chib and Jeliazkov marginal likelihood
%                            identity based estimator of the log marginal likelihood.
%
% USAGE:
%
%       PlotMargLikeChibJeliazkov(PostSample,LogPost,NumBurnin,CurrChain,DSGEModel,CurrINI)
%
% REQUIRED INPUT: PostSample (matrix) of dimension NumIterxNumParam with the posterior sample draws
%                            of the phi parameters.
%
%                 LogPost (vector) of dimension (NumIterx1) with all the log posterior values for the
%                            non-discarded posterior draws. These posterior values are calculated
%                            using the "logPosteriorPhiDSGE" function.
%
%                 NumBurnin (integer), with the number of burn-in draws.
%
%                 CurrChain (integer) with the current MCMC chain number.
%
%                 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.
%
%
%                       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: December 13, 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:
%
% * 29-03-2007: Updated the documentation.
%
% * 10-04-2007: Took the name change RWMH to RWM into account.
%
% * 12-04-2007: Made sure that the MargLikeChibJeliazkov function received a unit value for SequentialML.
%
% * 17-04-2007: Removed a SetLogText command.
%
% * 09-11-2007: Updated the documentation.
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 11-03-2008: Added code for writing output to file.
%
% * 12-03-2008: Random number generators are now initialized properly.
%
% * 23-05-2008: Updated the documentation.
%
% * 28-05-2008: Changed name text on figure window.
%
% * 03-07-2008: Made sure that the measurement equation receives only the active state variable names.
%
% * 19-12-2008: 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 function called additional functions
%               so that the dimension of F could be retrieved.
%
% * 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 copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% locate the controls structure
%
maingui = findobj('Type','figure','Tag','YADA');
controls = get(maingui,'UserData');
%
% Number of posterior draws
%
DrawsStr = get(controls.posterior.draws,'String');
NumPostDraws = str2num(StringTrim(strrep(DrawsStr(DSGEModel.PosteriorDrawsValue,:),',','')));
%
% 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,:)));
%
% 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,:)),',',''));
%
% Try to 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)==0;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
phiMode = ModeData.phiMode;
SigmaMode = ModeData.InverseHessian;
%
% 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;
logPostMode = ModeData.LogPostDensity;
%
% sampling data
%
DSGEModel.FirstPeriod = ModeData.FirstPeriod;
DSGEModel.LastPeriod = ModeData.LastPeriod;
DSGEModel.KalmanFirstObservation = ModeData.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;
%
% 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;
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 DSGEModel.RunInitializeFirst==1;
      %
      % check if the file with additional parameters that only needs to be updated is OK
      % 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);
         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;
         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);
         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;
%
% check if the model can be solved
%
if isempty(ErrorStr)==1;
   try;
      [mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
   catch;
      ErrorStr = ['Error found in "compute_aim_matrices.m". Message caught is: ' lasterr];
      ErrorFile = [GetPath(DSGEModel.AIMDataFile) 'compute_aim_matrices.m'];
   end;
end;
%
if (isempty(ErrorStr)==1)&(mcode==1);
   try;
      [F,B0] = AiMtoStateSpace(BMatrix,SMatrix,DSGEModel.StateVariablePositions,DSGEModel.StateShockPositions);
   catch;
      ErrorStr = ['Error found when using the output from "compute_aim_matrices.m". Message caught is: ' lasterr];
      ErrorFile = [GetPath(DSGEModel.AIMDataFile) 'compute_aim_matrices.m'];
   end;
elseif (isempty(ErrorStr)==1)&(mcode~=1);
   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;
   ErrorStr = ['The DSGE model doesn''t have a unique covergent solution. AiM returns the status: ' mcodeStr];
   ErrorFile = [GetPath(DSGEModel.AIMDataFile) 'compute_aim_matrices.m'];
end;
%
% check the measurement equation file for internal validity
%
if isempty(ErrorStr)==1;
   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;
   About(ErrorStr,'error',['Error - ' GetFilename(ErrorFile)],120,500,CurrINI);
   return;
end;
%
% initialize the random number generators
%
%
if DSGEModel.RandomNumberValue==1;
   rand('state',0);
   randn('state',0);
else;
   rand('state',sum(100*clock));
   randn('state',sum(100*clock));
end;
%
% Open up a wait dialog.
%
txt = ['Please wait while YADA computes Chib and Jeliazkov''s marginal likelihood identity based estimator of the log marginal likelihood. Estimation started at: ' StringTrim(datestr(now,14)) '. This will take some time!'];
WaitHandle = WaitDLG(txt,'information','Sequential Marginal Likelihood',500,150,CurrINI,1);
WaitControls = get(WaitHandle,'UserData');
%
% run the Chib and Jeliazkov estimator
%
LogMargsCJ = MargLikeChibJeliazkov(PostSample,LogPost,phiMode,SigmaMode,c,logPostMode,NumBurnin,1,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData,WaitControls.cancel);
%
% shut down the wait dialog
%
set(WaitControls.text,'UserData','done');
delete(WaitHandle);
drawnow;
pause(0.02);
%
if isempty(LogMargsCJ)==0;
   %
   % ask if we write to file
   %
   txt = 'Would you like to save the sequential estimation data to file?';
   answer = Query(txt,'question',200,'Save - Chib and Jeliazkov',500,'yes',CurrINI);
   if strcmp(lower(answer),'yes')==1;
      Xvalues = LogMargsCJ(:,1)+NumBurnin;
      LogMargLikeCJ = LogMargsCJ(:,2:3);
      %
      file = [DSGEModel.OutputDirectory '\SeqChibJeliazkov-' DSGEModel.NameOfModel '-' int2str(NumPostDraws) '-' int2str(CurrChain) '.mat'];
      save(file,'Xvalues','LogMargLikeCJ');
      %
      mfile = [DSGEModel.OutputDirectory '\SeqCJ' strrep(DSGEModel.NameOfModel,'-','') int2str(NumPostDraws) int2str(CurrChain) '.m'];
      fid = fopen(mfile,'wt');
      fprintf(fid,['%%\n%% load the data in ' GetFilename(file) '\n%%\n']);
      fprintf(fid,'DataStructure = load(''%s'');\n',GetFilename(file));
      fprintf(fid,'%%\n%% Number of posterior draws used for sequential estimation (N x 1) vector\n%%\n');
      fprintf(fid,'Xvalues = DataStructure.Xvalues;\n');
      fprintf(fid,'%%\n%% N x 2 matrix with log marginal likelihood values and numerical standard\n%% error using the Chib and Jeliazkov estimator\n%%\n');
      fprintf(fid,'LogMargLikeCJ = DataStructure.LogMargLikeCJ;\n');
      fprintf(fid,'%%\n%% Add your own commands below\n%%\n\n\n');
      %
      % end the script file
      %
      fprintf(fid,'%%\n%% Created by YADA on %s\n%%\n',datestr(now,0));
      fclose(fid);
      %
      txt = ['The data on sequential estimation of the log marginal likelihood has been saved to the file "' GetFilename(file) '" in the directory "' GetPath(file) '". The file contains 2 entries: Xvalues (Nx1 vector with the number of posteror draws used for the sequential estimation), and LogMargLikeCJ (Nx2 matrix with estimates of the log marginal likelihood and the numerical standard errors with the Chib and Jeliazkov approach). To access this data you may run the matlab script file "' GetFilename(mfile) '".'];
      About(txt,'information','Chib and Jeliazkov',200,500,CurrINI);
   end;
   %
   GraphDirectory = [DSGEModel.OutputDirectory '\graphs'];
   %
   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'];
      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;

%
% end of PlotMargLikeChibJeliazkov.m
%
