function [CorrData,status] = CalculateDSGECorrelations(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI,FirstPeriod,LastPeriod)
% CalculateDSGECorrelations: Computes the state-space model determined autocorrelations for the obserserved
%                            variables.
%
% USAGE:
%
%       [CorrData,status] = CalculateDSGECorrelations(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI,FirstPeriod,LastPeriod)
%
% 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.
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure) with initialization information.
%
%                  FirstPeriod (integer), with the first period in the sample when the training sample size
%                                 for the Kalman filter is not taken into account.
%
%                  LastPeriod (integer), with the last period in the sample.
%
% REQUIRED OUTPUT: CorrData (structure) with main fields "StateSpace" and "Data". These fields have sub-fields
%                                 "Mean", "Cov", and "Corr". The main field "StateSpace" refers to the
%                                 results based on the state-space model, while the main field "Data" refers to
%                                 the results based on the data. Notice that the autocovariance and autocorrelation
%                                 results using the data are calculated by de-meaning using the mean from the
%                                 state-space model. This should make the results across the two approaches
%                                 easier to compare.
%
% OPTIONAL OUTPUT: status (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
%
% NOTE: Two output arguments should only be used when theta has the initial values.
%
%
%                       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: October 25, 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:
%
% * 26-10-2006: Added the data determined autocovariances and autocorrelations. The inputs "FirstPeriod" and "LastPeriod"
%               were also added to make sure that the data determined results are based on the same sample as the state-space
%               model determined results.
%
% * 31-10-2006: Made sure that the value of DSGEModel.DAToleranceValue is sent in the appropriate form to the
%               "DoublingAlgorithmLyapunov" function. Added the output field "StateVars" with sub-fields "Variance"
%               and "Correlations". The data stored here is the vector of state variable variances and the matrix with
%               contemporaneous state variable correlations.
%
% * 07-11-2006: The data covariances are now based on the data mean rather than the model mean.
%
% * 17-11-2006: Fixed a bug in Y and YL when YL is constructed for the autocovariances and there are exogenous variables.
%
% * 23-11-2006: Added try-catch statement around eval call for measurement equation function.
%
% * 12-12-2006: Made sure that the the same sample mean is removed from Y and the lagged data matrix YL.
%
% * 13-11-2007: Updated the documentation.
%
% * 29-02-2008: The autocorrelations are now written to a 2D matrix via instead of a 3D matrix. The autocovariances
%               are no longer sent to the output structure.
%
% * 01-04-2008: Change the AutoCovHorizon value to "min([SampleFreq*3 DSGEModel.T-20 18])".
%
% * 16-04-2008: The code now checks if the model has a unique convergent solution also at the posterior mode
%               estimates.
%
% * 16-05-2008: Took the "DSGEModel.AIMTolerance" input for AiMSolver into account.
%
% * 23-05-2008: Updated the documentation.
%
% * 27-05-2008: The mcode value 7 is now supported.
%
% * 03-07-2008: Made sure that the measurement equation receives only the active state variable names.
%               Added try-catch code for AiMSolver.
%
% * 09-07-2008: Added code to allow for a time-varying measurement matrix H.
%
% * 22-01-2009: Updated the documentation.
%
% * 30-01-2009: Made sure that status is only set to mcode when the latter exists.
%
% * 09-02-2009: Changed the input "length(DSGEModel.StateVariablePositions)" to "size(F,1)" when calling the
%               measurement equation function.
%
% * 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
%
CorrData = [];
%
% set up the ModelParameters structure based on theta and thetaPositions
%
if nargout==1;
   %
   % Make sure that the ModelParameter structure takes the correct values.
   % This is not necessary to do when initial values are given to the function.
   % In that case, more substantial checking is required!
   %
   ModelParameters = ThetaToModelParameters(ModelParameters,theta,thetaPositions);
   %
   % 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;
   %
   % Update any parameteters that need updating
   %
   if FileExist(DSGEModel.UpdateParameterFile)==1;
      eval(['ModelParameters = ' GetName(DSGEModel.UpdateParameterFile) '(ModelParameters);']);
   end;
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;
   %
   % See if we can solve the model
   %
   if nargout==1;
      [mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
   else;
      try;
         [mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
      catch;
         if exist('mcode')>0;
            status = mcode;
         else;
            status = 0;
         end;
         ErrorStr = ['There was an error running AiM. YADA received the following message: ' lasterr];
         About(ErrorStr,'error','AiM Error',300,500,CurrINI);
         drawnow;
         return;
      end;
   end;
   if nargout==2;
      status = mcode;
      if status~=1;
         return;
      end;
   else;
      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','YADA - AiM Solver Error',200,500,CurrINI);
         return;
      end;
   end;
   %
   % Convert to the state equation matrices
   %
   [F,B0] = AiMtoStateSpace(BMatrix,SMatrix,DSGEModel.StateVariablePositions,DSGEModel.StateShockPositions);
   %
   % obtain the measurement equation matrices
   %
   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;
      %
      % display an error message
      %
      if nargout==2;
         status = 0;
      end;
      ErrorStr = ['Error found in "' GetName(DSGEModel.MeasurementEquationFile) '.m". Message caught is: ' lasterr];
      About(ErrorStr,'error','Error - Measurement Equation',150,500,CurrINI);
      return;
   end;
   %
   % Determine the autocovariance horizon
   %
   if strcmp(DSGEModel.DataFrequency,'q')==1;
      SampleFreq = 4;
   elseif strcmp(DSGEModel.DataFrequency,'m')==1;
      SampleFreq = 12;
   else;
      SampleFreq = 1;
   end;
   DSGEModel.Y = DSGEModel.Y(:,FirstPeriod-1+DSGEModel.KalmanFirstObservation:LastPeriod);
   if min(size(DSGEModel.X))>0;
      DSGEModel.X = DSGEModel.X(:,FirstPeriod-1+DSGEModel.KalmanFirstObservation:LastPeriod);
   end;
   DSGEModel.T = size(DSGEModel.Y,2);
   AutoCovHorizon = min([SampleFreq*3 DSGEModel.T-20 18]);
   %
   % Now we can compute the autocovariances
   %
   if min(size(DSGEModel.X))>0;
      CorrData.StateSpace.Mean = A';
   else;
      CorrData.StateSpace.Mean = zeros(DSGEModel.n,1);
   end;
   %
   % calculate the unconditional variance of the state variables
   %
   [Sigma,stat] = DoublingAlgorithmLyapunov(F,B0*B0',DSGEModel.DAMaximumIterationsValue*100,10^(-(DSGEModel.DAToleranceValue+1)));
   %
   r = size(H,1);
   n = size(H,2);
   CorrData.StateVars.StDev = sqrt(diag(Sigma));
   CorrData.StateVars.Correlations = Sigma./(CorrData.StateVars.StDev*CorrData.StateVars.StDev');
   CorrData.StateSpace.Corr = zeros((n*(n+1))/2,(2*AutoCovHorizon)+1);
   if length(size(H))==2;
      AutoCov = (H'*Sigma*H)+R;
      CorrVec = sqrt(diag(AutoCov));
      CorrData.StateSpace.Cov = AutoCov;
      CorrData.StateSpace.Corr(:,AutoCovHorizon+1) = vech(AutoCov./(CorrVec*CorrVec'));
      CorrData.AutoCovHorizon = (-AutoCovHorizon:1:AutoCovHorizon);
      FPower = F;
      for i=1:AutoCovHorizon;
         Corr = (H'*FPower*Sigma*H)./(CorrVec*CorrVec');
         CorrData.StateSpace.Corr(:,AutoCovHorizon+1+i) = vech(Corr');
         CorrData.StateSpace.Corr(:,AutoCovHorizon+1-i) = vech(Corr);
         FPower = FPower*F;
      end;
   else;
      AutoCov = zeros(n,n);
      Sample = FirstPeriod-1+DSGEModel.KalmanFirstObservation:LastPeriod;
      for t=Sample;
         AutoCov = AutoCov+((1/length(Sample))*(H(:,:,t)'*Sigma*H(:,:,t))+R);
      end;
      CorrVec = sqrt(diag(AutoCov));
      CorrData.StateSpace.Cov = AutoCov;
      CorrData.StateSpace.Corr(:,AutoCovHorizon+1) = vech(AutoCov./(CorrVec*CorrVec'));
      CorrData.AutoCovHorizon = (-AutoCovHorizon:1:AutoCovHorizon);
      FPower = F;
      for i=1:AutoCovHorizon;
         AutoCov = zeros(n,n);
         for t=FirstPeriod-1+DSGEModel.KalmanFirstObservation:LastPeriod-i;
            AutoCov = AutoCov+((1/(length(Sample)-i))*(H(:,:,t+i)'*FPower*Sigma*H(:,:,t)));
         end;
         Corr = AutoCov./(CorrVec*CorrVec');
         CorrData.StateSpace.Corr(:,AutoCovHorizon+1+i) = vech(Corr');
         CorrData.StateSpace.Corr(:,AutoCovHorizon+1-i) = vech(Corr);
         FPower = FPower*F;
      end;
   end;
   %
   % Do the same with the actual data
   %
   Y = DSGEModel.Y;
   if min(size(DSGEModel.X))>0;
      X = DSGEModel.X;
      CorrData.Data.Mean = (Y*X')/(X*X');
      %
      % de-mean Y using estimate of the mean from the data rather than A'.
      %
      Y = Y-(CorrData.Data.Mean*X);
   else;
      CorrData.Data.Mean = (1/T)*sum(Y,2);
      Y = Y-(CorrData.Data.Mean*ones(1,T));
   end;
   T = length(Y);
   %
   AutoCov = (1/T)*Y*Y';
   CorrVec = sqrt(diag(AutoCov));
   CorrData.Data.Cov = AutoCov;
   CorrData.Data.Corr = zeros((n*(n+1))/2,1+(2*AutoCovHorizon));
   CorrData.Data.Corr(:,AutoCovHorizon+1) = vech(AutoCov./(CorrVec*CorrVec'));
   %
   for i=1:AutoCovHorizon;
      Y = DSGEModel.Y(:,1+i:length(DSGEModel.Y));
      YL = DSGEModel.Y(:,1:length(DSGEModel.Y)-i);
      T = length(Y);
      if min(size(DSGEModel.X))>0;
         X = DSGEModel.X(:,1+i:length(DSGEModel.Y));
         XL = DSGEModel.X(:,1:length(DSGEModel.Y)-i);
         NewMeanY = (Y*X')/(X*X');
         Y = Y-(NewMeanY*X);
         YL = YL-(NewMeanY*XL);
      else;
         Y = Y-(CorrData.Data.Mean*ones(1,T));
         YL = YL-(CorrData.Data.Mean*ones(1,T));
      end;
      %
      AutoCov = (1/T)*Y*YL';
      CorrVec = sqrt(diag((1/T)*Y*Y'));
      CorrVecL = sqrt(diag((1/T)*YL*YL'));
      Corr = AutoCov./(CorrVec*CorrVecL');
      CorrData.Data.Corr(:,AutoCovHorizon+1+i) = vech(Corr');
      CorrData.Data.Corr(:,AutoCovHorizon+1-i) = vech(Corr);
   end;
else;
   if nargout==2;
      status = 0;
   end;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
end;

%
% end of CalculateDSGECorrelations.m
%
