function VDs = DSGEVarDecDistribution(DSGEModel,CurrINI,LevelsData,UseGroups,NumDraws,TotalDraws,CurrChain)
% DSGECVDDistribution: Computes the mean value of conditional variance decompositions as well as
%                      quantiles by loading data from disk.
%
% USAGE:
%
%       VDs = DSGEVarDecDistribution(DSGEModel,CurrINI,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.
%
%                  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: VDs (structure) with fields "FEVD", "LRVD", "Econ", "EconLR", and "Quantiles". The "FEVD" field holds a
%                                 3D matrix of dimension n x (n+q+1) 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
%                                 forecast error variance decomposition. The "LRVD" field holds a matrix of dimension nx(n+q)
%                                 with the posterior mean long-run forecast error variance decomposition. The "Econ" field holds
%                                 an nxh matrix with the posterior mean of forecast error variance decomposition for all
%                                 economic (structural) shocks. Similarly, "EconLR" is an nx1 vector with the posterior mean
%                                 long-run forecast error variance decomposition of all economic shocks. The field "Quantiles"
%                                 is a vector structure with fields "FEVD", "LRVD", "Econ", "EconLR", and "percent", where "FEVD"
%                                 holds a 3D matrix of dimension n x (n+q+1) x h with values from the percent quantile of the
%                                 marginal distribution of the forecast error variance decomposition, "LRVD" holds a matrix of
%                                 dimension nx(n+q) with values from the percent quantile of the marginal distribution of the
%                                 long-run forecast error variance decomposition, "Econ" is a matrix of dimension nxh with values
%                                 from the percent quantile of the marginal distribution of the forecast error variance
%                                 decomposition of all the economci shocks, "EconLR" is a matrix of dimension nx1 with values
%                                 from the percent quantile of the marginal distribution of the long-run forecast error variance
%                                 decomposition of all the economic shocks, while "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 23, 2008.
%                        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:
%
% * 25-04-2008: Changed shock to source on the wait dialog.
%
% * 23-05-2008: Updated the documentation.
%
% * 23-01-2009: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize
%
VDs = [];
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;
n = DSGEModel.n;
TypeStr = 'Original';
%
if LevelsData==1;
   TypeStr = 'Levels';
end;
%
if UseGroups==0;
   n1 = n;
   q = length(DSGEModel.StateShockPositions);
   ShockNames = '';
   for i=1:n;
      ShockNames = strvcat(ShockNames,['measurement error in ' StringTrim(DSGEModel.VariableNames(i,:)) ' equation']);
   end;
   ShockNames = strvcat(ShockNames,DSGEModel.StateShockNames,'signal extraction error');
else;
   IMatrix = SubIdentityMatrix(DSGEModel.ShockGroups);
   IMatrixSR = [ones(n,1) zeros(n,size(IMatrix,2)+1);zeros(size(IMatrix,1),1) IMatrix zeros(size(IMatrix,1),1);zeros(1,size(IMatrix,2)+1) 1];
   IMatrixLR = [ones(n,1) zeros(n,size(IMatrix,2));zeros(size(IMatrix,1),1) IMatrix];
   n1 = 1;
   q = max(DSGEModel.ShockGroups);
   ShockNames = char(zeros(size(DSGEModel.ShockGroupNames)));
   for i=1:q;
      ShockNames(i,:) = lower(strrep(DSGEModel.ShockGroupNames(i,:),'_',' '));
   end;
   ShockNames = strvcat('all measurement errors',ShockNames,'signal extraction error');
end;
r = length(DSGEModel.StateVariablePositions);
NumQuants = length(DSGEModel.Percentiles);
%
% initialize output
%
VDs.FEVD = zeros(n,n1+q+1,h);
VDs.LRVD = zeros(n,n1+q);
VDs.Econ = zeros(n,h);
VDs.EconLR = zeros(n,1);
for i=1:NumQuants;
   VDs.Quantiles(i).FEVD = zeros(n,n1+q+1,h);
   VDs.Quantiles(i).LRVD = zeros(n,n1+q);
   VDs.Quantiles(i).Econ = zeros(n,h);
   VDs.Quantiles(i).EconLR = zeros(n,1);
   VDs.Quantiles(i).percent = DSGEModel.Percentiles(i);
end;
%
% set up a wait dialog
%
txt = ['Please wait while YADA computes the posterior mean forecast error 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,source): (' StringTrim(DSGEModel.VariableNames(1,:)) ',' StringTrim(ShockNames(1,:)) ').'];
WaitHandle = WaitDLG(txt,'information',['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:n1+q+1;
      %
      % update the wait dialog
      %
      txt = ['Please wait while YADA computes the posterior mean forecast error 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,source): (' StringTrim(DSGEModel.VariableNames(i,:)) ',' StringTrim(ShockNames(j,:)) ').'];
      set(WaitControls.text,'String',txt);
      drawnow;
      %
      TempVar = zeros(NumDraws,h);
      if j<n1+q+1;
         TempLR = zeros(NumDraws,1);
      else;
         TempLR = [];
      end;
      if j==1;
         TempEcon = zeros(NumDraws,h);
         TempEconLR = zeros(NumDraws,1);
      end;
      CurrSaves = 0;
      for ThisSave=1:NumFiles;
         file = [DSGEModel.OutputDirectory '\vardec\DSGE-VD-' 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'];
         VDData = load(file);
         VD = VDData.VarDec;
         NumSaves = length(VD);
         %
         % 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 = VD(s).FEVD;
               LRVD = VD(s).LRVD;
            else;
               FEVD = zeros(n,n1+q+1,h);
               for k=1:h;
                  FEVD(:,:,k) = VD(s).FEVD(:,:,k)*IMatrixSR;
               end;
               LRVD = VD(s).LRVD*IMatrixLR;
            end;
            TempVar(CurrSaves+s,:) = FEVD(i,j,:);
            if j<n1+q+1;
               TempLR(CurrSaves+s,1) = LRVD(i,j);
            end;
            if j==1;
               for j1=1:h;
                  TempEcon(CurrSaves+s,j1) = sum(FEVD(i,n1+1:n1+q,j1),2);
               end;
               TempEconLR(CurrSaves+s,1) = sum(LRVD(i,n1+1:n1+q),2);
            end;
         end;
         CurrSaves = CurrSaves+NumSaves;
         drawnow;
      end;
      %
      % now we sort and compute the mean the data
      %
      if j<n1+q+1;
         TempLR = sort(TempLR);
         if max(abs(TempLR))>1.0e-10;
            VDs.LRVD(i,j) = mean(TempLR);
         end;
         if j==1;
            TempEconLR = sort(TempEconLR);
            if max(abs(TempEconLR))>1.0e-10;
               VDs.EconLR(i,1) = mean(TempEconLR);
            end;
         end;
         %
         % deal with quantiles
         %
         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(TempLR(QuantVal,1))>1.0e-10;
               VDs.Quantiles(l).LRVD(i,j) = TempLR(QuantVal,1);
            end;
            if j==1;
               if abs(TempEconLR(QuantVal,1))>1.0e-10;
                  VDs.Quantiles(l).EconLR(i,1) = TempEconLR(QuantVal,1);
               end;
            end;
         end;
      end;
      %
      for k=1:h;
         TempVar(:,k) = sort(TempVar(:,k));
         if max(abs(TempVar(:,k)))>1.0e-10;
            VDs.FEVD(i,j,k) = mean(TempVar(:,k));
         end;
         if j==1;
            TempEcon(:,k) = sort(TempEcon(:,k));
            if max(abs(TempEcon(:,k)))>1.0e-10;
               VDs.Econ(i,k) = mean(TempEcon(:,k));
            end;
         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;
               VDs.Quantiles(l).FEVD(i,j,k) = TempVar(QuantVal,k);
            end;
            if j==1;
               if abs(TempEcon(QuantVal,k))>1.0e-10;
                  VDs.Quantiles(l).Econ(i,k) = TempEcon(QuantVal,k);
               end;
            end;
         end;
         drawnow;
      end;
   end;
end;
%
% close the wait dialog
%
set(WaitControls.text,'UserData','done');
delete(WaitHandle);
drawnow;
pause(0.02);

%
% end of DSGEVarDecDistribution.m
%
