function LogMargs = MargLikeChibJeliazkov(PostSample,LogPost,phiMode,SigmaMode,c,logPostMode,NumBurnin,ComputeSequential,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData,WaitHandle)
% MargLikeChibJeliazkov: Computes the marginal likelihood from a MCMC sample using
%                        the Chib and Jeliazkov estimator along with numerical standard
%                        errors based on Newey-West correction.
%
% USAGE:
%
%       LogMargs = MargLikeChibJeliazkov(PostSample,LogPost,phiMode,SigmaMode,c,logPostMode,NumBurnin,ComputeSequential,...
%                                        thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,...
%                                        ModelParameters,DSGEModel,AIMData,WaitHandle)
%
% REQUIRED INPUT:  PostSample (matrix) of dimension (NumIter x NumParam) with posterior draws of parameters.
%
%                  LogPost (vector) of dimension (NumIter x 1) with log posterior values for the parameters.
%
%                  phiMode (vector) of dimension NumParam with the posterior mode estimate of the parameters.
%
%                  SigmaMode (matrix) of dimension (NumParam x NumParam) with the inverse Hessian evaluated at
%                            the posterior mode.
%
%                  c (real) with the scale factor for the proposal density
%
%                  logPostMode (real) with the value of the log of the posterior of phi at the mode.
%
%                  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
%
%                  thetaIndex (vector) indicating the type of distribution the parameters in theta are drawn from.
%                            If 1, then the distribution is gamma, inverse gamma, or left truncated normal.
%                            If 2, then the distribution is standardized beta. If 0, then the distribution
%                            is normal.
%
%                  UniformBounds (matrix) with the same number of rows as parameters and 2 columns. The columns
%                            are 0 1 except for the uniform and the beta distribution, where they represent
%                            the upper and lower bounds.
%
%                  LowerBound (vector) with the lower bounds for the distributions.
%
%                  thetaPositions (structure) with fields parameter and distribution. This function
%                            needs the former field to know which parameter name corresponds to
%                            position i in theta.
%
%                  thetaDist (vector) with integers 1 if gamma, 2 if beta, 3 if normal, 4 if invgamma, and 5 if
%                            truncnormal.
%
%                  PriorDist (structure) with fields "beta", "gamma", "normal", "invgamma", and "truncnormal", 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). Finally, for the left truncated normal there's a third parameter, the lower bound.
%
%                  ModelParameters (structure) whose fields have names given by the parameter names in the
%                            DSGE model.
%
%                  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.
%
%                  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.
%
% OPTIONAL INPUT:  WaitHandle (handle) to a wait dialog. Allows the user to end the function prematurely.
%
% REQUIRED OUTPUT: LogMargs (matrix) with the calculated marginal likelihoods. The sample size
%                            for the sequential estimate is located in the first column, while
%                            the estimated marginal likelihood is given in the second column, and
%                            the numerical Newey-West standard errors in the third column.
%
%
%                       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: November 29, 2006.
%                        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:
%
% * 30-11-2006: Finished the coding of the function.
%
% * 13-12-2006: Added a drawnow command. The optional input "WaitHandle" was also added.
%
% * 14-12-2006: Updated the documentation.
%
% * 8-1-2007:   If "SeqGrid" is empty, it is now set to NumIter.
%
% * 12-4-2007:  Changed the BarN value for the Newey-West function to floor(N^(1/2.01)).
%
% * 13-11-2007: Updated the documentation.
%
% * 12-3-2008:  Added the code for updating the text string on the wait handle.
%
% * 27-3-2008:  Made sure that the sequential grid has at least 2 values.
%
% * 23-5-2008:  Updated the documentation.
%
% * 23-7-2008:  Updated the documentation.
%
% * 19-12-2008: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if (exist('WaitHandle')==0)|(nargin==17);
   WaitHandle = [];
   WaitTxt = '';
else;
   WaitControls = get(get(WaitHandle,'Parent'),'UserData');
   WaitTxt = get(WaitControls.text,'String');
   k = size(WaitTxt,1);
   kLine = WaitTxt{k};
end;
%
% initializing elements
%
LogMargs = [];
InvSigma = (c^(-2))*inv(SigmaMode);
LogDetSigma = log(det((c^2)*SigmaMode));
CholeskiSigma = chol(SigmaMode)';
NumIter = length(LogPost)-NumBurnin;
NumParm = length(phiMode);
%
if ComputeSequential==1;
   Origin = DSGEModel.SequentialStartIterationValue;
   %
   % make sure that Origin is less than NumIter
   %
   if Origin==NumIter;
      Origin = floor(NumIter/10);
   end;
   StepLength = DSGEModel.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;
%
% Evaluate all the accepted draws using the proposal density
%
NumIter = length(LogPost)-NumBurnin;
lnq = zeros(NumIter,1);
for i=1:NumIter;
   lnq(i,1) = logMultiNormalPDF(PostSample(i+NumBurnin,:)',phiMode,LogDetSigma,InvSigma);
   drawnow;
end;
%
% The numerator of eq. (9) in Chib and Jeliazkov is given by the mean of exp(lnq) since the alpha weights
% are always 1 when the log of the posterior at the mode is used in the numerator of the definition of alpha.
%
StopEstimation = 0;
AcceptFcn = zeros(NumIter,1);
for i=1:NumIter;
   if isempty(WaitHandle)==0;
      %
      % check if the cancel button was pressed
      %
      if strcmp(get(WaitHandle,'UserData'),'cancel')==1;
         StopEstimation = 1;
         break;
      end;
   end;
   %
   % check if we should update the wait dialog
   %
   if (isempty(WaitTxt)==0)&(sum(SeqGrid==i)>0);
      WaitTxt{k} = [kLine ' Acceptance probabilities for sample size: ' int2str(i) ' (' int2str(NumIter) ').'];
      set(WaitControls.text,'String',WaitTxt);
   end;
   %
   % draw from proposal density
   %
   phi = phiMode+(c*CholeskiSigma*randn(NumParm,1));
   %
   % evaluate the log posterior
   %
   logProp = -logPosteriorPhiDSGE(phi,thetaIndex,UniformBounds,LowerBound,thetaPositions,thetaDist,PriorDist,ModelParameters,DSGEModel,AIMData);
   AcceptFcn(i,1) = min(1,exp(logProp-logPostMode));
   drawnow;
end;
%
% check if the estimation was canceled
%
if StopEstimation==1;
   return;
end;
%
% now we can do the sequential estimation
%
LogMargs = zeros(length(SeqGrid),1);
StdErrorVec = zeros(length(SeqGrid),1);
for i=1:length(SeqGrid);
   %
   % check if we should update the wait dialog
   %
   if isempty(WaitTxt)==0;
      WaitTxt{k} = [kLine ' Log marginal likelihood for sample size: ' int2str(SeqGrid(i)) ' (' int2str(NumIter) ').'];
      set(WaitControls.text,'String',WaitTxt);
   end;
   h = [(1/SeqGrid(i))*sum(exp(lnq(1:SeqGrid(i),1)));(1/SeqGrid(i))*sum(AcceptFcn(1:SeqGrid(i),1))];
   LogMargs(i,1) = logPostMode-(log(h(1))-log(h(2)));
   %
   % compute numerical standard errors.
   %
   StdErrMat = NeweyWestCovMat([exp(lnq(1:SeqGrid(i),1)) AcceptFcn(1:SeqGrid(i),1)],floor(SeqGrid(i)^(1/2.01)));
   lnhPartial = [1/h(1) -1/h(2)];
   StdErrorVec(i,1) = sqrt(lnhPartial*StdErrMat*lnhPartial');
   drawnow;
end;
LogMargs = [SeqGrid' LogMargs StdErrorVec];

%
% end of MargLikeChibJeliazkov.m
%
