function [par,v,isf,errc] = NAGpest(mr,xdro,xdv,errc)
%    PREDBEZNY ODHAD PARAMETRU SARIMA MODELU REALNE 1-ROZMERNE CASOVE RADY
%         (Gateway to NAG Foundation toolbox routine G13ADF)
% VOLANI: 1) [par,v,isf,errc] = NAGpest(mr,xdro,xdv,errc)
%         2) [par,v,isf,errc] = NAGpest(mr,x,[],errc) nebo
%            [par,v,isf,errc] = NAGpest(mr,x,errc)
% VSTUP:
% mr   ... =[p,d,q,P,D,Q,s] ... vektor radu modelu SARIMA(p,d,q,P,D,Q,s)
%          vytvareny pomocnymi procedurami AR,MA,ARMA,ARIMA,SARIMA:
%          p ... rad nesezonni AR casti modelu (p>=0, cele cislo)
%          d ... rad nesezonniho diferencovani (d>=0, cele cislo)
%          q ... rad nesezonni MA casti modelu (q>=0, cele cislo)
%          P ... rad sezonni AR casti modelu (P>=0, cele cislo)
%          D ... rad sezonniho diferencovani (D>=0, cele cislo)
%          Q ... rad sezonni MA casti modelu (Q>=0, cele cislo)
%          s ... perioda (pocet bodu) sezonni komponenty (s>=0, cele cislo)
%          Podminky korektnosti (jsou kontrolovany - viz errc=1):
%          p,d,q,P,D,Q,s>=0; p+q+P+Q>0; s~=1; s==0 => P+D+Q=0; s>1 => P+D+Q>0.
% xdro ... odhad autokorelacni funkce (vektor nq x 1):
%          xdro(h) = cov(xd(t+h),xd(t))/var(xd(t)) pro h=1,2,...,nq, kde
%          xd = NAGdiff(mr,x) je rada 'x' po odpovidajicim diferencovani a
%          nq >= nqmin:=max(p+q,s*(P+Q)).
% xdv  ... odhad rozptylu xd: xdv = std(xd)^2
% x    ... vektor pozorovani casove rady X(t) delky n.
%          Zadava se ve variante 2), kdy procedura 'NAGpest' spocte
%          'xdro' a 'xdv' sama pomoci 'NAGdiff' a 'NAGacf' pro
%          optimalni delku nqopt, nqopt=max(nqmin,length(xd)/2):
%          [xdro,xdv] = NAGacf(NAGdiff(mr,x),nqopt,errc).
%          Pokud length(xd)=n-(d+D*s)<=nqmin, pak 'NAGacf' hlasi chybu
%          errc=1 a 'NAGpest' chybu errc=2.1. Toto ve variante 2) signalizuje
%          prilis mnoho parametru v modelu vzhledem k delce casove rady 'n'.
%          Problem je mozno odstranit pouze zvysenim delky casove rady nebo
%          snizenim hodnot nekterych radu modelu v 'mr'.
% errc ... rizeni vypisu chyboveho hlaseni (nepovinny znakovy parametr):
%          '+' ... povoleni vypisu vysvetlujiciho chyboveho hlaseni
%                  pri vyskytu chyby (implicitni hodnota pri vynechani)
%          '-' ... potlaceni vypisu chyboveho hlaseni
%          Pouziti jakehokoli jineho retezce ma stejny efekt jako '+'.
%
% VYSTUP:
% par   ... bunecne pole (1 x 4) vsech ctyr typu predbezne odhadnutych
%           parametru (vektory po rade 1 x p, 1 x q, 1 x P, 1 x Q):
%           par{1}=phi   ... nesezonni AR parametry (vektor 1 x p)
%           par{2}=theta ... nesezonni MA parametry (vektor 1 x q)
%           par{3}=PHI   ... sezonni   AR parametry (vektor 1 x P)
%           par{4}=THETA ... sezonni   MA parametry (vektor 1 x Q)
% v     ... rezidualni rozptyl (odhad rozptylu generujiciho bileho sumu Z(t))
% isf   ... indikator uspesnosti (success/failure) predbeznych odhadu
%           jednotlivych typu parametru (vektor 1 x 4), pro k=1,2,3,4:
%           isf(k) = -1 ... uspokojive predbezne odhady parametru k-ho typu
%                           nebylo mozno ziskat (nutna korekce modelu).
%           isf(k) = 0  ... zadny parametr k-ho typu v modelu
%           isf(k) = 1  ... ziskany uspokojive predbezne odhady parametru
%                           k-ho typu
% errc  ... navratovy kod uspesnosti:
%           = 0 ... uspesne provedeni
%           ~=0 ... neuspesne provedeni, kod chyby udava pricinu
%           1   ... chybny 'mr': d<0 nebo D<0 nebo s<0 nebo (s=0 pro D>0)
%           1.x ... chyba x v 'NAGdiff' ve variante 2)
%           2   ... xdro je prilis kratka: nq < max(p+q,s*(P+Q))
%           2.x ... chyba x z 'NAGacf' ve variante 2)
%           3   ... |xdro(h)|>1 pro nektere h=1,2,...,nq
%           4   ... xdv <= 0
%           5,6 ... blize nespecifikovana chyba rozhrani MATLAB->NAG
%           7   ... uspokojive odhady nebyly obdrzeny pro vsechny typy
%                   parametru (isv(k)=-1 pro nektere k=1,2,3,4).
%           Je-li procedura volana z programu, mel by byt vystupni kod vzdy
%           testovan, nebot nasledny vypocet by jinak byl nekorektni.

% (C) 1999, V. Vesely, Masaryk University of Brno (Czech Republic)
%
error(nargchk(2,4,nargin));
if nargin<4,  errc = []; end
if nargin<3,  xdv = []; end
if nargin==3 & ischar(xdv), errc=xdv; xdv=[]; end
%
if any(mr<0), error('errc=1: invalid model range ''mr'''); end
p=mr(1); q=mr(3); P=mr(4); Q=mr(6); s=mr(7);
nqmin = max(p+q,s*(P+Q));
%
% xdv empty => variant 2)
ifail = 0;
if isempty(xdv)
   [xd,xr,ifail] = NAGdiff(mr,xdro,'-');
   ifail = (ifail~=0)+0.1*ifail;
   if ~ifail
      nqopt = max(nqmin,ceil(length(xd)/2)-1);
      [xdro,xdv,mi,stat,ifail] = NAGacf(xd,nqopt,'-');
      ifail = 2*(ifail~=0)+0.1*ifail;
   end
end
%
errc = ~strcmp(errc,'-');
if ifail
   par=[]; v=[]; isf=[];
else
   [par,v,isf,ifail] = g13adf(mr,xdro,xdv,1);
   par = par.';
   par = {par(1:p),-par(p+1:p+q),par(p+q+1:p+q+P),-par(p+q+P+1:p+q+P+Q)};
   isf = isf.';
end
clear mex;
if errc
   switch ifail
   case 1, error('errc=1: invalid model range ''mr''');
   case 2, error(['errc=2: autocorrelation function ''xdro'' ',...
                  'of differenced time series is too short']);
   case 2.1, error('errc=2.1: too many parameters in the model');
   case 3, error('errc=3: invalid ''xdro'' due to |xdro(h)|>1 for some h');
   case 4, error('errc=4: nonpositive varianace ''xdv <= 0''');
   case 7, disp([char(7),char(7),'??? Warning using ==> nagpest']);
           disp(['errc=7: satisfactory parameter estimates could not be ',...
                  'obtained for all parameter types']);
   otherwise
     if ifail~=0
        error(['errc=',num2str(ifail),': other error']);
     end
   end
end
errc = ifail;
