function LogMargs = MargLikeChib(LogPiDensity,LogPsiDensity,LogPosterior,LogOmegaDensity,q,NumBurnin,ComputeSequential,SequentialStartIterationValue,SequentialStepLengthValue,WaitHandle)
% MargLikeChib: Computes the log marginal likelihood using Chibs marginal likelihood identity.
%
% USAGE:
%
%       LogMargs = MargLikeChib(LogPiDensity,LogPsiDensity,LogPosterior,LogOmegaDensity,q,NumBurnin,...
%                               ComputeSequential,SequentialStartIterationValue,SequentialStepLengthValue,WaitHandle)
%
% REQUIRED INPUT:  LogPiDensity (NumIter) vector with values of log p(PiMode|Psi(i),Omega(i),Data), i=1,...,NumIter
%
%                  LogPsiDensity (NumIter) vector with values of log p(PsiMode|PiMode,Omega(j),Data), j=1,...,NumIter
%                                provided that q>0 where q is the number of exogenous variables.
%
%                  LogPosterior (scalar) with the values of the sum of the log likelihood and the log prior at the
%                                posterior mode, i.e., log p(Data|PsiMode,PiMode,OmegaMode) + log p(PsiMode,PiMode,OmegaMode)
%
%                  LogOmegaDensity (scalar) with the value of log p(OmegaMode|PsiMode,PiMode,Data)
%
%                  q (interger) with the number of exogenous variables. If q=0, then PsiMode is an empty matrix.
%
%                  NumBurnin (integer) with the number of burn-in draws.
%
%                  ComputeSequential (boolean), which determines if the marginal likelihood should be calculated
%                               sequentially or just for the full sample.
%
%                  SequentialStartIterationValue (integer) determines the sample size where the sequential estimation
%                               should begin.
%
%                  SequentialStepLengthValue (integer) that gives the sample size increment for the sequential estimation
%                               of the log marginal likelihood.
%
% OPTIONAL INPUT:  WaitHandle (handle) to a wait dialog. Allows the user to end the function prematurely.
%
% REQUIRED OUTPUT: LogMargs (matrix) where the first column gives the sample size, the second column the value of
%                               the log marginal likelihood, and the third column the numerical standard error of
%                               the estimate.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2008 European Central Bank.
%
%                       First version: February 1, 2007.
%                        This version: December 19, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:
%
% * 2-2-2007:   Fixed a bug in the information sent to the NeweyWestCovMat function.
%
% * 8-2-2007:   Fixed a bug when using the cancel button.
%
% * 12-4-2007:  Changed the BarN value for the Newey-West function to floor(N^(1/2.01)).
%
% * 13-11-2007: Updated the documentation.
%
% * 27-3-2008:  Made sure that the sequential grid has at least 2 values.
%
% * 23-5-2008:  Updated the documentation.
%
% * 19-12-2008: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if (exist('WaitHandle')==0)|(nargin==9);
   WaitHandle = [];
end;
NumIter = length(LogPiDensity)-NumBurnin;
%
if ComputeSequential==1;
   Origin = SequentialStartIterationValue;
   %
   % make sure that Origin is less than NumIter
   %
   if Origin==NumIter;
      Origin = floor(NumIter/10);
   end;
   StepLength = SequentialStepLengthValue;
else;
   Origin = NumIter;
   StepLength = 1;
end;
%
% Create the sequential grid over which we will compute the marginal likelihood
%
SeqGrid = (Origin:StepLength:NumIter);
if isempty(SeqGrid)==1;
   SeqGrid = NumIter;
elseif SeqGrid(length(SeqGrid))<NumIter;
   SeqGrid = [SeqGrid NumIter];
end;
LogMargs = [SeqGrid' zeros(length(SeqGrid),2)];
%
% compute the max values
%
MaxLogPiDensity = max(LogPiDensity);
if q>0;
   MaxLogPsiDensity = max(LogPsiDensity);
end;
N = SeqGrid(1);
TempLogPiDensity = 0;
TempLogPsiDensity = 0;
LastN = 0;
NumCount = 1;
StopEstimation = 0;
for N=SeqGrid;
   if isempty(WaitHandle)==0;
      %
      % check if the cancel button was pressed
      %
      if strcmp(get(WaitHandle,'UserData'),'cancel')==1;
         StopEstimation = 1;
         break;
      end;
   end;
   %
   % First we compute the sequential value of the log density of Pi
   %
   TempLogPiDensity = TempLogPiDensity+sum(exp(LogPiDensity(1+LastN+NumBurnin:N+NumBurnin,1)-MaxLogPiDensity));
   LogPiDensityEst = MaxLogPiDensity+log(TempLogPiDensity/N);
   if q>0;
      TempLogPsiDensity = TempLogPsiDensity+sum(exp(LogPsiDensity(1+LastN+NumBurnin:N+NumBurnin,1)-MaxLogPsiDensity));
      LogPsiDensityEst = MaxLogPsiDensity+log(TempLogPsiDensity/N);
      %
      StdErrMat = NeweyWestCovMat([exp(LogPiDensity(1+NumBurnin:N+NumBurnin,1)) exp(LogPsiDensity(1+NumBurnin:N+NumBurnin,1))],floor(N^(1/2.01)));
      lnhPartial = [1/exp(LogPiDensityEst) 1/exp(LogPsiDensityEst)];
   else;
      LogPsiDensityEst = 0;
      %
      StdErrMat = NeweyWestCovMat(exp(LogPiDensity(1+NumBurnin:N+NumBurnin,1)),floor(SeqGrid(i)^(1/2.01)));
      lnhPartial = 1/exp(LogPiDensityEst);
   end;
   LogMargs(NumCount,2) = LogPosterior-(LogOmegaDensity+LogPsiDensityEst+LogPiDensityEst);
   if isempty(StdErrMat)==0;
      LogMargs(NumCount,3) = sqrt(lnhPartial*StdErrMat*lnhPartial');
   end;
   LastN = N;
   NumCount = NumCount+1;
   drawnow;
end;

if StopEstimation==1;
   LogMargs = [];
end;

%
% end of MargLikeChib.m
%
