function InformationMatrix = DSGEInformationMatrix(theta,thetaPositions,ModelParameters,ParameterNames,DSGEModel,CurrINI)
% DSGEInformationMatrix: Computes Fisher's information matrix for the DSGE model.
%
% USAGE:
%
%       InformationMatrix = DSGEInformationMatrix(theta,thetaPositions,ModelParameters,ParameterNames,DSGEModel,CurrINI)
%
% REQUIRED INPUT:  theta (vector) with values of the original DSGE model parameters
%
%                  thetaPositions (structure) with fields parameter and distribution. This function
%                                 needs the former field to know which parameter name corresponds to
%                                 position i in theta.
%
%                  ModelParameters (structure) whose fields have names given by the parameter names in the
%                                 DSGE model.
%
%                  ParameterNames (string matrix) with the names of the parameters.
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure) with initialization information.
%
% REQUIRED OUTPUT: InformationMatrix (matrix) of the same dimension as theta.
%
%
%                       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: February 12, 2008.
%                        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:
%
% * 23-05-2008: Updated the documentation.
%
% * 29-05-2008: Added a test to check if the model can be solved at the mode.
%
% * 02-07-2008: Added the variable "DSGEModel.AllowUnitRoot" for the Kalman filters.
%
% * 07-07-2008: Added code to allow for the case when the measurement matrix H is time-varying.
%
% * 22-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 copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize the output
%
InformationMatrix = [];
%
% copy files to the tmp folder
%
[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;
%
% load the AIMData file
%
ErrorStr = '';
try;
   AIMData = load(DSGEModel.AIMDataFile);
catch;
   ErrorStr = ['Unable to load the file "' DSGEModel.AIMDataFile '". Message caught is: ' lasterr];
end;
if isempty(ErrorStr)==1;
   %
   % check if we can solve the model
   %
   ErrorStr = '';
   try;
      [A,H,R,F,B0,mcode] = SolveDSGEModel(DSGEModel,theta,thetaPositions,ModelParameters,AIMData);
   catch;
      ErrorStr = ['YADA caught an error when trying to solve the DSGE model. Message caught is:' lasterr];
   end;
   if isempty(ErrorStr)==0;
      About(ErrorStr,'information','YADA - Error Message',200,500,CurrINI);
      return;
   end;
   if 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;
      txt = ['The AiM solver provided the return code: ' num2str(mcode) ', i.e., "' mcodeStr '"'];
      About(txt,'information','Information Matrix',120,500,CurrINI);
      return;
   end;
else;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% get the sample
%
[DSGEModel.FirstPeriod,DSGEModel.LastPeriod] = CreateSubSample(DSGEModel);
%
if length(size(H))==3;
   if size(H,3)<DSGEModel.LastPeriod;
      %
      % too few data points in the 3D H matrix.
      %
      ErrorStr = ['The number of time periods for the time-varying measurement matrix H (' int2str(size(H,3)) ') is lower than the number of observations (T = ' int2str(DSGEModel.LastPeriod) '). YADA has therefore aborted from the information matrix computation.'];
      About(ErrorStr,'error','Error - Measurement Equation',200,500,CurrINI);
      return;
   end;
end;
%
% set the InitialStateVector value
%
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;
%
% run the Kalman filter
%
if length(size(H))==2;
   if isempty(DSGEModel.UnitRootStates)==1;
      if isempty(DSGEModel.X)==0;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilterHt(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);
      else;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilterHt(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);
      end;
   else;
      if isempty(DSGEModel.X)==0;
         [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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;
%
% get the default values of ytilde and B
%
Y = DSGEModel.Y(:,DSGEModel.FirstPeriod+DSGEModel.KalmanFirstObservation-1:DSGEModel.LastPeriod);
[n,T] = size(Y);
p = length(theta);
Ytilde = Y-Yhat;
B = MSEY;
vecB = zeros(n*n,T);
for t=1:T;
   vecB(:,t) = vec(MSEY(:,:,t));
end;
%
% check what kind of changes to make for the parameters
%
thetaChange = 0.001*abs(theta);
thetaAltChange = 0.0001*abs(theta);
IdentMat = eye(p);
dYtilde = zeros(n,p,T);
dvecB = zeros(n*n,p,T);
StopCalculation = 0;
for NumParam=1:p;
   %
   % create a new parameter vector
   %
   ParamDiff = thetaChange(NumParam);
   thetaNew = theta+(ParamDiff*IdentMat(:,NumParam));
   %
   % try to solve the model at this value
   %
   [A,H,R,F,B0,mcode] = SolveDSGEModel(DSGEModel,thetaNew,thetaPositions,ModelParameters,AIMData);
   if mcode~=1;
      %
      % try the alternative 
      %
      ParamDiff = thetaAltChange(NumParam);
      thetaNew = theta+(ParamDiff(NumParam)*IdentMat(:,NumParam));
      [A,H,R,F,B0,mcode] = SolveDSGEModel(DSGEModel,thetaNew,thetaPositions,ModelParameters,AIMData);
   end;
   if mcode==1;
      %
      % run the Kalman filter
      %
      if length(size(H))==2;
         if isempty(DSGEModel.UnitRootStates)==1;
            if isempty(DSGEModel.X)==0;
               [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilterHt(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);
            else;
               [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilterHt(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);
            end;
         else;
            if isempty(DSGEModel.X)==0;
               [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = 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;
      %
      % add values to dYtilde
      %
      dYtilde(:,NumParam,:) = (1/ParamDiff)*((Y-Yhat)-Ytilde);
      vecBNew = zeros(n*n,T);
      for t=1:T;
         vecBNew(:,t) = vec(MSEY(:,:,t));
      end;
      dvecB(:,NumParam,:) = (1/ParamDiff)*(vecBNew-vecB);
   else;
      StopCalculation = 1;
      txt = ['YADA was unable to compute the numerical partial derivatives of the log-likelihood function since the DSGE model couldn''t be solved for the tested new values of the "' StringTrim(ParameterNames(NumParam,:)) '" parameter.'];
      About(ErrorStr,'information','YADA - Bad Numerical Partial Derivatives',120,500,CurrINI);
      break;
   end;
end;
%
% check if the numerical partial derivatives could be calculated.
%
if StopCalculation==1;
   return;
end;
%
% we can now compute the information matrix
%
InformationMatrix = zeros(p,p);
for t=1:T;
   dYtildet = dYtilde(:,:,t);
   invBt = inv(B(:,:,t));
   dvecBt = dvecB(:,:,t);
   dYtildePart = dYtildet'*invBt*dYtildet;
   dvecBpart = (1/2)*dvecBt'*kron(invBt,invBt)*dvecBt;
   InformationMatrix = InformationMatrix+dYtildePart+dvecBpart;
end;

%
% end of DSGEInformationMatrix.m
%
