function [ForecastResults,status,kalmanstatus] = DSGEMultistepForecast(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI)
% DSGEMultistepForecast: Computes multistep forecasts and MSE matrices.
%
% USAGE:
%
%       [ForecastResults,status,kalmanstatus] = DSGEMultistepForecast(theta,thetaPositions,ModelParameters,DSGEModel,CurrINI)
%
% 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.
%
% REQUIRED OUTPUT: ForecastResults (structure) of dimension equal to the frequency of parameter re-estimation
%                                 when forecasting or to the number of times we can compute 1-step ahead forecast
%                                 errors (whichever is the smallest). The fields in this vector structure are
%                                 given by the matrices "Yhat", the forecasts Y(T+h|T), "Y", the observed
%                                 variables Y(T+h), "Error", the forecast errors, "RMSE", the theoretical root
%                                 mean squared errors for Y(T+h), "AnnYhat", the annualized forecasts, "AnnError"
%                                 the annualized forecast errors, and "DateVector" containing the dates. The
%                                 number of columns of the matrices and vectors are equal to the maximum possible
%                                 desired forecast 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.
%                                      mocde = 8:      The function "compute_aim_matrices" returns complex
%                                                      numbers.
%                                      else            Return code not properly specified
%
%                  kalmanstatus (boolean) indicating the covariance stationarity of the state variables. If all the
%                                 eigenvalues of F are inside the unit circle, then kalmanstatus=0 and the log-likelihood function
%                                 is calculated. If not, then status=1 and initP is set to 0.
%
% NOTE: Three 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: January 5, 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:
%
% * 08-01-2007: Finished the function.
%
% * 15-01-2007: Updated the documentation.
%
% * 22-01-2007: Added the "AnnY", "AnnYhat", and "AnnError" entries.
%
% * 25-01-2007: Made sure that the "Annualize" function receives the relevant terms of the DSGEModel.annualscale
%               vector.
%
% * 09-03-2007: Missing "CurrINI" input argument was added.
%
% * 20-04-2007: The fields "YAll", "AnnYAll", and "DateVectorAll" have been added to the output structure.
%
% * 13-11-2007: Updated the documentation.
%
% * 16-04-2008: The code now checks if the model has a unique convergent solution also at the posterior mode
%               estimates.
%
% * 16-05-2008: 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.
%
% * 02-07-2008: Added the variable "DSGEModel.AllowUnitRoot" for the Kalman filters.
%
% * 03-07-2008: Made sure that the measurement equation receives only the active state variable names.
%               Added the try-catch code for AiMSolver.
%
% * 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
%
ForecastResults = [];
%
% 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;
else;
   status = 1;
   kalmanstatus = 0;
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;
   if nargout>1;
      status = 0;
      kalmanstatus = 0;
   end;
   About(ErrorStr,'information','YADA - Bad MAT File',120,500,CurrINI);
   return;
end;
%
% 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
   %
   ErrorStr = ['Error found in "' GetName(DSGEModel.MeasurementEquationFile) '.m". Message caught is: ' lasterr];
   About(ErrorStr,'error','Error - Measurement Equation',150,500,CurrINI);
   return;
end;
%
% set the InitialStateVector value
%
if DSGEModel.UseOwnInitialState==1;
   if length(DSGEModel.InitialStateValues)==size(AIMData.endog,1);
      DSGEModel.InitialStateVector = DSGEModel.InitialStateValues(DSGEModel.StateVariablePositions);
   else;
      DSGEModel.InitialStateVector = zeros(length(DSGEModel.StateVariablePositions),1);
   end;
else;
   DSGEModel.InitialStateVector = zeros(length(DSGEModel.StateVariablePositions),1);
end;
%
% Compute the max length of forecasts based on the starting position
%
MaxFH = CalcMaxForecastHorizon(DSGEModel.MaxForecastHorizon,DSGEModel.DataFrequency,DSGEModel.EndYear,DSGEModel.EndPeriod,DSGEModel.SubEndYear,DSGEModel.SubEndPeriod);
%
% Get the sub-sample
%
[FirstPeriod,LastPeriod] = CreateSubSample(DSGEModel);
%
% set up a vector with date value for the full data sample
%
DateVector = MakeDateVector(DSGEModel.BeginYear,DSGEModel.BeginPeriod,DSGEModel.EndYear,DSGEModel.EndPeriod,DSGEModel.DataFrequency);
%
% the loop over the DSGEModel.FreqParameterReest parameter should begin here (number of times to use re-estimated parameters 
%
for i=1:DSGEModel.FreqParameterReest;
   if LastPeriod+i<=size(DSGEModel.Y,2);
      %
      % We can now run the Kalman filter
      %
      [logLikeValue,stat,lnLt,Yhat,MSEY,Ksitt1,Ptt1] = KalmanFilter(DSGEModel.Y(:,FirstPeriod:LastPeriod+i),DSGEModel.X(:,FirstPeriod:LastPeriod+i),A,H,F,B0*B0',R,DSGEModel.InitialStateVector,DSGEModel.UseDoublingAlgorithm,DSGEModel.DAMaximumIterationsValue,DSGEModel.DAToleranceValue,DSGEModel.KalmanFirstObservation,DSGEModel.AllowUnitRoot,DSGEModel.StateCovConstValue);
      %
      % we should adjust everything based on DSGEModel.KalmanFirstObservation
      %
      if nargout>1;
         kalmanstatus = stat;
         if kalmanstatus~=0;
            return;
         end;
      end;
      %
      % compute the possible forecast horizon
      %
      PosFH = min(size(DSGEModel.Y,2)-LastPeriod+1-i,MaxFH);
      Y = zeros(DSGEModel.n,PosFH);
      RMSE = Y;
      Error = Y;
      Y(:,1) = Yhat(:,size(Yhat,2));
      Error(:,1) = DSGEModel.Y(:,LastPeriod+i)-Y(:,1);
      PttEst = Ptt1(:,:,size(Yhat,2));
      RMSE(:,1) = sqrt(diag((H'*PttEst*H)+R));
      Fi = F;
      for h=2:PosFH;
         %
         % compute y(T+h|T)
         %
         Y(:,h) = (A'*DSGEModel.X(:,LastPeriod+i-1+h))+(H'*Fi*Ksitt1(:,size(Ksitt1,2)));
         Error(:,h) = DSGEModel.Y(:,LastPeriod+i-1+h)-Y(:,h);
         PttEst = (F*PttEst*F')+(B0*B0');
         RMSE(:,h) = sqrt(diag((H'*PttEst*H)+R));
         Fi = Fi*F;
      end;
      ForecastResults(i).Yhat = Y;
      ForecastResults(i).Y = DSGEModel.Y(:,LastPeriod+i:LastPeriod+i-1+PosFH);
      ForecastResults(i).YAll = DSGEModel.Y(:,FirstPeriod:LastPeriod+i-1+PosFH);
      ForecastResults(i).Error = Error;
      ForecastResults(i).RMSE = RMSE;
      ForecastResults(i).DateVector = DateVector(LastPeriod+i:LastPeriod+i-1+PosFH);
      ForecastResults(i).DateVectorAll = DateVector(FirstPeriod:LastPeriod+i-1+PosFH);
      if max(DSGEModel.annual)>1;
         ForecastResults(i).AnnYhat = Annualize([DSGEModel.Y(:,LastPeriod+i+1-max(DSGEModel.annual):LastPeriod+i-1) Y],DSGEModel.annual,DSGEModel.annualscale);
         ForecastResults(i).AnnError = Annualize([zeros(DSGEModel.n,max(DSGEModel.annual)-1) Error],DSGEModel.annual,DSGEModel.annualscale);
         ForecastResults(i).AnnY =  Annualize([DSGEModel.Y(:,LastPeriod+i+1-max(DSGEModel.annual):LastPeriod+i-1+PosFH)],DSGEModel.annual,DSGEModel.annualscale);
         ForecastResults(i).AnnYAll = Annualize(DSGEModel.Y(:,max(FirstPeriod+1-max(DSGEModel.annual),1):LastPeriod+i-1+PosFH),DSGEModel.annual,DSGEModel.annualscale);
      else;
         ForecastResults(i).AnnYhat = Y;
         ForecastResults(i).AnnError = Error;
         ForecastResults(i).AnnY = ForecastResults(i).Y
         ForecastResults(i).AnnYAll = ForecastResults(i).YAll;
      end;
   else;
      break;
   end;
end;

%
% end of DSGEMultistepForecast.m
%
