function ForecastResults = BVARMultistepForecast(Psi,Pi,Omega,DSGEModel)
% BVARMultistepForecast: Computes multistep ahead forecasts from the BVAR model with fixed parameter
%                        values.
%
% USAGE:
%
%       ForecastResults = BVARMultistepForecast(Psi,Pi,Omega,DSGEModel)
%
% REQUIRED INPUT:  Psi (pxq) matrix with steady state parameters.
%
%                  Pi (pxpk) matrix with the parameters on lagged enogenous variables.
%
%                  Omega (pxp) covariance matrix of the VAR residuals.
%
%                  DSGEModel (structure) with fields containing all the information about a DSGE Model for YADA.
%
% 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), and "DateVector" containing the dates. The
%                                  number of columns of the matrices and vectors are equal to the maximum possible
%                                  desired forecast horizon. Next, the field "rwY" contains the forecasts based on
%                                  the assumption that Y is a random walk, while the field "rwError" holds the
%                                  forecast errors based on these forecasts.
%
%
%                       Written by: Anders Warne
%                                   New Area Wide Model Project
%                                   DG-R/EMO
%                                   European Central Bank (ECB)
%                                   Email: anders.warne@ecb.europa.eu
%                                   Copyright © 2006-2008 European Central Bank.
%
%                       First version: January 15, 2007.
%                        This version: December 19, 2008.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 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:
%
% * 16-1-2007:  Added the random walk forecasts.
%
% * 22-1-2007:  Added the entries "AnnY", "AnnYhat", "AnnError", "AnnrwY", and "AnnrwError".
%
% * 23-1-2007:  Fixed a bug where "DSGEModel.BVARY" was missing in two places when picking entries
%               from "DSGEModel.Y", thus resulting in concatenation errors.
%
% * 24-1-2007:  Fixed bugs related to the using the vector DSGEModel.annual when the BVAR has fewer
%               variables than the DSGE model.
%
% * 25-1-2007:  Made sure that the "Annualize" function receives the relevant terms of the DSGEModel.annualscale
%               vector.
%
% * 8-5-2007:   Updated the documentation.
%
% * 13-11-2007: Updated the documentation.
%
% * 2-1-2008:   Removed code using "DSGEModel.InitializeKalman".
%
% * 23-5-2008:  Updated the documentation.
%
% * 19-12-2008: Updated the documentation.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

ForecastResults = [];
%
% set up the data matrices
%
p = length(DSGEModel.BVARY);
k = DSGEModel.BVARLags;
q = length(DSGEModel.BVARX);
%
% Determine the sample
%
[FirstPeriod,LastPeriod] = CreateSubSample(DSGEModel);
%
% Determine the first observation to use when computing log-likelihood with the Kalman filter
%
KalmanFirstObservation = DSGEModel.KalmanFirstObservation;
%
% Check if we should adjust this based on k, the lag order
%
if FirstPeriod+KalmanFirstObservation-1<k+1;
   KalmanFirstObservation = k+2-FirstPeriod;
end;
%
% rewrite the VAR on VAR form
%
if k>1;
   B = [Pi;eye(p*(k-1)) zeros(p*(k-1),p)];
   J = [eye(p) zeros(p,p*(k-1))];
else;
   B = Pi;
   J = eye(p);
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);
%
% 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);
      %
      % Setup all the basic data matrices
      %
      x = DSGEModel.Y(DSGEModel.BVARY,FirstPeriod+KalmanFirstObservation-1:LastPeriod+i-1);
      T = size(x,2);
      X = zeros(p*k,T);
      for j=1:k;
         X(1+(p*(j-1)):p*j,:) = DSGEModel.Y(DSGEModel.BVARY,FirstPeriod+KalmanFirstObservation-1-j:LastPeriod-j+i-1);
      end;
      if isempty(DSGEModel.BVARX)==0;
         d = DSGEModel.X(DSGEModel.BVARX,FirstPeriod+KalmanFirstObservation-1:LastPeriod+i-1);
         dLag = zeros(q*k,T);
         for j=1:k;
            dLag(1+(q*(j-1)):q*j,:) = -DSGEModel.X(DSGEModel.BVARX,FirstPeriod+KalmanFirstObservation-1-j:LastPeriod-j+i-1);
         end;
      else;
         d = [];
         dLag = [];
      end;
      %
      % Set up data matrices for computing the log-likelihood
      %
      if q>0;
         y = x-(Psi*d);
         Y = X+(kron(eye(k),Psi)*dLag);
      else;
         y = x;
         Y = X;
      end;

      %
      % compute the possible forecast horizon
      %
      PosFH = min(size(DSGEModel.Y,2)-LastPeriod+1-i,MaxFH);
      Ymat = zeros(p,PosFH);
      RMSE = Ymat;
      Bh = eye(size(B,1));
      for h=1:PosFH;
         %
         % compute a forecast of y, then add the deterministic term later
         %
         Bh = B*Bh;
         Ymat(:,h) = J*Bh*Y(:,T);
         %
         % compute the MSE
         %
         MSE = Omega;
         for i1=1:h-1;
            Bi1 = J*(B^i1)*J';
            MSE = MSE+(Bi1*Omega*Bi1');
         end;
         RMSE(:,h) = sqrt(abs(diag(MSE)));
      end;
      ForecastResults(i).Yhat = Ymat+(Psi*DSGEModel.X(DSGEModel.BVARX,LastPeriod+i:LastPeriod+i-1+PosFH));
      ForecastResults(i).Y = DSGEModel.Y(DSGEModel.BVARY,LastPeriod+i:LastPeriod+i-1+PosFH);
      ForecastResults(i).rwY = DSGEModel.Y(DSGEModel.BVARY,LastPeriod+i-1)*ones(1,PosFH);
      ForecastResults(i).rwError = ForecastResults(i).Y-ForecastResults(i).rwY;
      ForecastResults(i).Error = ForecastResults(i).Y-ForecastResults(i).Yhat;
      ForecastResults(i).RMSE = RMSE;
      ForecastResults(i).DateVector = DateVector(LastPeriod+i:LastPeriod+i-1+PosFH);
      if max(DSGEModel.annual)>1;
         ForecastResults(i).AnnYhat = Annualize([DSGEModel.Y(DSGEModel.BVARY,LastPeriod+i+1-max(DSGEModel.annual):LastPeriod+i-1) ForecastResults(i).Yhat],DSGEModel.annual(DSGEModel.BVARY),DSGEModel.annualscale(DSGEModel.BVARY));
         ForecastResults(i).AnnError = Annualize([zeros(length(DSGEModel.BVARY),max(DSGEModel.annual)-1) ForecastResults(i).Error],DSGEModel.annual(DSGEModel.BVARY),DSGEModel.annualscale(DSGEModel.BVARY));
         ForecastResults(i).AnnrwY = Annualize(DSGEModel.Y(DSGEModel.BVARY,LastPeriod+i-max(DSGEModel.annual):LastPeriod+i-1),DSGEModel.annual(DSGEModel.BVARY),DSGEModel.annualscale(DSGEModel.BVARY))*ones(1,PosFH);
         ForecastResults(i).AnnY = Annualize(DSGEModel.Y(DSGEModel.BVARY,LastPeriod+i+1-max(DSGEModel.annual):LastPeriod+i-1+PosFH),DSGEModel.annual(DSGEModel.BVARY),DSGEModel.annualscale(DSGEModel.BVARY));
         ForecastResults(i).AnnrwError = Annualize(DSGEModel.Y(DSGEModel.BVARY,LastPeriod+i+1-max(DSGEModel.annual):LastPeriod+i-1+PosFH),DSGEModel.annual(DSGEModel.BVARY),DSGEModel.annualscale(DSGEModel.BVARY))-ForecastResults(i).AnnrwY;
      else;
         ForecastResults(i).AnnYhat = ForecastResults(i).Yhat;
         ForecastResults(i).AnnError = ForecastResults(i).Error;
         ForecastResults(i).AnnrwY = ForecastResults(i).rwY;
         ForecastResults(i).AnnrwError = ForecastResults(i).rwError;
         ForecastResults(i).AnnY = ForecastResults(i).Y;
      end;
   else;
      break;
   end;
end;

%
% end of BVARMultistepForecast.m
%
