function HFED = HFEDDistribution(DSGEModel,CurrINI,h,NumDraws,TotalDraws,CurrChain)
% HFEDDistribution: Computes the mean value of the historical forecast error decompositions as well as
%                   quantiles by loading data from disk.
%
% USAGE:
%
%       HFED = HFEDDistribution(DSGEModel,CurrINI,h,NumDraws,TotalDraws,CurrChain)
%
% REQUIRED INPUT:  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
%                  CurrINI (structure) with intialization information.
%
%                  h (positive integer), with the forecast horizon.
%
%                  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.
%
% REQUIRED OUTPUT: HFED (structure) with fields "epstth", an nxT matrix with the mean h-steps ahead historical
%                           forecast errors, "wthT", an nxT matrix with the mean h-steps ahead smooth estimate
%                           of the measurement error, "KsiErr", nxT matrix with the mean of the share of the
%                           historical forecast error due to the state prediction error, "etathT", nxTxq matrix
%                           with the mean of the shares of the historical forecast error due to the structura
%                           shocks. Furthermore, "R" (nxn) holds the posterior mode estimate of the measurement
%                           error covariance matrix, "KeepVar" is a qx1 vector with intgeres indicating the
%                           nonzero structural shocks, "forecastHorizon" postive integer with the forecast
%                           horizon, and "Quantiles" a vector structura with the quantile values for "epstth",
%                           "wthT", "KsiErr", and "etathT".
%
%
%                       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: May 7, 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:
%
% * 08-05-2008: Finished the function.
%
% * 13-05-2008: Fixed a bug when updating the CurrSaves variable.
%
% * 23-05-2008: Updated the documentation.
%
% * 22-01-2009: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize the output
%
HFED.ForecastHorizon = h;
NumFiles = ceil(NumDraws/min(NumDraws,100));
NumQuants = length(DSGEModel.Percentiles);
%
% load one file to record the size of the data
%
file = [DSGEModel.OutputDirectory '\hfed\DSGE-HFED-' int2str(h) '-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-1-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '.mat'];
FEData = load(file);
[n,T] = size(FEData.HFED(1).HistDecompStructure.epstth);
q = length(FEData.KeepShocks);
HFED.R = FEData.Rmode;
HFED.KeepVar = FEData.KeepShocks;
ShockNames = DSGEModel.StateShockNames(FEData.KeepShocks,:);
%
HFED.epstth = zeros(n,T);
HFED.wthT = zeros(n,T);
HFED.KsiErr = zeros(n,T);
HFED.etathT = zeros(n,T,q);
for i=1:NumQuants;
   HFED.Quantiles(i).epstth = zeros(n,T);
   HFED.Quantiles(i).wthT = zeros(n,T);
   HFED.Quantiles(i).KsiErr = zeros(n,T);
   HFED.Quantiles(i).etathT = zeros(n,T,q);
   HFED.Quantiles(i).percent = DSGEModel.Percentiles(i);
end;
%
% start with the nxT matrices
%
%
% set up a wait dialog
%
txt = ['Please wait while YADA computes the mean historical forecast error decompsotions as well as percentile values. 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','Forecast Error Decompositions - Distribution',500,150,CurrINI,0);
WaitControls = get(WaitHandle,'UserData');
for i=1:n;
   txt = ['Please wait while YADA computes the mean historical forecast error decompsotions as well as percentile values. The data are loaded sequentially from disk and it may therefore take some time. Current (variable,shock): (' StringTrim(DSGEModel.VariableNames(i,:)) ',' StringTrim(ShockNames(1,:)) ').'];
   set(WaitControls.text,'String',txt);
   drawnow;
   %
   % initialize the data
   %
   TmpEps = zeros(NumDraws,T);
   TmpKsi = zeros(n,T);
   TmpShock = zeros(n,T);
   if rank(HFED.R)>0;
      TmpwthT = zeros(NumDraws,T);
   end;
   CurrSaves = 0;
   %
   % load the data
   %
   for ThisSave=1:NumFiles;
      file = [DSGEModel.OutputDirectory '\hfed\DSGE-HFED-' int2str(h) '-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' num2str(ThisSave,'%0.0f') '-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '.mat'];
      HFEDData = load(file);
      FED = HFEDData.HFED;
      NumSaves = length(FED);
      for s=1:NumSaves;
         TmpEps(CurrSaves+s,:) = FED(s).HistDecompStructure.epstth(i,:);
         TmpKsi(CurrSaves+s,:) = FED(s).HistDecompStructure.KsiErr(i,:);
         TmpShock(CurrSaves+s,:) = FED(s).HistDecompStructure.etathT(i,:,1);
         if rank(HFED.R)>0;
            TmpwthT(CurrSaves+s,:) = FED(s).HistDecompStructure.wthT(i,:);
         end;
      end;
      CurrSaves = CurrSaves+NumSaves;
      drawnow;
   end;
   %
   % compute the mean
   %
   HFED.epstth(i,:) = (1/NumDraws)*sum(TmpEps);
   HFED.KsiErr(i,:) = (1/NumDraws)*sum(TmpKsi);
   HFED.etathT(i,:,1) = (1/NumDraws)*sum(TmpShock);
   if rank(HFED.R)>0;
      HFED.wthT(i,:) = (1/NumDraws)*sum(TmpwthT);
   end;
   %
   % sort the data
   %
   TmpEps = sort(TmpEps);
   drawnow;
   TmpKsi = sort(TmpKsi);
   drawnow;
   TmpShock = sort(TmpShock);
   drawnow;
   if rank(HFED.R)>0;
      TmpwthT = sort(TmpwthT);
      drawnow;
   end;
   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;
      HFED.Quantiles(l).epstth(i,:) = TmpEps(QuantVal,:);
      HFED.Quantiles(l).KsiErr(i,:) = TmpKsi(QuantVal,:);
      HFED.Quantiles(l).etathT(i,:,1) = TmpShock(QuantVal,:);
      if rank(HFED.R)>0;
         HFED.Quantiles(l).wthT(i,:) = TmpwthT(QuantVal,:);
      end;
   end;
   drawnow;
end;
%
% now we turn to the structural shocks
%
for i=1:n;
   for j=2:q;
      txt = ['Please wait while YADA computes the mean historical forecast error decompsotions as well as percentile values. 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;
      %
      % initialize
      %
      TmpShock = zeros(n,T);
      CurrSaves = 0;
      %
      % load the data
      %
      for ThisSave=1:NumFiles;
         file = [DSGEModel.OutputDirectory '\hfed\DSGE-HFED-' int2str(h) '-' DSGEModel.NameOfModel '-' num2str(CurrChain,'%0.0f') '-' num2str(ThisSave,'%0.0f') '-' num2str(NumFiles,'%0.0f') '.' num2str(NumDraws,'%0.0f') '-' num2str(TotalDraws,'%0.0f') '.mat'];
         HFEDData = load(file);
         FED = HFEDData.HFED;
         NumSaves = length(FED);
         for s=1:NumSaves;
            TmpShock(CurrSaves+s,:) = FED(s).HistDecompStructure.etathT(i,:,j);
         end;
         CurrSaves = CurrSaves+NumSaves;
         drawnow;
      end;
      %
      % compute the mean
      %
      HFED.etathT(i,:,j) = (1/NumDraws)*sum(TmpShock);
      %
      % sort the data
      %
      TmpShock = sort(TmpShock);
      drawnow;
      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;
         HFED.Quantiles(l).etathT(i,:,j) = TmpShock(QuantVal,:);
      end;
      drawnow;
   end;
end;
%
% close the wait dialog
%
set(WaitControls.text,'UserData','done');
delete(WaitHandle);
drawnow;
pause(0.02);

%
% end of HFEDDistribution.m
%
