function [fe] = forecast_errors(y,yhat)
% calculates matrix of forecast errors

% inputs
% y: (Tx1) vector of data y(t), t= 1,2,...,T
% yhat: (TxH) matrix of forecasts y(t+h|t), t= 1,2,...,T, h=1,2,...,H
%        i.e. element yhat(t,h) contains forecast  y(t+h|t)

% output
% fe: (TxH) matrix of forecast errors FOR a certain period t
%      i.e. element fe(t,h) contains forecast error y(t) - y(t|t-h)
%      elements in upper right corner are NaNs
%      for instance, fe(1,1) is NaN, because forecast y(1|0) is unknown
%                    fe(1,2) and fe(2,2) are NaN, because forecasts y(1|-1) and y(2|0) are unknown
%                    etc.

T = numel(y);
Tcheck = size(yhat,1);
if T~= Tcheck
    error('Length of observation vector does not correspond to length of forecast matrix! ');
end
H = size(yhat,2);

fe = NaN(T,H);
for h = 1:H
    fe(h+1:T,h) = y(h+1:T) - yhat(1:T-h,h);
end
    