function CVDs = DSGECVDDistribution(DSGEModel,CurrINI,AnnualizeData,LevelsData,UseGroups,NumDraws,TotalDraws,CurrChain)
% DSGECVDDistribution: Computes the mean value of conditional variance decompositions as well as
%                      quantiles by loading data from disk.
%
% USAGE:
%
%       CVDs = DSGECVDDistribution(DSGEModel,CurrINI,AnnualizeData,LevelsData,UseGroups,NumDraws,TotalDraws,CurrChain)
%
% REQUIRED INPUT:  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure) with initialization information.
%
%                  AnnualizeData (boolean) that takes on the value 1 if the impulse responses should be annualized
%                                 and 0 otherwise.
%
%                  LevelsData (boolean) that takes on the value 1 if the impulse responses should be computed for
%                                 levels variables and 0 otherwise.
%
%                  UseGroups (boolean) that is 1 if shock groups should be used and 0 otherwise.
%
%                  NumDraws (integer) with the number of effective draws from the posterior distribution of the
%                                 parameters.
%
%                  TotalDraws (integer) with the total 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.
%
% REQUIRED OUTPUT: CVDs (structure) with fields "FEVD" and "Quantiles". The "FEVD" field holds a 3D matrix of dimension
%                                 n x q x h, where n is the number of observed variables, q the number of shocks,
%                                 and h the impulse response horizon. The values are the posterior mean conditional
%                                 variance decomposition. The field "Quantiles" is a vector structure with fields "FEVD" and
%                                 "percent", where "FEVD" holds a 3D matrix of dimension n x q x h with values from
%                                 the percent quantile of the marginal distribution of the conditional variance decomposition
%                                 and "percent" holds the percentile value in percent.
%
%
%                       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: April 4, 2007.
%                        This version: January 23, 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:
%
% * 13-11-2007: Updated the documentation.
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 14-04-2008: The code takes into account that the levels decompositions have one additional decomposition,
%               the long run.
%
% * 17-04-2008: Added the input variable "UseGroups" as well as the code for sorting results for shock groups
%               rather than for the individual shocks.
%
% * 23-04-2008: Updated the documentation.
%
% * 23-05-2008: Updated the documentation.
%
% * 23-01-2009: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize
%
CVDs = [];
NumFiles = ceil(NumDraws/min(NumDraws,1000));
%
% Determine the impulse response horizon
%
if strcmp(DSGEModel.DataFrequency,'q')==1;
   SampleFreq = 4;
elseif strcmp(DSGEModel.DataFrequency,'m')==1;
   SampleFreq = 12;
else;
   SampleFreq = 1;
end;
h = SampleFreq*DSGEModel.IRHorizon;
TypeStr = 'Original';
Addh = 0;
if AnnualizeData==1;
   TypeStr = 'Annualized';
end;
if LevelsData==1;
   TypeStr = 'Levels';
   Addh = 1;
end;
%
if UseGroups==0;
   q = length(DSGEModel.StateShockPositions);
   ShockNames = DSGEModel.StateShockNames;
   IMatrix = eye(q);
else;
   q = max(DSGEModel.ShockGroups);
   ShockNames = char(zeros(size(DSGEModel.ShockGroupNames)));
   for i=1:q;
      ShockNames(i,:) = lower(strrep(DSGEModel.ShockGroupNames(i,:),'_',' '));
   end;
   IMatrix = SubIdentityMatrix(DSGEModel.ShockGroups);
end;
n = DSGEModel.n;
r = length(DSGEModel.StateVariablePositions);
NumQuants = length(DSGEModel.Percentiles);
%
% initialize output
%
CVDs.FEVD = zeros(n,q,h+Addh);
for i=1:NumQuants;
   CVDs.Quantiles(i).FEVD = zeros(n,q,h+Addh);
   CVDs.Quantiles(i).percent = DSGEModel.Percentiles(i);
end;
%
% set up a wait dialog
%
txt = ['Please wait while YADA computes the posterior mean conditional variance decompositions as well as percentile values for the observed variables. The data are loaded sequentially from disk and it may therefore take some time. Current (variable,shock): (' StringTrim(DSGEModel.VariableNames(1,:)) ',' StringTrim(ShockNames(1,:)) ').'];
WaitHandle = WaitDLG(txt,'information',['Conditional Variance Decompositions - ' TypeStr ' Data'],500,150,CurrINI,0);
WaitControls = get(WaitHandle,'UserData');
%
% now we can start loading the data
%
for i=1:n;
   for j=1:q;
      %
      % update the wait dialog
      %
      txt = ['Please wait while YADA computes the posterior mean conditional variance decompositions as well as percentile values for the observed variables. The data are loaded sequentially from disk and it may therefore take some time. Current (variable,shock): (' StringTrim(DSGEModel.VariableNames(i,:)) ',' StringTrim(ShockNames(j,:)) ').'];
      set(WaitControls.text,'String',txt);
      drawnow;
      %
      TempVar = zeros(NumDraws,h+Addh);
      CurrSaves = 0;
      for ThisSave=1:NumFiles;
         file = [DSGEModel.OutputDirectory '\cvd\DSGE-CVD-' TypeStr '-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' num2str(h,'%0.0f') '-' num2str(ThisSave,'%0.0f') '-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '.mat'];
         CVDData = load(file);
         CVD = CVDData.CVD;
         NumSaves = length(CVD);
         %
         % load data from CurrSaves to NumSaves+1-CurrSaves
         %
         for s=1:NumSaves;
            %
            % check if we should study data for shock groups or not.
            %
            if UseGroups==0
               FEVD = CVD(s).FEVD;
            else;
               FEVD = zeros(n,q,h+Addh);
               for k=1:h+Addh;
                  FEVD(:,:,k) = CVD(s).FEVD(:,:,k)*IMatrix;
               end;
            end;
            TempVar(CurrSaves+s,:) = FEVD(i,j,:);
         end;
         CurrSaves = CurrSaves+NumSaves;
         drawnow;
      end;
      %
      % now we sort and compute the mean the data
      %
      for k=1:h+Addh;
         TempVar(:,k) = sort(TempVar(:,k));
         if abs(max(TempVar(:,k)))>1.0e-10;
            CVDs.FEVD(i,j,k) = mean(TempVar(:,k));
         end;
         %
         % now we put the quantile information in its correct place
         %
         for l=1:NumQuants;
            if DSGEModel.Percentiles(l)<50;
               QuantVal = ceil((DSGEModel.Percentiles(l)/100)*NumDraws);
            else;
               QuantVal = floor((DSGEModel.Percentiles(l)/100)*NumDraws);
            end;
            if abs(TempVar(QuantVal,k))>1.0e-10;
               CVDs.Quantiles(l).FEVD(i,j,k) = TempVar(QuantVal,k);
            end;
         end;
         drawnow;
      end;
   end;
end;
%
% close the wait dialog
%
set(WaitControls.text,'UserData','done');
delete(WaitHandle);
drawnow;
pause(0.02);

%
% end of DSGECVDDistribution.m
%
