function [ar,ma,ASAsellog,ASAcontrol]=armasel_rs(varargin)
%ARMASEL_RS ARMAsel model identification
%   [AR,MA,SELLOG] = ARMASEL_RS(ARH,N_OBS) estimates autoregressive, moving 
%   average, and autoregressive moving average models from the high-order 
%   AR model ARH and selects the model with optimal predictive qualities.
%   ARH has been estimated from N_OBS observations. The AR and MA parts of
%   the selected model, each possibly of order 0, are returned in the
%   parameter vectors AR and MA. The structure SELLOG provides additional
%   information on the selection process.
%   
%   SELLOG contains the fields 'ar', 'ma' and 'arma', in which SELLOG 
%   structures are nested, as returned by the functions ARH2AR, ARH2MA 
%   and ARH2ARMA, invoked by ARMASEL. In the field 'armasel' a structure 
%   is nested that reports information about the final stage of model 
%   selection, where the preselected AR, MA and ARMA models are compared.
%   
%   ARMASEL(ARH,N_OBS,CAND_AR_ORDER,CAND_MA_ORDER,CAND_ARMA_ORDER)
%   narrows the selection to candidate models with orders provided by the 
%   rows CAND_AR_ORDER, CAND_MA_ORDER and CAND_ARMA_ORDER AR, MA and
%   ARMA(r,r-1) models independently.
%   For any of these arguments it is allowed to pass an 
%   empty array. Alternatively, additional arguments may be omitted from 
%   the input list. In both cases, default values are automatically 
%   determined and substituted for the missing arguments. The functions 
%   ARH2AR, ARH2MA and ARH2ARMA provide additional information on 
%   defining candidate orders. Note in this respect, that the candidate 
%   AR orders of the ARMA model are called CAND_ARMA_ORDER in this help 
%   text, while in ARH2ARMA they are called CAND_AR_ORDER.
%
%   The selection of MA and ARMA models can be conditioned to the 
%   selection of AR models from a specific set of candidate orders 
%   CAND_AR_ORDER. See ASAGLOB_AR_COND for more information.
%   
%   ARMASEL_RS is an ARMASA_RS main function.
%   ARMASEL_RS can also be used to selected models based on an ARH model
%   that has been estimated from segments of data.
%   
%   See also: ARMASEL, ARH2AR, ARH2MA, ARH2ARMA, DATA_SEGMENTS.

%   Reference:  P. M. T. Broersen and S. de Waele, Selection of Order and
%               Type of Time Series Models Estimated from Reduced Statistics
%               Proceedings of SYSID 2002, May 2002.

%Header
%===================================================================================================

%Declaration of variables
%------------------------

%Declare and assign values to local variables
%according to the input argument pattern
[ar_rs,n_obs,cand_ar_order,cand_ma_order,cand_arma_order,arma_order_diff,ASAcontrol] = ASAarg(varargin, ...
{'ar_rs'       ;'n_obs'      ;'cand_ar_order';'cand_ma_order';'cand_arma_order';'arma_order_diff';'ASAcontrol'}, ...
{'isnumeric'   ;'isnumeric'  ;'isnumeric'    ;'isnumeric'    ;'isnumeric'      ;'isnumeric'      ;'isstruct'  }, ...
{'ar_rs'       ;'n_obs'                                                                                       }, ...
{'ar_rs'       ;'n_obs'      ;'cand_ar_order'                                                                 }, ...
{'ar_rs'       ;'n_obs'      ;'cand_ar_order';'cand_ma_order'                                                 }, ...
{'ar_rs'       ;'n_obs'      ;'cand_ar_order';'cand_ma_order';'cand_arma_order'                               }, ...
{'ar_rs'       ;'n_obs'      ;'cand_ar_order';'cand_ma_order';'cand_arma_order';'arma_order_diff'             });

%Declare ASAglob variables 
ASAglob = {'ASAglob_subtr_mean';'ASAglob_mean_adj';'ASAglob_rc';'ASAglob_ar';'ASAglob_final_f'; ...
      'ASAglob_final_b';'ASAglob_ar_cond'};

%Assign values to ASAglob variables by screening the
%caller workspace
for ASAcounter = 1:length(ASAglob)
   ASAvar = ASAglob{ASAcounter};
   eval(['global ' ASAvar]);
   if evalin('caller',['exist(''' ASAvar ''',''var'')'])
      eval([ASAvar '=evalin(''caller'',ASAvar);']);
   else
      eval([ASAvar '=[];']);
   end
end

%ARMASA-function version information
%-----------------------------------

%This ARMASA-function is characterized by
%its current version,
ASAcontrol.is_version = [2000 12 30 20 0 0];
%and its compatability with versions down to,
ASAcontrol.comp_version = [2000 12 30 20 0 0];

%This function calls other functions of the ARMASA
%toolbox. The versions of these other functions
%must be greater than or equal to:
ASAcontrol.req_version.sig2ar = [2000 12 30 20 0 0];
ASAcontrol.req_version.sig2ma = [2000 12 30 20 0 0];
ASAcontrol.req_version.sig2arma = [2000 12 30 20 0 0];

%Checks
%------

if ~isfield(ASAcontrol,'error_chk') | ASAcontrol.error_chk
      %Perform standard error checks
   %Input argument format checks
   ASAcontrol.error_chk = 1;
   if ~isnum(ar_rs)
      error(ASAerr(11,'ar_rs'))
   elseif ~isvector(ar_rs)
      error([ASAerr(14) ASAerr(15,'ar_rs')])
   end
   if ~isempty(cand_ar_order)
      if ~isnum(cand_ar_order) | ~isintvector(cand_ar_order) |...
            cand_ar_order(1)<0 | ~isascending(cand_ar_order)
         error(ASAerr(12,{'candidate';'cand_ar_order'}))
      elseif size(cand_ar_order,1)>1
         cand_ar_order = cand_ar_order';
         warning(ASAwarn(25,{'column';'cand_ar_order';'row'},ASAcontrol))
      end
   end
   if ~isempty(cand_ma_order)
      if ~isnum(cand_ma_order) | ~isintvector(cand_ma_order) |...
            cand_ma_order(1)<0 | ~isascending(cand_ma_order)
         error(ASAerr(12,{'candidate';'cand_ma_order'}))
      elseif size(cand_ma_order,1)>1
         cand_ma_order = cand_ma_order';
         warning(ASAwarn(25,{'column';'cand_ma_order';'row'},ASAcontrol))
      end
   end
   if ~isempty(cand_arma_order)
      if ~isnum(cand_arma_order) | ~isintvector(cand_arma_order) |...
            cand_arma_order(1)<0 | ~isascending(cand_arma_order)
         error(ASAerr(12,{'candidate';'cand_arma_order'}))
      elseif size(cand_arma_order,1)>1
         cand_arma_order = cand_arma_order';
         warning(ASAwarn(25,{'column';'cand_arma_order';'row'},ASAcontrol))
      end
   end
   if ~isempty(arma_order_diff) & ...
         (~isnum(arma_order_diff) | ...
         ~isintscalar(arma_order_diff) |...
         arma_order_diff<0)
      error(ASAerr(17,'arma_order_diff'))
   end

   %Input argument value checks
   if ~isreal(ar_rs)
      error(ASAerr(13))
   end
   if max(cand_ar_order) > length(ar_rs)-1
      error(ASAerr(37,'cand_ar_order'))
   end
   if max(cand_ma_order) > length(ar_rs)-1
      error(ASAerr(37,'cand_ma_order'))
   end
   if ~isempty(cand_arma_order) & ...
         ~isempty(arma_order_diff)
      if cand_arma_order(1)~=0 & ...
            (arma_order_diff < 1 | ...
            arma_order_diff > cand_arma_order(1))
         error(ASAerr(18,{'arma_order_diff';'1';...
               num2str(cand_arma_order(1))}))
      elseif length(cand_arma_order)>1 & ...
            (arma_order_diff < 1 | ...
            arma_order_diff > cand_arma_order(2))
         error(ASAerr(18,{'arma_order_diff';'1';...
               num2str(cand_arma_order(2))}))
      end
   end
end

if ~isfield(ASAcontrol,'version_chk') | ...
      ASAcontrol.version_chk %Perform version check
   ASAcontrol.version_chk = 1;
      
   %Make sure the requested version of this function
   %complies with its actual version
   ASAversionchk(ASAcontrol);
   
   %Make sure the requested versions of the called
   %functions comply with their actual versions
   arh2ar(ASAcontrol);
   arh2ma(ASAcontrol);
   arh2arma(ASAcontrol);
end

if ~isfield(ASAcontrol,'run') | ASAcontrol.run
   ASAcontrol.run = 1;
   ASAdate = now;
end

if ASAcontrol.run %Run the computational kernel
   ASAcontrol.version_chk = 0;
   ASAcontrol.error_chk = 0;

%Main   
%=====================================================
   
%Initialization of variables
%---------------------------

if ASAglob_ar_cond
   ASAglob_ar_cond = 1;
else
   ASAglob_ar_cond = 0;
end

n_obs_tot = sum(n_obs);

%AR-, MA- and ARMA-model identification
%--------------------------------------

[ar_ar,ar_sellog] = arh2ar...
   (ar_rs,n_obs,cand_ar_order,ASAcontrol);
[ma_ma,ma_sellog] = arh2ma...
   (ar_rs,n_obs,cand_ma_order,ASAcontrol);
[arma_ar,arma_ma,arma_sellog] = arh2arma...
   (ar_rs,n_obs,cand_arma_order,arma_order_diff,ASAcontrol);

%Selection of the ARMAsel model
%------------------------------

%Asess the selected model orders
order = ...
   [length(ar_ar)-1;length(ma_ma)-1;length(arma_ar)-1];

%Asess the corresponding prediction error estimates
sel_location = ...
   [find(order(1) == ar_sellog.cand_order);...
    find(order(2) == ma_sellog.cand_order);...    
    find(order(3) == arma_sellog.cand_ar_order)];
pe_est = [ar_sellog.pe_est(sel_location(1));...
      ma_sellog.pe_est(sel_location(2));...
      arma_sellog.pe_est(sel_location(3))];

%Select the model with the smallest prediction error
%estimate
[sel_pe_est,model] = min(pe_est);

%Arranging output arguments
%--------------------------

%Retrieve the parameters of the selected model
ar = 1;
ma = 1;
switch model
case 1 %The AR model has been selected
   ar = ar_ar;
case 2 %The MA model has been selected
   ma = ma_ma;
case 3 %The ARMA model has been selected
   ar = arma_ar;
   ma = arma_ma;
end

%Gernerate a structure variable ASAsellog to report
%the selection process
ASAsellog.funct_name = mfilename;
ASAsellog.funct_version = ASAcontrol.is_version;
ASAsellog.date_time = ...
   [datestr(ASAdate,8) 32 datestr(ASAdate,0)];
ASAsellog.comp_time = ar_sellog.comp_time+...
   ma_sellog.comp_time+...
   arma_sellog.comp_time;
ASAsellog.armasel.ar = ar;
ASAsellog.armasel.ma = ma;
ASAsellog.armasel.ar_pe_est = pe_est(1);
ASAsellog.armasel.ma_pe_est = pe_est(2);
ASAsellog.armasel.arma_pe_est = pe_est(3);
ASAsellog.armasel.ar_cond = ASAglob_ar_cond;
ASAsellog.armasel.mean_adj = ASAglob_mean_adj;
ASAsellog.ar = ar_sellog;
ASAsellog.ma = ma_sellog;
ASAsellog.arma = arma_sellog;

%Footer
%=====================================================

else %Skip the computational kernel
   %Return ASAcontrol as the first output argument
   if nargout>1
      warning(ASAwarn(9,mfilename,ASAcontrol))
   end
   ar = ASAcontrol;
   ASAcontrol = [];
end

%Program history
%======================================================================
%
% Version                Programmer(s)          E-mail address
% -------                -------------          --------------
% former versions        P.M.T. Broersen        broersen@tn.tudelft.nl
%                        S. de Waele            waele@tn.tudelft.nl
% [2000 11  1 12 0 0]    W. Wunderink           wwunderink01@freeler.nl
% [2000 12 30 20 0 0]         ,,                          ,,
