function DoneCalc = DSGECondCorrEstimation(theta,thetaPostSample,thetaPositions,ModelParameters,SelectedParameters,TotalDraws,CurrChain,DSGEModel,CurrINI)
% DSGECondCorrEstimation: Computes conditional correlations for a given parameter vector based on population
%                         moments.
%
% USAGE:
%
%       DoneCalc = DSGECondCorrEstimation(theta,thetaPostSample,thetaPositions,ModelParameters,SelectedParameters,TotalDraws,CurrChain,DSGEModel,CurrINI)
%
% REQUIRED INPUT:  theta (vector) with the posterior mode estimates of the parameters.
%
%                  thetaPostSample (matrix) of dimension NumDrawsxNumParam with draws from the posterior distribution
%                                  of the original DSGE model 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.
%
%                  SelectedParameters (vector) of dimension NumParam which indicates which parameters can vary from
%                                 one draw to the next (1) and which parameters are fixed at posterior mode (0).
%
%                  TotalDraws (integer) with the total number of parameter draws that could have been used.
%
%                  CurrChain (integer) with the number for the current Markov chain.
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure), with initialization information.
%
% REQUIRED OUTPUT: DoneCalc (boolean) that is 1 if the calculations were completed and 0 otherwise.
%
%
%                       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 13, 2009.
%                        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:
%
% * 22-01-2009: Updated the documentation.
%
% * 13-02-2009: The mcode values 7 and 8 are now properly mentioned when encountered.
%
% * 14-02-2009: Fixed the name string for the progress dialog.
%
% * 17-02-2009: Made sure that the variable "KeepShocks" is stored in the mat-files.
%
% * 09-03-2009: Added quotation marks around paths for the dos copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% intialize
%
DoneCalc = 0;
%
% Prepare variables for selecting the parameters that are updated through the
% posterior draws
%
ScenarioParameters = SelectedParameters .* (1:length(theta));
ScenarioParameters = ScenarioParameters(ScenarioParameters>0);
SelParamStr = strrep(num2str(SelectedParameters),' ','');
thetaMode = theta;
%
% 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)==0;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% create the output directory
%
stat = MakeDir(DSGEModel.OutputDirectory,'condcorr');
if stat~=1;
   txt = ['YADA was for some reason unable to create the directory "' DSGEModel.OutputDirectory '\condcorr". The computation of conditional correlation distribution has therefore been aborted.'];
   About(txt,'information','YADA - Directory Creation Problem',180,500,CurrINI);
   return;
end;
%
% initialize index for saving
%
NumDraws = size(thetaPostSample,1);
SaveAfterDraws = min(NumDraws,1000);
NumFiles = ceil(NumDraws/SaveAfterDraws);
%
% solve the model for the posterior mode estimate
%
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;
if length(size(H))==3;
   %
   % Get the sub-sample
   %
   [DSGEModel.FirstPeriod,DSGEModel.LastPeriod] = CreateSubSample(DSGEModel);
   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 computation of the ' lower(EstType) ' distribution.'];
      About(ErrorStr,'error','Error - Measurement Equation',200,500,CurrINI);
      return;
   end;
end;
[B0,KeepShocks] = RemoveRedundantColumns(B0);
[r,q] = size(B0);
n = DSGEModel.n;
%
% begin the loop
%
abort = '';
LastDraw = 0;
MeanEstimationTime = 0;
ThisSave = 0;
CC = [];
%
% check if we should setup a wait handle
%
if DSGEModel.ShowProgress==1;
   ProgressStructure.title = ['Progress for ' num2str(NumDraws,'%0.0f') ' Parameter Draws'];
   ProgressStructure.facecolor = CurrINI.progress_facecolor;
   ProgressStructure.startfacecolor = CurrINI.progress_startfacecolor;
   ProgressStructure.edgecolor = CurrINI.progress_edgecolor;
   ProgressStructure.bgcolor = CurrINI.progress_bgcolor;
   ProgressStructure.stop = 0;
   ProgressStructure.clock = DSGEModel.ShowProgressClock;
   ProgressStructure.label = 'Mean Estimation Time:';
   %
   ProgressStructure.name = 'Population-Based Conditional Correlation Distribution';
   ProgressStructure.CurrINI = CurrINI;
   WaitHandle = ProgressDLG(0,ProgressStructure);
   set(WaitHandle,'Color',get(0,'defaultuicontrolbackgroundcolor'));
   drawnow;
else;
   txt = ['Please wait while YADA computes the population based conditional correlation distribution. Computations started at: ' StringTrim(datestr(now,14))];
   WaitHandle = WaitDLG(txt,'information','Conditional Correlation Distribution',500,150,CurrINI,0);
   WaitControls = get(WaitHandle,'UserData');
end;
%
for it=1:NumDraws;
   LastDraw = LastDraw+1;
   if DSGEModel.ShowProgress==1;
      abort = get(WaitHandle,'UserData');
      if strcmp(abort,'cancel')==1;
         break;
      else;
         ProgressDLG([it/NumDraws MeanEstimationTime]);
      end;
   end;
   %
   % Solve the model for the current parameter vector
   %
   theta(ScenarioParameters) = thetaPostSample(it,ScenarioParameters)';
   tic;
   [A,H,R,F,B0,status] = SolveDSGEModel(DSGEModel,theta,thetaPositions,ModelParameters,AIMData);
   if status~=1;
      %
      % use posterior mode estimate as a backup
      %
      [A,H,R,F,B0,status] = SolveDSGEModel(DSGEModel,thetaMode,thetaPositions,ModelParameters,AIMData);
   end;
   B0 = B0(:,KeepShocks);
   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 = DSGEModel.FirstPeriod-1+DSGEModel.KalmanFirstObservation:DSGEModel.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;
      CC(LastDraw).Shock(j).CondCorr = SigmaYj./(CorrVec*CorrVec');
      drawnow;
   end;
   %
   MeanEstimationTime = (((it-1)/it)*MeanEstimationTime)+((1/it)*toc);
   %
   % check if we should save to disk
   %
   if (LastDraw==SaveAfterDraws)|(it==NumDraws);
      ThisSave = ThisSave+1;
      LastDraw = 0;
      %
      % set up the file name and then save
      %
      file = [DSGEModel.OutputDirectory '\condcorr\DSGE-CondCorr-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' SelParamStr '-' num2str(ThisSave,'%0.0f') '-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '.mat'];
      save(file,'CC','KeepShocks');
      %
      % restore the data structure
      %
      CC = [];
   end;
   drawnow;
end;
%
% close the wait dialog
%
if DSGEModel.ShowProgress==1;
   if ishandle(WaitHandle)==1;
      set(WaitHandle,'UserData','done');
      close(WaitHandle);
      drawnow;
      pause(0.02);
   end;
else;
   set(WaitControls.text,'UserData','done');
   delete(WaitHandle);
   drawnow;
   pause(0.02);
end;
%
% check if we have computed all output
%
if strcmp(abort,'cancel')==0;
   DoneCalc = 1;
end;

%
% end of DSGECondCorrEstimation.m
%
