function [y_ff,rss]=var_ols(lags,hmax,m,TT,y_l,x_l)

% /* INPUT:  y is the data (TTxn matrix)
big_X = x_l'*x_l;
big_Y = x_l'*y_l;
b     = big_X\big_Y;

%      RSS
rss  = diag((y_l-x_l*b)'*(y_l-x_l*b)/(TT));
%      now perform recursive forecasts
y_f = y_l(TT,:);
x_f = x_l(TT,:);
y_ff = zeros(hmax,m);
for llags = 1:hmax
    x_f  = [y_f x_f(1,1:m*lags-m) 1];
    y_f  = x_f*b;
    y_ff(llags,:) = y_f;
end;    
% %%%%%%%%%%%%%%%%