function [pqo,v,vv,ke,pq,rp] = NAGpqest(mrmax,x,xv,fig)
%         ODHAD RADU p,q SARIMA MODELU REALNE 1-ROZMERNE CASOVE RADY
%   (Using routines: 'pacf', 'NAGacf', 'NAGpest', 'NAGest' and 'armasim')
% VOLANI: [pqo,v,vv,ke,pq,rp] = NAGpqest(mrmax,x,xv,fig)
% VSTUP:
% mrmax... =[pmax,d,qmax,P,D,Q,s], kde
%          pmax ... horni mez pro p, pmax>=0 (pmax=0 => jen MA(q) modely),
%          d    ... rad nesezonniho diferencovani (d>=0, cele cislo)
%          qmax ... horni mez pro q, qmax>=0 (qmax=0 => jen AR(p) modely),
%          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)
%          0 < pmax + qmax + P + Q + 1 < n-d-D*s, kde n=length(x).
%          s~=1; s==0 => P+D+Q=0; s>1 => P+D+Q>0.
% x    ... usek pozorovani casove rady X(t) delky n, jenz bude pouzit
%          pro odhad parametru pro nasledujici kombinace radu (p,q)
%          usporadane vzestupne dle poctu parametru p+q
%          (poradovy index znacen dale jako 'k'):
%          pri qmax=0: q=0, p = 1,2,...,pmax
%          pri qmax>0: p=0,1...,pmax; q=1,2,...,qmax
% xv   ... validacni usek pozorovani navazujici na 'x' (nepovinny nebo [])
% fig  ... cislo obrazku ('handle') pro vykresleni prubehu ztratovych
%          funkci: v,FPE,AIC,MDL a pripadne vv (v samostatnem obrazku)
%          (nepovinny: implicitne se grafy nekresli)
%
% VYSTUP:
% pqo  ... p=pqo(1,j), q=pqo(2,j), j=1,2,3,resp.4 udava optimalni rady urcene
%          j-tou ztratovou funkci (kriteriem) v tomto poradi:
%          a) na 'x':
%             j=1: FPE odpovidajici minimu vfpe(k)=v(k)*(1+rp(k))/(1-rp(k))
%             j=2: AIC odpovidajici minimu vaic(k)=v(k)*(1+2*rp(k))
%             j=3: MDL odpovidajici minimu vmdl(k)=v(k)*(1+log(n)*rp(k))
%          b) na 'xv' (pokud bylo 'xv' zadano)
%             j=4: minimum vv(k) - rezidualni rozptyl
%          pro k=1,2,...,ke.
% v  ... v(k)  = rezidualni rozptyl na diferencovane rade 'x'  pro k-ty model
% vv ... vv(k) = rezidualni rozptyl na diferencovane rade 'xv' pro k-ty model
%        vv==[], pokud 'xv' nebylo zadano.
% ke ... k=1,...,ke udava pocatecni usek posloupnosti modelu, jez mohly
%        byt uspesne odhadnuty, pro (ke+1)-ty model odhad poprve selhal
%        i pri pouziti 'NAGest' misto 'NAGpest';
%        pro takova 'k' je v(k)=vv(k)=NaN neplatna hodnota.
% pq ... pq(:,k) udava dvojici radu (p,q) v k-tem modelu, tj.:
%        p=pq(1,k) a q=p(2,k).
% rp ... rp(k) = (p+q)/n pocet parametru pripadajici na jedno pozorovani
%        v k-tem modelu.
%
% POZNAMKA: 'pq' a 'rp' umoznuji spolu s 'ke','v', resp. 'vv' spocist
%           dodatecne dalsi kriteria, napr. typu v*(1+c*rp) pro vhodnou
%           vahovou konstantu 'c'.

% (C) 1999, V. Vesely, Masaryk University of Brno (Czech Republic)
%
error(nargchk(2,4,nargin));
if nargin<4, fig=[]; end
if nargin<3, xv =[]; end

mrmax = round(mrmax);
pmax = mrmax(1); d = mrmax(2); qmax  = mrmax(3);
P = mrmax(4); D = mrmax(5); Q = mrmax(6); s = mrmax(7);
n = length(x); nv = length(xv);
x = reshape(x,n,1); xv = reshape(xv,nv,1);
nd = n-d-D*s;
if any(mrmax<0)|s==1|(s==0 & P+D+Q>0)|(s>1 & P+D+Q==0)|pmax+qmax+P+Q+1>=nd
   error('Invalid range limits');
end
usexv = ~isempty(xv);           % & P+Q==0;
if usexv
   xd = NAGdiff([d,D,s],[x;xv]);
   xdv = xd(nd+1:length(xd));
   xd  = xd(1:nd);
else
   xd = NAGdiff([d,D,s],x);
end
nd = length(xd);
[xdro,xdvar,xdmi0] = NAGacf(xd,nd-1);
%
vv = [];
if qmax+P+Q==0       % AR model diferencovane rady
   ke = pmax;
   Npq = 1:pmax; pq =[Npq;zeros(size(Npq))];
   xpacf=[]; phi=[]; vp=[];
   for p=Npq
       [xpacf,phi,phiE,v(p)] = pacf(xdvar*[1;xdro(1:p)],n,[],xpacf,phi,vp);
       vp = v(p);
       if ~isempty(xv)
          [xvp,z,vv(p)] = armasim(xv,phi,[],x);
       end
   end
else              % MA model nebo smiseny ARMA model
%-----------------------------------------------------------------------------
% Vypocet dvojic [p,q] v rozsahu p=0,1,...,pmax; q=1,2,...,qmax:
  if qmax==0   % AR model + sezonni AR model (pro P+Q>0)
     Npq = 1:pmax; pq =[Npq;zeros(size(Npq))];
     phi0 = 0; theta0 = [];
     p0 = 1; q0 = 0;
  else
     pq  = [kron(0:pmax,ones(1,qmax));kron(ones(1,pmax+1),1:qmax)];
     Npq = sum(pq);  % pocty odhadovanych parametru p+q pro kazdou dvojici [p,q]
     [Npq,I] = sort(Npq); % pq setridime vzestupne podle poctu parametru
     pq = pq(:,I);
     phi0 = []; theta0 = 0;
     p0 = 0; q0 = 1;
  end
  kmax = size(pq,2);   % pocet dvojic
% Odhad v(k) z useku x a vv(k) z useku xv pro kazdou dvojici pq(:,k):
  ke = kmax;  % nejvyssi index verohodne dvojice (p,q)
  PHI0 = zeros(1,P); THETA0 = zeros(1,Q);
  mr = mrmax;
  for k=1:kmax
     p = pq(1,k); q = pq(2,k); mr(1)=p; mr(3)=q;
     [par,v(k),isf,errc] = NAGpest(mr,xdro,xdvar,'-');
     xdmi = xdmi0;
     if errc
%        disp(['with p=',num2str(p),' q=',num2str(q)]);
%        disp('''NAGest'' will be used instead of ''NAGpest''');
        phi = phi0(1:min(p,p0)); theta = theta0(1:min(q,q0));
        if p>p0, phi(p)   = 0; end
        if q>q0, theta(q) = 0; end
%     else
%        p0 = p; q0 = q;
%        phi0 = par{1}; theta0 = par{2}; phi = phi0; theta = theta0;
%        PHI0 = par{3}; THETA0 = par{4};
%     end
        [par,parS,v(k),estINFO,errc] = ...
             NAGest(mr,x,{phi,theta,PHI0,THETA0},[],xdmi0,[],[],[],'-');
        if errc
          [par,parS,v(k),estINFO,errc] = ...
               NAGest(mr,x,{phi,theta,PHI0,THETA0},'MLE',xdmi0,[],[],[],'-');
        end
        xdmi = estINFO.xdmi;
        if errc
           disp(['The loss function could not be evaluated ',...
                 'with p=',num2str(p),', q=',num2str(q)]);
           disp('--------------------------------------------------------');
           v(k) = NaN; if ~isempty(xv), vv(k) = NaN; end
           if ke==kmax, ke = k-1; end
        end
     end
     if usexv & ~isnan(v(k))
% Puvodni varianta:
%        [ans,zd] = armasim(xd-xdmi,par{1},par{2});
%        [xdvp,ans,vv(k)] = armasim(xdv-xdmi,par{1},par{2},xd-xdmi,zd);
        [phix,thetax,xdc,phid]=ARMApar(mr,par,xdmi);
% Varianta s diferencovanou radou:
        pd = length(phid);
        [ans,zd] = armasim(xd(pd+1:nd),phid,thetax,xd(1:pd));
        [xdvp,ans,vv(k)] = armasim(xdv,phid,thetax,xd,zd);
% Varianta s nediferencovanou radou:
%        px = length(phix);
%        [ans,z] = armasim(x(px+1:n),phix,thetax,x(1:px));
%        [xvp,ans,vv(k)] = armasim(xv,phix,thetax,x,z);
     end
     if ~isnan(v(k)) & p>p0, phi0   = par{1}; p0=p; end
     if ~isnan(v(k)) & q>q0, theta0 = par{2}; q0=q; end
     if ~isnan(v(k)), PHI0   = par{3}; end
     if ~isnan(v(k)), THETA0 = par{4}; end
  end
end
%
% Vykreslime prubeh v(k) a prislusnych kriterialnich funkci FPE, AIC a MDL:
rp   = Npq./n; AIC = 2; MDL = log(n);
vfpe = v.*(1+rp)./(1-rp);
vaic = v.*(1+AIC*rp);
vmdl = v.*(1+MDL*rp);
if ~isempty(fig)
   figure(fig);
   plot(v,'w-');
   hold on; plot(vfpe,'r--'); plot(vaic,'g--'); plot(vmdl,'g:');
   title('Loss function (-), FPE (r--), AIC (g--) a MDL (g:)');
   xlabel('index of (p,q) sorted in ascending order of p+q');
   hold off;
   if ~isempty(vv)
      figure;
      figp = get(gcf,'Position');
      figp(2) = max(1,figp(2)-figp(4));
      set(gcf,'Position',figp);
      plot(vv,'w');
      title('Loss function from independent observation');
      xlabel('index of (p,q) sorted in ascending order of p+q');
   end
end
% a spocteme minima
[vfpe ,k(1)] = min(vfpe(1:ke)); % p,q kde FPE nabyva minima
[vaicm,k(2)] = min(vaic(1:ke)); % p,q kde AIC nabyva minima
[vmdlm,k(3)] = min(vmdl(1:ke)); % p,q kde MDL nabyva minima
if ~isempty(vv)
  [vvm,k(4)] = min(vv)  ; % p,q kde vv  nabyva minima
end
pqo = pq(:,k);
