function [IRStructure,status] = CalculateDSGEImpulseResponse(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI,UseLevels)
% CalculateDSGEImpulseResponse: Computes the impulse responses for observed variables and state variables for a given
%                               value of the parameters.
%
% USAGE:
%
%       [IRStructure,status] = CalculateDSGEImpulseResponse(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI,UseLevels)
%
% REQUIRED INPUT:  theta (vector) with values 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.
%
%                  UseLevels (interger) which indicates if variables in first differences should be
%                                 accumulated (1), annualized (2), or left alone (0).
%
% REQUIRED OUTPUT: IRStructure (structure) with fields Ksi (rxqxh+1) and Y (nxqxh+1). The responses in the invidual
%                                 variables are located in the rows, while the columns give the shocks. The third
%                                 dimension is the response horizon.
%
% OPTIONAL OUTPUT: status (integer) Indicates whether the solution of the model is unique or not.
%                                   mcode is 1 when unique, 0 if not computed, and various other values when
%                                   there are problems with the solution. Specifically:
%
%                                      mcode = 2:      Roots not correctly computed by real_schur
%                                      mcode = 3:      Too many big roots
%                                      mcode = 35:     Too many big roots, and q(:,right) is singular
%                                      mcode = 4:      Too few big roots
%                                      mcode = 45:     Too few big roots, and q(:,right) is singular
%                                      mcode = 5:      q(:,right) is singular
%                                      mcode = 61:     Too many exact shiftrights
%                                      mcode = 62:     Too many numeric shiftrights
%                                      mcode = 7:      The a matrix has infinite or NaN values. Eigenvalues
%                                                      cannot be calculated.
%                                      mcode = 8:      The function "compute_aim_matrices" returns complex
%                                                      numbers.
%                                      else            Return code not properly specified
%
% NOTE: Two output arguments should only be used when theta has the initial values.
%
%
%                       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: August 15, 2006.
%                        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:
%
% * 17-08-2006: Added the two new inputs for the measurement equation file with string matrices containing the names of
%               the observed variables and of the exogenous variables.
%
% * 21-08-2006: Updated the documentation.
%
% * 25-08-2006: The compute_aim_matrices, measurement equation file, and the file with parameters to update are now
%               always copied to the tmp directory to make sure that the correct version is located there.
%
% * 14-09-2006: Added the changes to the measurement equation call.
%
% * 17-10-2006: Added a try, catch call when attempting to load mat-files. This is needed if matlab can't
%               read the mat-file, e.g., because it was create with a later version of matlab.
%
% * 07-12-2006: Took the change from CurrINI.IRHorizon to DSGEModel.IRHorizon into account.
%
% * 11-12-2006: The input variable "UseLevels" is now accepted.
%
% * 18-12-2006: Save the results in IRStructure to file.
%
% * 02-01-2007: Changed the Query figure name for saving to file.
%
% * 03-01-2007: When saving data to disk, the function now writes a matlab file to disk for reading the data.
%
% * 23-01-2007: The "UseLevels" variable can now take on the value 2, implying that annualized impulse responses
%               are calculated.
%
% * 25-01-2007: Made sure that the "Annualize" function receives the relevant terms of the DSGEModel.annualscale
%               vector.
%
% * 26-02-2007: Changed the output into a structure with fields "Ksi" and "Y", both being 3D matrices.
%
% * 13-11-2007: Updated the documentation.
%
% * 07-01-2008: Added drawnow calls after Query and About calls.
%
% * 16-04-2008: The code now checks if the model has a unique convergent solution also at the posterior mode
%               estimates.
%
% * 16-05-2088: Took the "DSGEModel.AIMTolerance" input for AiMSolver into account.
%
% * 23-05-2008: Updated the documentation.
%
% * 27-05-2008: The mcode value 7 is now supported.
%
% * 03-07-2008: Made sure that the measurement equation receives only the active state variable names.
%               Added try-catch code for AiMSolver.
%
% * 09-07-2008: Added code to handle time-varying measurement matrix H.
%
% * 22-01-2009: Updated the documentation.
%
% * 30-01-2009: Made sure that status is only set to mcode when the latter exists.
%
% * 09-02-2009: Changed the input "length(DSGEModel.StateVariablePositions)" to "size(F,1)" when calling the
%               measurement equation function.
%
% * 13-02-2009: The mcode values 7 and 8 are now properly mentioned when encountered.
%
% * 09-03-2009: Added quotation marks around paths for the dos copy commands.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

%
% initialize the output
%
IRStructure = [];
%
% set up the ModelParameters structure based on theta and thetaPositions
%
if nargout==1;
   %
   % Make sure that the ModelParameter structure takes the correct values.
   % This is not necessary to do when initial values are given to the function.
   % In that case, more substantial checking is required!
   %
   ModelParameters = ThetaToModelParameters(ModelParameters,theta,thetaPositions);
   %
   % 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;
   %
   % Update any parameteters that need updating
   %
   if FileExist(DSGEModel.UpdateParameterFile)==1;
      eval(['ModelParameters = ' GetName(DSGEModel.UpdateParameterFile) '(ModelParameters);']);
   end;
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)==1;
   %
   % See if we can solve the model
   %
   if nargout==1;
      [mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
   else;
      try;
         [mcode,MaxAbsError,BMatrix,SMatrix] = AiMSolver(ModelParameters,AIMData.neq,AIMData.nlag,AIMData.nlead,DSGEModel.AIMTolerance);
      catch;
         if exist('mcode')>0;
            status = mcode;
         else;
            status = 0;
         end;
         ErrorStr = ['There was an error running AiM. YADA received the following message: ' lasterr];
         About(ErrorStr,'error','AiM Error',300,500,CurrINI);
         drawnow;
         return;
      end;
   end;
   if nargout==2;
      status = mcode;
      if status~=1;
         return;
      end;
   else;
      if mcode~=1;
         if mcode==2;
            mcodeStr = 'Roots are not correctly computed by real_schur.'
         elseif mcode==3;
            mcodeStr = 'Too many big roots.';
         elseif mcode==35;
            mcodeStr = 'Too many big roots, and q(:,right) is singular.';
         elseif mcode==4;
            mcodeStr = 'Too few big roots.';
         elseif mcode==45;
            mcodeStr = 'Too few big roots, and q(:,right) is singular.';
         elseif mcode==5;
            mcodeStr = 'q(:,right) is singular.';
         elseif mcode==61;
            mcodeStr = 'Too many exact shiftrights.';
         elseif mcode==62;
            mcodeStr = 'Too many numeric shiftrights.';
         elseif mcode==7;
            mcodeStr = 'Infinite or NaN values detected.';
         elseif mcode==8;
            mcodeStr = 'The function "compute_aim_matrices" returns complex numbers.';
         else;
            mcodeStr = 'Return code not properly specified.';
         end;
         txt = ['The AiM solver provided the return code: ' num2str(mcode) ', i.e., "' mcodeStr '"'];
         About(txt,'information','YADA - AiM Solver Error',200,500,CurrINI);
         return;
      end;
   end;
   %
   % Convert to the state equation matrices
   %
   [F,B0] = AiMtoStateSpace(BMatrix,SMatrix,DSGEModel.StateVariablePositions,DSGEModel.StateShockPositions);
   %
   % obtain the measurement equation matrices
   %
   try;
      eval(['[A,H,R] = ' GetName(DSGEModel.MeasurementEquationFile) '(ModelParameters,AIMData.endog(DSGEModel.StateVariablePositions,:),DSGEModel.VariableNames,DSGEModel.XVariableNames,DSGEModel.n,size(F,1),DSGEModel.k);']);
   catch;
      %
      % display an error message
      %
      if nargout==2;
         status = 0;
      end;
      ErrorStr = ['Error found in "' GetName(DSGEModel.MeasurementEquationFile) '.m". Message caught is: ' lasterr];
      About(ErrorStr,'error','Error - Measurement Equation',150,500,CurrINI);
      drawnow;
      return;
   end;
   if length(size(H))==3;
      %
      % check if H has the required size
      %
      [FirstPeriod,LastPeriod] = CreateSubSample(DSGEModel);
      if size(H,3)>=LastPeriod;
         if nargout==1;
            TypeStr = 'Posterior Mode';
         else;
            TypeStr = 'Initial Values';
         end;
         %
         % 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],'','',CurrINI);
         %
         SelPos = FirstPeriod+DSGEModel.KalmanFirstObservation-1+(SelectedPosition-1);
         if strcmp(lower(action),'cancel')==1;
            return;
         end;
         H = H(:,:,SelPos);
      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;
   %
   % 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;
   IRHorizon = SampleFreq*DSGEModel.IRHorizon;
   %
   % Now we can compute the impulse responses
   %
   if UseLevels==0;
      IRStructure = DSGEImpulseResponseFcn(H,F,B0,IRHorizon);
   elseif UseLevels==1;
      IRStructure = DSGELevImpulseResponseFcn(H,F,B0,diag(1-DSGEModel.levels),IRHorizon);
   else;
      IRStructure = DSGEImpulseResponseFcn(H,F,B0,IRHorizon);
      if isempty(IRStructure)==0;
         IRFY = zeros(size(H,2),IRHorizon+1,size(B0,2));
         for i=1:IRHorizon+1;
            for j=1:size(B0,2);
               IRFY(:,i,j) = IRStructure.Y(:,j,i);
            end;
         end;
         %
         % we can now annualize
         %
         for j=1:size(B0,2);
            IRFY(:,:,j) = Annualize([zeros(DSGEModel.n,max(DSGEModel.annual)-1) IRFY(:,:,j)],DSGEModel.annual,DSGEModel.annualscale);
         end;
         %
         % now we save this to the structure IRStrcuture again
         %
         for i=1:IRHorizon+1;
            for j=1:size(B0,2);
               IRStructure.Y(:,j,i) = IRFY(:,i,j);
            end;
         end;
      end;
   end;
   if isempty(IRStructure)==0;
      txt = 'Would you like to save the impulse response functions to file?';
      answer = Query(txt,'question',140,'Save - Impulse Responses',500,'no',CurrINI);
      drawnow;
      if strcmp(lower(answer),'yes')==1;
         NumResp = length(IRStructure);
         IRFState = zeros(size(F,1),size(B0,2),NumResp);
         IRFY = zeros(size(H,2),size(B0,2),NumResp);
         YNames = DSGEModel.VariableNames;
         StateVariableNames = DSGEModel.StateVariableNames;
         ShockNames = DSGEModel.StateShockNames;
         IRFState = IRStructure.Ksi;
         IRFY = IRStructure.Y;
         file = [DSGEModel.OutputDirectory '\ImpulseResponses-'];
         if nargout==1;
            file = [file 'PosteriorMode-'];
         else;
            file = [file 'InitialValues-'];
         end;
         if UseLevels==0;
            file = [file 'Original.mat'];
         elseif UseLevels==1;
            file = [file 'Levels.mat'];
         else;
            file = [file 'Annual.mat'];
         end;
         save(file,'YNames','StateVariableNames','ShockNames','IRFState','IRFY');
         %
         % create a matlab file to read the data
         %
         mfile = strrep(strrep(strrep(file,'.mat','.m'),'-',''),'ImpulseResponses','IR');
         fid = fopen(mfile,'wt');
         fprintf(fid,['%%\n%% load the data in ' GetFilename(file) '\n%%\n']);
         fprintf(fid,'DataStructure = load(''%s'');\n',GetFilename(file));
         fprintf(fid,'%%\n%% string matrix with the names of the observed variables\n%%\n');
         fprintf(fid,'YNames = DataStructure.YNames;\n');
         fprintf(fid,'%%\n%% string matrix with the names of the state variables\n%%\n');
         fprintf(fid,'StateVariableNames = DataStructure.StateVariableNames;\n');
         fprintf(fid,'%%\n%% string matrix with the names of the economic shocks\n%%\n');
         fprintf(fid,'ShockNames = DataStructure.ShockNames;\n');
         fprintf(fid,'%%\n%% a 3D matrix of dimension rxqxh with the responses for the r state variables from the q shocks from period 0 to period h-1\n%%\n');
         fprintf(fid,'IRFState = DataStructure.IRFState;\n');
         fprintf(fid,'%%\n%% a 3D matrix of dimension nxqxh with the responses for the n observed variables\n%%\n');
         fprintf(fid,'IRFY = DataStructure.IRFY;\n');
         fprintf(fid,'%%\n%% Add your own commands below\n%%\n\n\n');
         %
         % end the script file
         %
         fprintf(fid,'%%\n%% Created by YADA on %s\n%%\n',datestr(now,0));
         fclose(fid);
         txt = ['The impulse response data have been saved to the file "' GetFilename(file) '" in the directory "' GetPath(file) '". The file contains 5 entries: YNames (string matrix with the names of the observed variables), StateVariableNames (string matrix with the names of the state variables), ShockNames (string matrix with the names of the economic shocks), IRFState (a 3D matrix of dimension rxqxh with the responses for the r state variables from the q shocks from period 0 to period h-1), and IRFY (a 3D matrix of dimension nxqxh with the responses for the n observed variables). To access this data you may run the matlab script file "' GetFilename(mfile) '".'];
         About(txt,'information','Impulse Responses',200,500,CurrINI);
         drawnow;
      end;
   end;
else;
   if nargout==2;
      status = 0;
   end;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   drawnow;
end;

%
% end of CalculateDSGEImpulseResponse.m
%
