function [m_oracle,m_aic,m_aicc,m_bic] = armax_with_oracle_and_with_AIC(Data,TestData,MaxOrder)
%  [m_oracle,m_aic,m_aicc,m_bic] = armax_with_oracle_and_with_AIC(Data,TestData,MaxOrder)
%
%  Estimate an ARMAX model with order selection using an oracle and
%        AIC-variants
%  Data: The estimation Data
%  TestData: The "Oracle" Data set, It is used to evaluate estimated
%        ARMAX-models of different orders. The best one is returned as
%        m_oracle
%  MaxOrder: The maximum ARMAX-order tested
%  m_oracle: The oracle model
%  m_aic:The model that minimizes AIC on estimation data
%  m_aicc:The model that minimizes AICc on estimation data
%  m_bic:The model that minimizes BIC on estimation data

% TestData is used only for the oracle.

[N,ny,nu]=size(Data);

%first estimate armax models of many orders:


Orders = 1:MaxOrder;
for Order=Orders
    m{Order} = armax(Data,[Order*ones(1,nu+2),ones(1,nu)]);
    V(Order) = m{Order}.Report.Fit.LossFcn;
end

npar = (nu+3)*Orders; % number of parameters = na+nu*nb+nc+na(initial conditions) 
%Then find the Order that gives the smallest AIC:
[~,aicord]= min(log(V)+2*npar/N); %AIC
[~,bicord]= min(log(V)+log(N)*npar/N); %BIC
[~,aiccord] =  min(log(V)+2*npar./(N-npar-1)); %AICc

m_aic = m{aicord};
m_aicc = m{aiccord};
m_bic = m{bicord};


%Now compute the prediction fits to the oracle TestData, all predictione
%horizons up to 20, and the average fit for all horizons

for k=1:MaxOrder
    for ph = 1:20
        [~,fitt]=compare(TestData,m{k},ph,'ini','z');
        fit(ph,k)=max(fitt,-100);
    end
    oraclefit(k) = mean(fit(:,k));
end

% And the oracle choice of order.

[~,ord_orac] = max(oraclefit);
m_oracle = m{ord_orac};
