function SVEs = DSGESVESimulatedDistribution(DSGEModel,CurrINI,Selector,SelectedParameters,NumDraws,TotalDraws,CurrChain,NumPaths)
% DSGESVESimulatedDistribution:  Computes the mean value of X as well as quantiles by loading data from disk.
%                                Based on simulated data, X is here sample based moments for the observed
%                                variables or the the state variables.
%
% USAGE:
%
%       SVEs = DSGESVESimulatedDistribution(DSGEModel,CurrINI,Selector,SelectedParameters,NumDraws,TotalDraws,CurrChain,NumPaths)
%
% REQUIRED INPUT:  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure) with initialization information.
%
%                  Selector (integer) that is:
%                                            5 - if the function should work with observed variable correlations.
%                                            6 - if the function should work with state variable correlations.
%
%                  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).
%
%                  NumDraws (integer) with the number of effective draws from the posterior distribution of the
%                                     parameters.
%
%                  TotalDraws (integer) with the number of draws from the posterior distribution of the parameters
%                                     that could have been used.
%
%                  CurrChain (integer) with the number for the current Markov chain.
%
%                  NumPaths (integer) with the number of data paths that have been simulated per parameter value.
%
% REQUIRED OUTPUT: SVEs (structure) with fields that depend on Selector. For the observed variable correlations the fields
%                                     are: "MeanY", "StdY", and "CorrY". In addition, it yields the fields "DateMeanY",
%                                     "DataStdY", and "DataCorrY". For the state variable correlations the fields are
%                                     "MeanKsi", "StdKsi" and "CorrKsi". In addition, the vector field "Quantiles" is used,
%                                     with subfields "percent" and field names determined by the Selector.
%
%
%                       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 28, 2008.
%                        This version: January 22, 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:
%
% * 29-02-2008: Added the code for "Selector" value 6.
%
% * 03-03-2008: Made sure that the number of exogenous variables is properly taken into account for the observed
%               variable correlations.
%
% * 01-04-2008: Took into account that "NumFiles" has changed. Moreover, the filename now depends on value for
%               "Selector".
%
% * 23-05-2008: Updated the documentation.
%
% * 22-01-2009: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize the output
%
SVEs = [];
%
% fix the information needed to setup the file names
%
SelParamStr = strrep(num2str(SelectedParameters),' ','');
%
%NumFiles = ceil(NumDraws/min(NumDraws,1000));
NumFiles = ceil(NumDraws/1);
NumQuants = length(DSGEModel.Percentiles);
%
% we need to first determine the size of the data
%
if Selector==5;
   file = [DSGEModel.OutputDirectory '\state\DSGE-OVE-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' SelParamStr '-1-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '-' num2str(NumPaths,'%0.0f') '.mat'];
else;
   file = [DSGEModel.OutputDirectory '\state\DSGE-SVE-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' SelParamStr '-1-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '-' num2str(NumPaths,'%0.0f') '.mat'];
end;
SVEData = load(file);
if Selector<=4;
   SVEs.DateVector = SVEData.DateVector;
end;
%
if Selector==5;
   %
   % we have the observed variable correlations
   %
   Y = SVEData.Y;
   [n T] = size(Y);
   %
   % compute the data based mean, standard deviation, and autocorrelation
   %
   X = SVEData.X;
   k = min(size(X));
   if isempty(X)==0;
      SVEs.Data.Mean = (Y*X')/(X*X');
      %
      % de-mean Y using estimate of the mean from the data rather than A'.
      %
      Y = Y-(SVEs.Data.Mean*X);
   else;
      SVEs.Data.Mean = (1/T)*sum(Y,2);
      Y = Y-(SVEs.Data.Mean*ones(1,T));
   end;
   SVEs.Data.Std = sqrt(diag((1/T)*Y*Y'));
   AutoCovHorizon = (size(SVEData.SVE.CorrData(1,1).CorrY,2)-1)/2;
   %
   SVEs.Data.Corr = zeros((n*(n+1))/2,1+(2*AutoCovHorizon));
   SVEs.Data.Corr(:,AutoCovHorizon+1) = vech(((1/T)*Y*Y')./(SVEs.Data.Std*SVEs.Data.Std'));
   for i=1:AutoCovHorizon;
      Y = SVEData.Y(:,1+i:size(SVEData.Y,2));
      YL = SVEData.Y(:,1:length(SVEData.Y)-i);
      if isempty(X)==0;
         X = SVEData.X(:,1+i:size(SVEData.X,2));
         Y = Y-(((Y*X')/(X*X'))*X);
         XL = SVEData.X(:,1:size(SVEData.X,2)-i);
         YL = YL-(((Y*X')/(X*X'))*XL);
      else;
         Y = Y-(SVEs.Data.Mean*ones(1,size(Y,2)));
         YL = YL--(SVEs.Data.Mean*ones(1,size(YL,2)));
      end;
      Corr = ((1/T)*Y*YL')./(SVEs.Data.Std*SVEs.Data.Std');
      %
      % first vector is lagged and second is leaded.
      %
      SVEs.Data.Corr(:,AutoCovHorizon+1+i) = vech(Corr');
      SVEs.Data.Corr(:,AutoCovHorizon+1-i) = vech(Corr);
   end;
   %
   % initialize model based mean, standard deviation, and autocorrelations
   %
   if k>0;
      SVEs.Mean = zeros(n,k);
   else;
      SVEs.Mean = zeros(n,1);
   end;
   SVEs.Std = zeros(n,1);
   SVEs.Corr = zeros((n*(n+1))/2,1+(2*AutoCovHorizon));
   for i=1:NumQuants;
      if k>0;
         SVEs.Quantiles(i).Mean = zeros(n,k);
      else;
         SVEs.Quantiles(i).Mean = zeros(n,1);
      end;
      SVEs.Quantiles(i).Std = zeros(n,1);
      SVEs.Quantiles(i).Corr = zeros((n*(n+1))/2,1+(2*AutoCovHorizon));
      SVEs.Quantiles(i).percent = DSGEModel.Percentiles(i);
   end;
   %
   % setup a wait dialog
   %
   txt = ['Please wait while YADA computes the distributions for the model based sample means and autocorrelations of the observed variables. The data are loaded sequentially from disk and it may therefore take some time. Current variable: ' StringTrim(DSGEModel.VariableNames(1,:)) '.'];
   WaitHandle = WaitDLG(txt,'information','Autocorrelation Distribution',500,200,CurrINI,0);
   WaitControls = get(WaitHandle,'UserData');
   %
   for i=1:n;
      %
      % setup a wait dialog
      %
      txt = ['Please wait while YADA computes the distributions for the model based sample means and autocorrelations of the observed variables. The data are loaded sequentially from disk and it may therefore take some time. Current variable: ' StringTrim(DSGEModel.VariableNames(i,:)) '.'];
      set(WaitControls.text,'String',txt);
      drawnow;
      %
      % set up a matrix for the data
      %
      if k>0;
         TempMean = zeros(NumDraws*NumPaths,k);
      else;
         TempMean = zeros(NumDraws*NumPaths,1);
      end;
      TempStd = zeros(NumDraws*NumPaths,1);
      CurrSaves = 0;
      for ThisSave=1:NumFiles;
         file = [DSGEModel.OutputDirectory '\state\DSGE-OVE-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' SelParamStr '-' num2str(ThisSave,'%0.0f') '-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '-' num2str(NumPaths,'%0.0f') '.mat'];
         SVEData = load(file);
         SVE = SVEData.SVE;
         NumSaves = size(SVE.CorrData,1);
         %
         % load data from CurrSaves to NumSaves+1-CurrSaves
         %
         for s=1:NumSaves;
            for p=1:NumPaths;
               TempMean(CurrSaves+p+((s-1)*NumPaths),:) = SVE.CorrData(s,p).MeanY(i,:);
               TempStd(CurrSaves+p+((s-1)*NumPaths),1) = SVE.CorrData(s,p).StdY(i,1);
            end;
            drawnow;
         end;
         CurrSaves = CurrSaves+(NumSaves*NumPaths);
         drawnow;
      end;
      %
      % we can now calculate the mean
      %
      SVEs.Mean(i,:) = (1/(NumDraws*NumPaths))*sum(TempMean,1);
      SVEs.Std(i,1) = (1/(NumDraws*NumPaths))*sum(TempStd,1);
      %
      % now we can sort these and create the quantiles from the sorted
      %
      for j=1:size(TempMean,2);
         TempMeanj = sort(TempMean(:,j));
         if j==1;
            TempStd = sort(TempStd);
         end;
         for l=1:NumQuants;
            if DSGEModel.Percentiles(l)<50;
               QuantVal = ceil((DSGEModel.Percentiles(l)/100)*NumDraws*NumPaths);
            else;
               QuantVal = floor((DSGEModel.Percentiles(l)/100)*NumDraws*NumPaths);
            end;
            SVEs.Quantiles(l).Mean(i,j) = TempMeanj(QuantVal,1);
            if j==1;
               SVEs.Quantiles(l).Std(i,1) = TempStd(QuantVal,1);
            end;
         end;
         drawnow;
      end;
   end;
   %
   % now we need to deal with the autocorrelations
   %
   SVEs.AutoCovHorizon = (-AutoCovHorizon:1:AutoCovHorizon);
   for j=1:(2*AutoCovHorizon)+1;
      %
      % setup a wait dialog
      %
      txt = ['Please wait while YADA computes the distributions for the model based sample means and autocorrelations of the observed variables. The data are loaded sequentially from disk and it may therefore take some time. Current autocorrelation horizon: ' int2str(j-(AutoCovHorizon+1)) '.'];
      set(WaitControls.text,'String',txt);
      drawnow;
      %
      TempCorr = zeros(NumDraws*NumPaths,(n*(n+1))/2);
      CurrSaves = 0;
      for ThisSave=1:NumFiles;
         file = [DSGEModel.OutputDirectory '\state\DSGE-OVE-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' SelParamStr '-' num2str(ThisSave,'%0.0f') '-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '-' num2str(NumPaths,'%0.0f') '.mat'];
         SVEData = load(file);
         SVE = SVEData.SVE;
         NumSaves = size(SVE.CorrData,1);
         %
         % load data from CurrSaves to NumSaves+1-CurrSaves
         %
         for s=1:NumSaves;
            for p=1:NumPaths;
               TempCorr(CurrSaves+p+((s-1)*NumPaths),:) = SVE.CorrData(s,p).CorrY(:,j)';
            end;
            drawnow;
         end;
         CurrSaves = CurrSaves+(NumSaves*NumPaths);
         drawnow;
      end;
      %
      % compute the mean
      %
      SVEs.Corr(:,j) = (1/(NumDraws*NumPaths))*sum(TempCorr,1)';
      %
      % sort the data and determine the quantile values
      %
      TempCorr = sort(TempCorr);
      for l=1:NumQuants;
         if DSGEModel.Percentiles(l)<50;
            QuantVal = ceil((DSGEModel.Percentiles(l)/100)*NumDraws*NumPaths);
         else;
            QuantVal = floor((DSGEModel.Percentiles(l)/100)*NumDraws*NumPaths);
         end;
         SVEs.Quantiles(l).Corr(:,j) = TempCorr(QuantVal,:)';
      end;
      drawnow;
   end;
elseif Selector==6;
   %
   % we have the state variable correlations
   %
   r = length(SVEData.SVE.CorrState(1,1).MeanKsi);
   SVEs.MeanKsi = zeros(r,1);
   SVEs.StdKsi = zeros(r,1);
   SVEs.CorrKsi = zeros(r,r);
   for i=1:NumQuants;
      SVEs.Quantiles(i).MeanKsi = zeros(r,1);
      SVEs.Quantiles(i).StdKsi = zeros(r,1);
      SVEs.Quantiles(i).CorrKsi = zeros(r,r);
      SVEs.Quantiles(i).percent = DSGEModel.Percentiles(i);
   end;
   %
   % setup a wait dialog
   %
   txt = ['Please wait while YADA computes the mean state variable correlations for model-based simulated data as well as percentile values for the state variable correlations. The data are loaded sequentially from disk and it may therefore take some time. Current variable: ' StringTrim(DSGEModel.StateVariableNames(1,:)) '.'];
   WaitHandle = WaitDLG(txt,'information','State Variable Distribution',500,200,CurrINI,0);
   WaitControls = get(WaitHandle,'UserData');
   %
   for i=1:r;
      %
      % update the wait dialog
      %
      txt = ['Please wait while YADA computes the mean state variable correlations for model-based simulated data as well as percentile values for the state variable correlations. The data are loaded sequentially from disk and it may therefore take some time. Current variable: ' StringTrim(DSGEModel.StateVariableNames(i,:)) '.'];
      set(WaitControls.text,'String',txt);
      drawnow;
      %
      % set up matrices for storing the data
      %
      TempMen = zeros(NumDraws*NumPaths,1);
      TempStd = zeros(NumDraws*NumPaths,1);
      TempCorr = zeros(NumDraws*NumPaths,r);
      CurrSaves = 0;
      for ThisSave=1:NumFiles;
         file = [DSGEModel.OutputDirectory '\state\DSGE-SVE-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' SelParamStr '-' num2str(ThisSave,'%0.0f') '-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '-' num2str(NumPaths,'%0.0f') '.mat'];
         SVEData = load(file);
         SVE = SVEData.SVE;
         NumSaves = size(SVE.CorrState,1);
         %
         % load data from CurrSaves to NumSaves+1-CurrSaves
         %
         for s=1:NumSaves;
            for p=1:NumPaths;
               TempMean(CurrSaves+p+((s-1)*NumPaths),:) = SVE.CorrState(s,p).MeanKsi(i,1);
               TempStd(CurrSaves+p+((s-1)*NumPaths),:) = SVE.CorrState(s,p).StdKsi(i,1);
               TempCorr(CurrSaves+p+((s-1)*NumPaths),:) = SVE.CorrState(s,p).CorrKsi(i,:);
            end;
            drawnow;
         end;
         CurrSaves = CurrSaves+(NumSaves*NumPaths);
         drawnow;
      end;
      %
      % we can now calculate the mean
      %
      SVEs.MeanKsi(i,1) = (1/(NumDraws*NumPaths))*sum(TempMean);
      SVEs.StdKsi(i,1) = (1/(NumDraws*NumPaths))*sum(TempStd);
      SVEs.CorrKsi(i,:) = (1/(NumDraws*NumPaths))*sum(TempCorr,1);
      %
      % now we can sort these and create the quantiles from the sorted
      %
      TempMean = sort(TempMean);
      TempStd = sort(TempStd);
      TempCorr = sort(TempCorr);
      for l=1:NumQuants;
         if DSGEModel.Percentiles(l)<50;
            QuantVal = ceil((DSGEModel.Percentiles(l)/100)*NumDraws*NumPaths);
         else;
            QuantVal = floor((DSGEModel.Percentiles(l)/100)*NumDraws*NumPaths);
         end;
         SVEs.Quantiles(l).MeanKsi(i,1) = TempMean(QuantVal,1);
         SVEs.Quantiles(l).StdKsi(i,1) = TempStd(QuantVal,1);
         SVEs.Quantiles(l).CorrKsi(i,:) = TempCorr(QuantVal,:);
      end;
      drawnow;
   end;
end;
%
% close the wait dialog
%
set(WaitControls.text,'UserData','done');
delete(WaitHandle);
drawnow;
pause(0.02);

%
% end of DSGESVESimulatedDistribution.m
%
