function DoneCalc = DSGEImpulseResponses(thetaPostSample,thetaPositions,ModelParameters,DSGEModel,CurrINI,AnnualizeData,LevelsData,TotalDraws,CurrChain)
% DSGEImpulseResponses: Calculates the impulse responses for the posterior draws of the DSGE model and stores them on disk.
%
% USAGE:
%
%       DoneCalc = DSGEImpulseResponses(thetaPostSample,thetaPositions,ModelParameters,DSGEModel,CurrINI,AnnualizeData,LevelsData,TotalDraws,CurrChain)
%
% REQUIRED INPUT:  thetaPostSample (matrix) of dimension NumDrawsxNumParam with draws from the posterior distribution
%                                  of the original DSGE model parameters.
%
%                  thetaPositions (structure) with fields parameter and distribution. This function
%                                 needs the former field to know which parameter name corresponds to
%                                 position i in theta.
%
%                  ModelParameters (structure) whose fields have names given by the parameter names in the
%                                 DSGE model.
%
%                  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.
%
%                  TotalDraws (integer) with the total number of parameter draws that could have been used.
%
%                  CurrChain (integer) with the number for the current Markov chain.
%
% 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: February 27, 2007.
%                        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:
%
% * 04-04-2007: Added the "TotalDraws" input variable.
%
% * 13-11-2007: Updated the documentation.
%
% * 21-11-2007: Fixed a spelling error.
%
% * 10-01-2008: Changed the last drawnow call to pause(0.02).
%
% * 23-05-2008: Updated the documentation.
%
% * 09-07-2008: Added code to handle a time-varying measurement matrix H.
%
% * 15-07-2008: Fixed a bug where "SelPos" was not applied to a 3D H-matrix.
%
% * 22-01-2009: Updated the documentation.
%
% * 09-03-2009: Added quotation marks around paths for the dos copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize
%
NumDraws = size(thetaPostSample,1);
DoneCalc = 0;
%
% 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,'irfs');
if stat~=1;
   txt = ['YADA was for some reason unable to create the directory "' DSGEModel.OutputDirectory '\irfs". The computation of impulse response distributions has therefore been aborted.'];
   About(txt,'information','YADA - Directory Creation Problem',180,500,CurrINI);
   return;
end;
%
% 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';
if AnnualizeData==1;
   InitIR = zeros(DSGEModel.n,max(DSGEModel.annual)-1);
   TypeStr = 'Annualized';
end;
if LevelsData==1;
   AccMat = diag(1-DSGEModel.levels);
   TypeStr = 'Levels';
end;
%
% initialize index for saving
%
SaveAfterDraws = min(NumDraws,1000);
NumFiles = ceil(NumDraws/SaveAfterDraws);
q = length(DSGEModel.StateShockPositions);
n = DSGEModel.n;
r = length(DSGEModel.StateVariablePositions);
%
% check if we have a time-varying measurement matrix H
%
[A,H,R,F,B0] = SolveDSGEModel(DSGEModel,thetaPostSample(1,:)',thetaPositions,ModelParameters,AIMData);
if length(size(H))==3;
   %
   % check if H has the required size
   %
   [FirstPeriod,LastPeriod] = CreateSubSample(DSGEModel);
   if size(H,3)>=LastPeriod;
      %
      % pick a certain period for the H matrix
      %
      YADAHandle = findobj('Type','figure','Tag','YADA');
      controls = get(YADAHandle,'UserData');
      SubSampleMat = get(controls.dsge.initializekalman,'String');
      SubSampleMat = SubSampleMat(DSGEModel.KalmanFirstObservation:size(SubSampleMat,1),:);
      [action,SelectedPosition] = SelectionDlg(SubSampleMat,size(SubSampleMat,1),'Select time period for the measurement matrix H','Measurement Matrix',['Impulse Responses - ' TypeStr ' Data'],'','',CurrINI);
      %
      SelPos = FirstPeriod+DSGEModel.KalmanFirstObservation-1+(SelectedPosition-1);
      if strcmp(lower(action),'cancel')==1;
         return;
      end;
   else;
      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(LastPeriod) '). YADA has therefore aborted from computing the impulse responses.'];
      About(ErrorStr,'error','Error - Measurement Equation',200,500,CurrINI);
      return;
   end;
end;
%
% begin the loop
%
abort = '';
MeanSolutionTime = 0;
LastDraw = 0;
ThisSave = 0;
IR = [];
it = 1;
%
% 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 DSGE model solution time:';
   %
   ProgressStructure.name = ['Impulse Responses - ' TypeStr ' Data'];
   ProgressStructure.CurrINI = CurrINI;
   WaitHandle = ProgressDLG(0,ProgressStructure);
   set(WaitHandle,'Color',get(0,'defaultuicontrolbackgroundcolor'));
   drawnow;
else;
   txt = ['Please wait while YADA computes impulse responses for ' lower(TypeStr) ' data. Computations started at: ' StringTrim(datestr(now,14))];
   WaitHandle = WaitDLG(txt,'information',['Impulse Responses - ' TypeStr ' Data'],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 MeanSolutionTime]);
      end;
   end;
   %
   % Solve the model for the current parameter vector
   %
   theta = thetaPostSample(it,:)';
   tic;
   [A,H,R,F,B0] = SolveDSGEModel(DSGEModel,theta,thetaPositions,ModelParameters,AIMData);
   MeanSolutionTime = (((it-1)/it)*MeanSolutionTime)+((1/it)*toc);
   IR(LastDraw).Y = zeros(n,q,h+1);
   if length(size(H))==3;
      H = H(:,:,SelPos);
   end;
   %
   % replace small numbers with zero
   %
   Bpos = (B0>1.0e-12);
   Bneg = (B0<1.0e-12);
   B0 = (B0 .* Bpos)+(B0 .* Bneg);
   if LevelsData==0;
      %
      % compute the basic impulse responses
      %
      IR(LastDraw).Y(:,:,1) = H'*B0;
      Fi = eye(r);
      for i=1:h;
         Fi = F*Fi;
         Y = H'*Fi*B0;
         Ypos = (Y>1.0e-12);
         Yneg = (Y<1.0e-12);
         IR(LastDraw).Y(:,:,i+1) = (Y .* Ypos) + (Y .* Yneg);
      end;
      %
      % check if we should annualize
      %
      if AnnualizeData==1;
         %
         % reorganize the 3D output
         %
         IRFY = zeros(n,h+1,q);
         for i=1:h+1;
            for j=1:q;
               IRFY(:,i,j) = IR(LastDraw).Y(:,j,i);
            end;
         end;
         %
         % now annualize per shock
         %
         for j=1:q;
            IRFY(:,:,j) = Annualize([InitIR IRFY(:,:,j)],DSGEModel.annual,DSGEModel.annualscale);
         end;
         %
         % reverse the 3D output ordering to the original
         %
         for i=1:h+1;
            for j=1:q;
               IR(LastDraw).Y(:,j,i) = IRFY(:,i,j);
            end;
         end;
      end;
   else;
      %
      % we have impulse responses for levels
      %
      IR(LastDraw).Y(:,:,1) = H'*B0;
      Fi = eye(r);
      for i=1:h;
         AddY = AccMat*IR(LastDraw).Y(:,:,i);
         Fi = F*Fi;
         Y = (H'*Fi*B0)+AddY;
         Ypos = (Y>1.0e-12);
         Yneg = (Y<1.0e-12);
         IR(LastDraw).Y(:,:,i+1) = (Y .* Ypos) + (Y .* Yneg);
      end;
   end;
   %
   % 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 '\irfs\DSGE-IR-' 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'];;
      save(file,'IR');
      %
      % restore the data structure
      %
      IR = [];
   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 DSGEImpulseResponses.m
%
