function [m_or1,m_or2,m_aic,m_aicc,m_bic] = oe_with_oracle_and_with_AIC(Data,TestData1,TestData2,MaxOrder)
%  [m_oracle1,m_oracle2,m_aic,m_aicc,m_bic] = oe_with_oracle_and_with_AIC(Data,TestData1,TestData2,MaxOrder)
%
%  Estimate an OE model with order selection using an oracle and
%  AIC-variants
%  Data: The estimation Data
%  TestData(1,2): The "Oracle" Data set, It is used to evaluate estimated
%        OE-models of different orders. The best one is returned as
%        m_oracle(1,2)
%  MaxOrder: The maximum OE-order tested
%  m_oracle(1,2): The oracle models for oracles (TestData) 1 and 2
%  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
%  The model calculations are done under the information that we start from
%         zero initial conditions


% TestData is used only for the oracle.

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

%first estimate oe models of many orders:

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

%Then find the Order that gives the smallest AIC:
NumPar=2*Orders;
[~,aicord]= min(log(V)+2*NumPar/N);
    [~,bicord]= min(log(V)+log(N)*NumPar/N);
    [~,aiccord] =  min(log(V)+2*NumPar./(N-NumPar-1));
m_aic = m{aicord};
m_aicc = m{aiccord};
m_bic = m{bicord};


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

for k=1:MaxOrder
     
        [~,fitt1]=compare(TestData1,m{k},1,'ini','z');
        oraclefit1(k)=max(fitt1,-100);
         [~,fitt2]=compare(TestData2,m{k},1,'ini','z');
        oraclefit2(k)=max(fitt2,-100);
end

% And the oracle choice of order.

[~,ord_orac1] = max(oraclefit1);
m_or1 = m{ord_orac1};
[~,ord_orac2] = max(oraclefit2);
m_or2 = m{ord_orac2};
 