function [par,parS,v,estINFO,errc] = ...
                         NAGest(mr,x,par,alg,xdmi,xdmiEST,verbose,Nit,errc)
%  FINALNI MNC-ODHAD PARAMETRU SARIMA MODELU REALNE 1-ROZMERNE CASOVE RADY
%  (Gateway to NAG Foundation toolbox routines G13AFF, G13BEF AND G13AJF)
%
% VOLANI: [par,parS,v,estINFO,errc] = ...
%                        NAGest(mr,x,par,alg,xdmi,xdmiEST,verbose,Nit,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.
% x    ... vektor pozorovani (nediferencovane) casove rady X(t) delky n,
%          xd znaci prislusnou diferencovanou radu delky nd:=n-d-D*s.
% par  ... bunecne pole (1 x 4) vsech ctyr typu pocatecnich odhadu
%          parametru (viz vystup procedury 'NAGpest'):
% alg  ... algoritmus pro vypocet odhadu parametru (nepovinny):
%          'LSQ' = metoda nejmensich ctvercu z 'G13BEF' (kef=1)
%          'MLE' = exaktni metoda maximalni verohodnosti z 'G13BEF' (kef=2)
%          'MLM' = marginalni metoda maximalni verohodnosti z 'G13BEF' (kef=3)
%          Jiny retezec (na velikosti pismen nezalezi) nebo []
%                = metoda nejmensich ctvercu z 'G13AFF'
%                = implicitni hodnota pri vynechani.
% xdmi ... (pocatecni) odhad stredni hodnoty xd: xdmi = mean(xd)
%          (nepovinny nebo []: implicitne 0)
% xdmiEST ... rizeni odhadu stredni hodnoty xdmi (nepovinny):
%             = 1: 'xdmi' pridat k 'par' jako dodatecny odhadovany parametr
%                  (implicitni pri vynechani,[] nebo jine nenulove hodnote)
%             = 0: 'xdmi' ponechat pevny
%             Celkovy pocet odhadovanych parametru je tedy
%             Npar := p + q + P + Q + xdmiEST.
% verbose ... monitorovani postupu optimalizacniho procesu na displeji
%             (nepovinny):
%             ~= 0 ... ano
%              = 0 ... ne  (implicitni pri vynechani nebo [])
% Nit  ... maximalni pocet iteraci optimalizacniho procesu, Nit>=0
%          (nepovinny nebo []: implicitne 50)
% 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 finalne odhadnutych
%            parametru (vektory po rade 1 x p, 1 x q, 1 x P, 1 x Q).
% parS   ... bunecne pole (1 x 4) smerodatnych odchylek odpovidajicich
%            parametrum z 'par': 'parS{k}' prislusi k 'par{k}', k=1,2,3,4.
% v      ... rezidualni rozptyl (odhad rozptylu generujiciho bileho sumu Z(t))
% estINFO... dalsich 9 informacnich vystupnich parametru,
%            struktura s polozkami:
%   estINFO.xdmi   ... finalni odhad stredni hodnoty xd.
%   estINFO.xdmiS  ... smerodatna odchylka odhadu stredni hodnoty xd
%                      pri xdmiEST=1, jinak [].
%   estINFO.corrmx ... korelacni matice vektoru vsech finalnich odhadu
%                      parametru (Npar x Npar)
%   estINFO.isf ... indikator uspesnosti (success/failure) finalniho odhadu
%            jednotlivych typu parametru (vektor 1 x 4), pro k=1,2,3,4:
%            isf(k) = -2 ... pocatecni odhad parametru nevyhovuje podminkam
%                            kauzality a stacionarity
%            isf(k) = -1 ... iteracni proces nekonverguje ke spravnemu reseni,
%                            nebot byl obdrzen neplatny odhad parametru
%                            (nevyhovuje podminkam kauzality a invertibility)
%            isf(k) = 0  ... zadny parametr k-ho typu v modelu
%            isf(k) = 1  ... ziskany uspokojive finalni odhady parametru
%                            k-ho typu
%   estINFO.stset ... stavova mnozina informaci pouzivana pro dalsi zpracovani
%                     v procedurach 'NAGupd' a 'NAGpred' (vektor 1 x Nst), kde
%                     Nst = P*s + D*s + d + q + max(p,Q*s).
%   estINFO.z   ... rezidua odvozena z diferencovane rady xd
%                   (odhad bileho sumu Z(t))
%   estINFO.rss ... rezidualni soucet ctvercu pouzity pri vypoctu 'v':
%                   v = rss/ndf.
%   estINFO.ndf ... pocet stupnu volnosti: ndf = n-d-D*s-p-q-P-Q-xdmiEST.
%   estINFO.itc ... pocet skutecne provedenych iteraci.
% errc  ... navratovy kod uspesnosti z procedury 'G13AFF':
%           = 0 ... uspesne provedeni
%           ~=0 ... neuspesne provedeni, kod chyby udava pricinu
%           1   ... chybny 'mr'
%           2   ... xd je prilis kratka: n-d-D*s <= length(par)+xdmiEST
%           3   ... Nit < 0
%           4   ... x je prilis kratka: n < Nst
%           5,6 ... blize nespecifikovana chyba rozhrani MATLAB->NAG
%           7   ... iteracni proces skoncil neuspesne v dusledku spatne
%                   podminenosti funkce pocitajici soucet ctvercu, nebo
%                   implicitni konvergencni kriterium je prilis striktni.
%           8   ... nelze invertovat matici nutnou k vypoctu 'corrmx'.
%           9   ... rovnice pro zpetnou predpoved nelze vyresit
%           10  ... uspokojive odhady nebyly obdrzeny pro vsechny typy
%                   parametru.
%           11  ... blize nespecifikovana chyba rozhrani MATLAB->NAG
%           #.1 ... navratovy kod # z procedury 'G13BEF'
%           #.2 ... navratovy kod # z procedury 'G13AJF'
%           Pri #.1, #.2 se vypisuje ten z informacnich textu pro chyby 1-10,
%           ktery je nejblize mozne pricine (podrobneji - viz manual k NAGu)
%           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(3,9,nargin));
if nargin<9, errc = []; end
if nargin<8, Nit = []; end
if nargin<7, verbose = []; end
if nargin<6, xdmiEST = []; end
if nargin<5, xdmi    = []; end
if nargin<4, alg     = []; end
%
if isempty(xdmi),    xdmi=0;    end
if isempty(xdmiEST), xdmiEST=1; end
if isempty(verbose), verbose=0; end
if isempty(Nit),     Nit=50;    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);
errc = ~strcmp(errc,'-');
alg = upper(alg);
%
par = [par{1},-par{2},par{3},-par{4}];
par = reshape(par,p+q+P+Q,1);
c   = xdmi;
kfc = (xdmiEST~=0)+0;
npar= p+q+P+Q+kfc;
kpiv= verbose;
nit = Nit;
kef = strcmp(alg,'LSQ') + 2*strcmp(alg,'MLE') + 3*strcmp(alg,'MLM');
if kef==0
   [par,s,ndf,sd,cm,st,nst,itc,isf,res,nres,c,ifail] = ...
                                    g13aff(mr,par,x,c,kfc,kpiv,nit,1);
v   = s/ndf;
else
   mt = [0;0;0;1];
   [par,ans,itc,sd,cm,s,d,ndf,res,sttf,nsttf,zsp,ifail1] = ...
             g13bef(mr,mt,[par;c],x,0,zeros(4,1),kfc,kpiv,kef,0,nit,1);
   c = par(p+q+P+Q+1); par(p+q+P+Q+1)=[];
   ifail = ifail1+0.1*(ifail1~=0);
   cm = cm(1:npar,1:npar);
%   if ifail1==0
   [v,st,nst,ans,ans,isf,ifail2] = g13ajf(mr,par,c,x,0,kfc,1);
   ifail = ifail2+0.2*(ifail2~=0);
% 'v' a 'st' pri 'kfc=0' a 'LSQ' maji jine hodnoty, nez z 'g13aff',
% pricemz ostatni veliciny souhlasi ? Asi chyba v 'g13ajf' ?
% V dusledku toho se pak lisi i predikce!
%   end
end
clear mex;
nres = length(res);
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)};
sd  = sd.';
parS= { sd(1:p),  sd(p+1:p+q), sd(p+q+1:p+q+P),  sd(p+q+P+1:p+q+P+Q)};
if kfc, xdmiS=sd(length(sd)); else xdmiS=[]; end
estINFO = struct('xdmi',c,'xdmiS',xdmiS,'corrmx',cm,'isf',isf.',...
                'stset',st(1:nst).','z',res(1:nres),'rss',s,...
                'ndf',ndf,'itc',itc);
%
if errc
   errn=['errc=',num2str(ifail),': '];
   switch ifail
   case {1,1.1,2.1,1.2}, error([errn,'invalid model range ''mr''']);
   case {3.1,10.1} , error([errn,'parameters do not satisfy causality or ',...
                                 'invertibility conditions']);
   case 4.1, error([errn,'some criteria controlling the iterative ',...
                         'process are invalid']);
   case {2,2.2}, error([errn,'the differenced time series is too short']);
   case 3, error([errn,'negative number of iterations supplied']);
   case 4, error([errn,'the time series is too short']);
   case {7,12.1}, error([errn,'the serach procedure in the algorithm ',...
                              'has failed']);
   case {8,9.1}, disp([char(7),char(7),'??? Warning using ==> nagest']);
                 disp([errn,'the covariance matrix cannot be computed']);
   case {9,8.1,5.2}, disp([char(7),char(7),'??? Warning using ==> nagest']);
              disp([errn,'the equations giving the latest backforecasts ',...
                         'cannot be solved']);
   case {10,6.2}, disp([char(7),char(7),'??? Warning using ==> nagest']);
              disp([errn,'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;
