function [m_CV,ord_valid] = armax_with_CV(Data,MaxOrder)
% [m_CV,ord_valid] = armax_with_CV(Data,MaxOrder)
%    Estimate armax models of several orders and choose the order by
%    Cross-Validation
%    Data: The input-output data
%    MaxOrder: All orders up to MaxOrder are tested
%    m_CV: the chosen model
%    ord_valid: the chosen model order

[N,ny,nu]=size(Data);
Nhalf = ceil(N/2);
% Estimation Data and Validation Data:

Est = Data(1:Nhalf);
Val = Data(Nhalf:end);

%first estimate armax models of many orders

for Order=1:MaxOrder
    try
    m{Order} = armax(Est,[Order*ones(1,nu+2),ones(1,nu)]);
    catch % Order was too big
        m{Order} = m{1};
    end
end

% Validate their average prediction performance on Validation Data:
for k=1:MaxOrder
    for ph = 1:20
        [~,fitt]=compare(Val,m{k},ph,'ini','z');
        fit(ph,k)=max(fitt,-100);
    end
    validfit(k) = mean(fit(:,k));
end

[~,ord_valid] = max(validfit);

m_CV = armax(Data,[ord_valid*ones(1,nu+2),ones(1,nu)]);