function [CondCorr,status] = DSGEConditionalCorrsTheta(theta,thetaPositions,ModelParameters,NumPaths,EstType,DSGEModel,CurrINI,SimulateData,FirstPeriod,LastPeriod)
% DSGEConditionalCorrsTheta: Computes conditional correlations for a given parameter vector based on population
%                                   moments or sample moments.
%
% USAGE:
%
%       [CondCorr,status] = DSGEConditionalCorrsTheta(theta,thetaPositions,ModelParameters,NumPaths,EstType,DSGEModel,CurrINI,SimulateData,FirstPeriod,LastPeriod)
%
% REQUIRED INPUT:  theta (vector) with the posterior mode estimates of the parameters.
%
%                  thetaPositions (structure) where thetaPositions(i).parameter is a string vector giving the name
%                                 of the model parameter in position i of theta. The field "distribution" has a string
%                                 vector giving the assumed prior distribution of the parameter.
%
%                  ModelParameters (structure) whose fields have names given by the parameter names in the
%                                 DSGE model.
%
%                  NumPaths (integer) with the number of data paths to simulate per parameter value.
%
%                  EstType (string vector) which determines the type of parameter values that are used. This is
%                                 either "Posterior Mode" or "Initial Values".
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure), with initialization information.
%
%                  SimulateData (boolean) that indicated if sample correlations should be computed (1) instead of
%                                 population correlations.
%
%                  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: CondCorr (structure) with fields Mean and Quantiles and ShockNames.
%
% OPTIONAL OUTPUT: status (integer) Indicates whether the solution of the model is unique or not.
%                                 status is 1 when unique, 0 if not computed, and various other values when
%                                 there are problems with the solution. Specifically:
%
%                                      status = 2:      Roots not correctly computed by real_schur
%                                      status = 3:      Too many big roots
%                                      status = 35:     Too many big roots, and q(:,right) is singular
%                                      status = 4:      Too few big roots
%                                      status = 45:     Too few big roots, and q(:,right) is singular
%                                      status = 5:      q(:,right) is singular
%                                      status = 61:     Too many exact shiftrights
%                                      status = 62:     Too many numeric shiftrights
%                                      status = 7:      The companion matrix has infinite or NaN entries.
%                                      status = 8:      The function "compute_aim_matrices" returns complex
%                                                       numbers.
%                                      else             Return code not properly specified
%
%
%                       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: January 9, 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:
%
% * 12-01-2009: Added the wait dialog when simulating sample based conditional correlations.
%
% * 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 output
%
CondCorr = [];
%
% 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)==0;
   if nargout==2;
      status = 0;
   end;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% now we try to solve the model
%
if nargout==1;
   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','YADA - AiM Solver Error',200,500,CurrINI);
      return;
   end;
else;
   ErrorStr = '';
   try;
      [A,H,R,F,B0,status] = 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 status~=1;
      return;
   end;
end;
%
if length(size(H))==3;
   if size(H,3)<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(LastPeriod) '). YADA has therefore aborted from the conditional correlations.'];
      About(ErrorStr,'error','Error - Measurement Equation',200,500,CurrINI);
      return;
   end;
end;
%
% make sure that only non-zero shocks are kept
%
[B0,KeepShocks] = RemoveRedundantColumns(B0);
n = DSGEModel.n;
k = DSGEModel.k;
[r,q] = size(B0);
CondCorr(1).ShockNames = DSGEModel.StateShockNames(KeepShocks,:);
%
% check if population or sample based correlations should be computed
%
if SimulateData==0;
   for j=1:q;
      %
      % compute the conditional state population covariance matrix when there is
      % only shock i using the doubling algorithm
      %
      B0j = B0(:,j);
      [SigmaKsij,stat] = DoublingAlgorithmLyapunov(F,B0j*B0j',DSGEModel.DAMaximumIterationsValue*100,10^(-(DSGEModel.DAToleranceValue+1)));
      if stat==1;
         %
         % we need to solve the problem analytically
         %
         vecSigmaKsij = (eye(r*r)-kron(F,F))\vec(B0j*B0j');
         SigmaKsij = reshape(vecSigmaKsij,r,r);
      end;
      %
      % compute the conditional covariance matrix for the observed variables
      %
      if length(size(H))==2;
         SigmaYj = H'*SigmaKsij*H;
      else;
         SigmaYj = zeros(n,n);
         Sample = FirstPeriod-1+DSGEModel.KalmanFirstObservation:LastPeriod;
         for t=Sample;
            SigmaYj = SigmaYj+((1/length(Sample))*(H(:,:,t)'*SigmaKsij*H(:,:,t)));
         end;
      end;
      CorrVec = sqrt(abs(diag(SigmaYj)));
      for i=1:n;
         if CorrVec(i)<1.0e-12;
            CorrVec(i) = 1;
         end;
      end;
      CondCorr(j).Mean = SigmaYj./(CorrVec*CorrVec');
   end;
else;
   %
   % initialize the random number generator
   %
   if DSGEModel.RandomNumberValue==1;
      randn('state',0);
      rand('state',0);
   else;
      randn('state',sum(100*clock));
      rand('state',sum(100*clock));
   end;
   T = LastPeriod-FirstPeriod-DSGEModel.KalmanFirstObservation+2;
   %
   % set the InitialStateVector value
   %
   if DSGEModel.UseOwnInitialState==1;
      if length(DSGEModel.InitialStateValues)==size(AIMData.endog,1);
         Ksi0 = DSGEModel.InitialStateValues(DSGEModel.StateVariablePositions);
      else;
         Ksi0 = zeros(length(DSGEModel.StateVariablePositions),1);
      end;
   else;
      Ksi0 = zeros(length(DSGEModel.StateVariablePositions),1);
   end;
   Dn = DuplicationMatrix(n);
   abort = '';
   MeanEstimationTime = 0;
   StopEstimation = 0;
   %
   % setup a wait dialog
   %
   ShockNames = DSGEModel.StateShockNames(KeepShocks,:);
   txt = ['Please wait while YADA computes the distributions for the model based sample conditional correlations of the observed variables. Current shock: ' StringTrim(ShockNames(1,:)) '.'];
   WaitHandle = WaitDLG(txt,'information',['Conditional Correlations - ' EstType],500,200,CurrINI,1);
   WaitControls = get(WaitHandle,'UserData');
   %
   for j=1:q;
      %
      % update the wait dialog
      %
      if isempty(WaitControls.cancel)==0;
         %
         % check if the cancel button was pressed
         %
         if strcmp(get(WaitControls.cancel,'UserData'),'cancel')==1;
            StopEstimation = 1;
            break;
         end;
      end;
      %
      txt = ['Please wait while YADA computes the distributions for the model based sample conditional correlations of the observed variables. Current shock: ' StringTrim(ShockNames(j,:)) '.'];
      set(WaitControls.text,'String',txt);
      drawnow;
      pause(0.02);
      %
      % compute the conditional state population covariance matrix when there is
      % only shock i using the doubling algorithm
      %
      B0j = B0(:,j);
      [SigmaKsij,stat] = DoublingAlgorithmLyapunov(F,B0j*B0j',DSGEModel.DAMaximumIterationsValue*100,10^(-(DSGEModel.DAToleranceValue+1)));
      if stat==1;
         %
         % we need to solve the problem analytically
         %
         vecSigmaKsij = (eye(r*r)-kron(F,F))\vec(B0j*B0j');
         SigmaKsij = reshape(vecSigmaKsij,r,r);
      end;
      [S,L] = SchurDecomposition(SigmaKsij);
      %
      % we can now simulate the paths
      %
      SimCondCorr = [];
      for p=1:NumPaths;
         tic;
         %
         % draw the economic shocks
         %
         etaSim = zeros(q,T);
         etaSim(j,:) = randn(1,T);
         %
         % draw an initial state
         %
         StateSim = zeros(r,T+1);
         StateSim(:,1) = Ksi0+(S*sqrt(L)*randn(size(L,2),1));
         %
         % generate the state matrix
         %
         for t=1:T;
            StateSim(:,t+1) = (F*StateSim(:,t))+(B0*etaSim(:,t));
         end;
         drawnow;
         %
         % remove the initial value of the state vector
         %
         StateSim = StateSim(:,2:T+1);
         %
         % compute the sample covariance matrix
         %
         if length(size(H))==2;
            SigmaYjs = (1/T)*H'*StateSim*StateSim'*H;
         else;
            SigmaYjs = zeros(n,n);
            for t=1:T;
               SigmaYjs = SigmaYjs+((1/T)*H(:,:,t+FirstPeriod+DSGEModel.KalmanFirstObservation-2)'*StateSim(:,t)*StateSim(:,t)'*H(:,:,t+FirstPeriod+DSGEModel.KalmanFirstObservation-2));
            end;
         end;
         CorrVec = sqrt(abs(diag(SigmaYjs)));
         for i=1:n;
            if CorrVec(i)<1.0e-12;
               CorrVec(i) = 1;
            end;
         end;
         SimCondCorr(p).Corr = SigmaYjs./(CorrVec*CorrVec');
         MeanEstimationTime = (((p-1)/p)*MeanEstimationTime)+((1/p)*toc);
         drawnow;
      end;
      %
      % compute the mean and the quantiles
      %
      if strcmp(abort,'cancel')==0;
         NumQuants = length(DSGEModel.Percentiles);
         %
         CondCorr(j).Mean = zeros(n,n);
         for p=1:NumPaths;
            CondCorr(j).Mean = CondCorr(j).Mean+((1/NumPaths)*SimCondCorr(p).Corr);
         end;
         %
         % collect all the unique data in a matrix
         %
         if n>1;
            QuantileData = zeros(p,n*(n+1)/2);
            for p=1:NumPaths;
               QuantileData(p,:) = vech(SimCondCorr(p).Corr)';
            end;
            %
            % sort each row
            %
            for i=1:n*(n+1)/2;
               QuantileData(:,i) = sort(QuantileData(:,i));
            end;
            %
            % deal with the quantiles
            %
            for i=1:NumQuants;
               %
               % initialize
               %
               CondCorr(j).Quantiles(i).percent = DSGEModel.Percentiles(i);
               if DSGEModel.Percentiles(i)<50;
                  QuantVal = ceil((DSGEModel.Percentiles(i)/100)*NumPaths);
               else;
                  QuantVal = floor((DSGEModel.Percentiles(i)/100)*NumPaths);
               end;
               CondCorr(j).Quantiles(i).Mean = reshape(Dn*QuantileData(QuantVal,:)',n,n);
            end;
         end;
      end;
   end;
   if StopEstimation==1;
      CondCorr = [];
   end;
   %
   % close the wait dialog
   %
   set(WaitControls.text,'UserData','done');
   delete(WaitHandle);
   drawnow;
   pause(0.02);
end;

%
% end of DSGEConditionalCorrsTheta.m
%
