function [logLikeValue,mcode,status] = logLikelihoodDSGE(ModelParameters,DSGEModel,AIMData)
% logLikelihoodDSGE: Calculates the value of the log-likelihood function at the current theta value for the DSGE model.
%
% USAGE:
%
%       [logLikeValue,mcode,status] = logLikelihoodDSGE(ModelParameters,DSGEModel,AIMData)
%
% REQUIRED INPUT:  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.
%
% REQUIRED OUTPUT: logLikeValue (scalar) the value of the log-likelihood function at the current theta value.
%
%                  mcode (integer) Indicates whether the solution of the model is unique or not.
%                                  mcode is 1 when unique, 0 if not computed, and various other values when
%                                  there are problems with the solution. Specifically:
%
%                                     mcode = 2:      Roots not correctly computed by real_schur
%                                     mcode = 3:      Too many big roots
%                                     mcode = 35:     Too many big roots, and q(:,right) is singular
%                                     mcode = 4:      Too few big roots
%                                     mcode = 45:     Too few big roots, and q(:,right) is singular
%                                     mcode = 5:      q(:,right) is singular
%                                     mcode = 61:     Too many exact shiftrights
%                                     mcode = 62:     Too many numeric shiftrights
%                                     mcode = 7:      The a matrix has infinite or NaN values. Eigenvalues
%                                                     cannot be calculated.
%                                     mcode = 8:      The function "compute_aim_matrices" returns complex
%                                                     numbers.
%                                     else            Return code not properly specified
%
%                  status (integer) Indicates if the log-likelihood was calculated (status=0), or not (status=1).
%
%
%                       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 8, 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:
%
% * 09-08-2006: Included the doubling algorithm options when calling the KalmanFilter function.
%
% * 14-08-2006: Conducted an output check. No changes made!
%
% * 15-08-2006: Added the sub-sample setting for the Y and X matrices when calling the Kalman filter.
%
% * 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.
%
% * 22-08-2006: Added the optional variable DSGEModel.KalmanFirstObservation for the Kalman filter. This variable
%               was previously always 1. It can now be any integer between 1 and ceil(0.5*T), where T is the sample size.
%
% * 25-08-2006: Removed the updating of additional parameters. This is handled by the function that evaluates the
%               log posterior instead.
%
% * 14-09-2006: Added the changes to the measurement equation call.
%
% * 15-09-2006: Made some input and output check. No changes to code.
%
% * 31-10-2006: Added the input "DAToleranceValue" for the call to the KalmanFilter function.
%
% * 17-11-2006: The function now allows for 3 choices of initializing the state covariance matrix for the
%               Kalman filter. The constant term for the alternative of a copnstant times the identity is
%               determined by the field "StateCovConstValue".
%
% * 20-11-2006: The initial value for the state vector is now user defined through the variable DSGEModel.InitialStateVector.
%
% * 12-11-2007: Updated the documentation.
%
% * 03-12-2007: Made sure the Kalman filter function includes the correct X matrix.
%
% * 21-12-2007: Added the code for dealing with unit-root state variables when running the Kalman filter.
%
% * 16-05-2008: The function now writes error data to the file "AiMErrorData.mat" when mcode is not unity.
%               Took the "DSGEModel.AIMTolerance" input for AiMSolver into account.
%
% * 23-05-2008: Updated the documentation.
%
% * 27-05-2008: The file "AiMErrorData.mat" is deleted if it exists and mcode is 1. The mcode value 7 is also
%               supported.
%
% * 02-07-2008: Added the variable "DSGEModel.AllowUnitRoot" for the Kalman filters.
%
% * 03-07-2008: Made sure that the measurement equation receives only the active state variables names.
%
% * 07-07-2008: Added code to allow for the case when the measurement matrix H is time-varying.
%
% * 21-01-2009: Updated the documentation.
%
% * 30-01-2009: Made sure that all outputs are initialized.
%
% * 09-02-2009: Changed the input "length(DSGEModel.StateVariablePositions)" to "size(F,1)" when calling the
%               measurement equation function.
%
% * 13-02-2009: Updated the documentation for mcode value 8.
%
% * 09-03-2009: Added quotation marks around paths for the dos del command.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

logLikeValue = 0;
mcode = 0;
status = 0;
%
% See if we can solve the model
%
[mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
if mcode~=1;
   [mcode,MaxAbsError,BMatrix,SMatrix,ModelEigenvalues] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
   save([pwd '\tmp\AiMErrorData.mat'],'mcode','ModelEigenvalues');
   return;
else;
   if FileExist([pwd '\tmp\AiMErrorData.mat'])==1;
      stat = dos(['del /Q "' pwd '\tmp\AiMErrorData.mat"']);
   end;
end;
%
% transform into state space
%
[F,B0] = AiMtoStateSpace(BMatrix,SMatrix,DSGEModel.StateVariablePositions,DSGEModel.StateShockPositions);
%
% Evaluate the measurement equation file
%
eval(['[A,H,R] = ' GetName(DSGEModel.MeasurementEquationFile) '(ModelParameters,AIMData.endog(DSGEModel.StateVariablePositions,:),DSGEModel.VariableNames,DSGEModel.XVariableNames,DSGEModel.n,size(F,1),DSGEModel.k);']);
%
% run the Kalman filter
%
if length(size(H))==2;
   if isempty(DSGEModel.UnitRootStates)==1;
      if isempty(DSGEModel.X)==0;
         [logLikeValue,status] = KalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.KalmanFirstObservation,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
      else;
         [logLikeValue,status] = KalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.KalmanFirstObservation,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
      end;
   else;
      if isempty(DSGEModel.X)==0;
         [logLikeValue,status] = UnitRootKalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.KalmanFirstObservation,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
      else;
         [logLikeValue,status] = UnitRootKalmanFilter(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.KalmanFirstObservation,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
      end;
   end;
else;
   if isempty(DSGEModel.UnitRootStates)==1;
      if isempty(DSGEModel.X)==0;
         [logLikeValue,status] = KalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.KalmanFirstObservation,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
      else;
         [logLikeValue,status] = KalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.KalmanFirstObservation,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
      end;
   else;
      if isempty(DSGEModel.X)==0;
         [logLikeValue,status] = UnitRootKalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),DSGEModel.X(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),A,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.KalmanFirstObservation,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
      else;
         [logLikeValue,status] = UnitRootKalmanFilterHt(DSGEModel.Y(:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),[],A,H(:,:,DSGEModel.FirstPeriod:DSGEModel.LastPeriod),F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.KalmanFirstObservation,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue,GetPosition(DSGEModel.StateVariablePositions,DSGEModel.UnitRootStates));
      end;
   end;
end;

%
% end of logLikelihoodDSGE.m
%
