function [p,yy_ff]=arlagselect_forecast(x,P,CC,hmax)
%
% Luca Benati/FabrizioVenditti
%
% function p=arlagselect(x,P,XX,CC,F)
% This program selects the optimal lag order for an AR(p) process based on AIC or SIC.
%                             Input of the program is:
% x     = a Tx1 vector
% P     = the maximum order of the autoregression the algorithm has to consider
% CC    = the criterion for selecting the lag order: must be set equal to 'AIC' or 'SIC'
% hmax  = steps ahead to forecast
%                             Output of the program is:
% p     = the optimal lag order
% yy_ff = forecast hmax steps ahead;
%

%
N=length(x);
AIC=[];
SIC=[];
for j=1:P
    Y=x(j+1:N);
    NY=length(Y);
    X(:,1)=ones(NY,1);
    for h=1:j
	    X(:,h+1)=x(j+1-h:N-h);
    end
    bhat=(X'*X)\(X'*Y);
    uhat=Y-X*bhat;
	s2=(uhat'*uhat)/(NY-length(bhat));
    S=s2*ones(NY,1);
	L=-0.5*log(S)-0.5*((uhat.^2)./S);
	aic=NY*log(s2)+2*j;
	sic=log(s2)+log(NY-P)/(NY-P);
    AIC=[AIC' aic]';
    SIC=[SIC' sic]';
    clear Y
    clear X
    clear L
end
%
if CC=='AIC'
    p=(1:1:P)*(AIC==min(AIC));
else
    p=(1:1:P)*(SIC==min(SIC));
end

% now re-estimate the model with optimal p and forecast hmax steps ahead
    j=p;
    Y=x(j+1:N);
    NY=length(Y);
    X(:,1)=ones(NY,1);
    for h=1:j
	    X(:,h+1)=x(j+1-h:N-h);
    end

    bhat=(X'*X)\(X'*Y);
    Xf = X(end,2:end); % cut out the constant and take the last obs
    Yf = Y(end);
    yy_ff = NaN(hmax,1);
    for jj=1:hmax
        Xf = [Yf Xf(1,1:p-1)]; % if  p=1 then Xf(1,1:p-1) is an empty matrix
        Yf = [1 Xf]*bhat;
        yy_ff(jj,1) = Yf;
    end
    
    
