function fcast = VARfcast(y,p,h,B)
% fcast = VARfcast(y,p,h)
% VARfcast uses equation-by-equation OLS to generate 1- through h-step-ahead
%   forecasts for a VAR(p) model of the matrix of time series y.
% The n columns of y are separate time series; the T rows of y are the
%   observations for these series.
% fcast is an h x n matrix of forecasts for each of the n variables at
%   horizons of 1 to h steps ahead.
% B is a matrix of VAR coefficients.  If B is empty or missing, then a
%   VAR of the appropriate lag lenth is estimated on the data y.
%
% Eric Swanson, 2005.

if (nargin<4); B = []; end ;
if (nargin<3); h = 1; end ;
if (nargin<2); p = 1; end ;
[T,n] = size(y) ;

for i=0:p-1 ;
  X(:,i*n+1:(i+1)*n) = y(p-i:T-i,:) ; % X is matrix of lags of data y
end
X(:,p*n+1) = ones(T-p+1,1) ;

if (isempty(B)); B = ols(y(p+1:T,:),X(1:T-p,:),0,0); end ; % eqn-by-eqn OLS, if needed

for i=1:h ;
  X(end+1,:) = [X(end,:)*B, X(end,1:(p-1)*n), 1] ; % iterate the VAR forecast
end ;

fcast = X(end-h+1:end,1:n) ;


