%========================================================================
%         BOX-JENKINSOVY ARMA MODELY PRO STACIONARNI CASOVE RADY
% Systematicky postup identifikace, odhadu parametru a overeni smiseneho
%                          modelu ARIMA(p,d,q)
%                  (vyzaduje NAG FOUNDATION TOOLBOX)
%========================================================================

% (C) 1999, V. Vesely, Masaryk University of Brno (Czech Republic)

format compact;
echo on
%========================================================================
%         BOX-JENKINSOVY ARMA MODELY PRO STACIONARNI CASOVE RADY
% Systematicky postup identifikace, odhadu parametru a overeni smiseneho
%                          modelu ARIMA(p,d,q)
%                  (vyzaduje NAG FOUNDATION TOOLBOX)
%========================================================================
%
% Zjisteni pritomnosti NAG FOUNDATION TOOLBOXu:
if ~exist('nagg')==3
   disp('CHYBA: NAG FOUNDATION TOOLBOX NENALEZEN');
   return;
end
% Libovolny gaussovsky proces ARIMA(p,d,q) vygenerujeme pomoci procedury
% 'armagen' tak, aby 1 byla d-nasobnym korenem PHI(z).
% Knihovna NAG nabizi funkce 'NAGacf' a 'NAGpacf' jako alternativu k
% funkcim 'acf' a 'pacf' z knihovny 'cas_rady', coz nam umozni porovnat
% vysledky.
% Z knihovny 'cas_rady' dale vyuzijeme take procedury 'armagen' a 'armasim'.
% Provedte help armagen, help armasim, help acf, help NAGacf,
%          help pacf, help NAGpacf (navod k pouziti),
% resp.
%          type armagen, type armasim, type acf, type NAGacf,
%          type pacf, type NAGpacf (algoritmus).
keyboard;
%------------------------------------------------------------
% DEMONSTRACE ODHADU PARAMETRU SIMULOVANEHO PROCESU ARMA(2,2)
%============================================================
% 1. GENEROVANI PROCESU
%    ------------------
  p=2; q=2; d=1; n=300;
% Rozptyl, resp. smerodatna odchylka generujiciho gaussovskeho bileho sumu
  sigma2=5e-8; sigma=sqrt(sigma2);
% Generovani ARIMA(2,1,2) procesu:
  PHI = [0.5;0.2]; THETA = [-0.6;0.3]; % Parametry procesu ARMA(p,q)
  phir = roots([-flipud(PHI);1]) % Koreny polynomu PHI(z) procesu ARMA(p,q)
% Pridanim d-nasobneho korene 1 ziskame koeficienty polynomu PHID(z) procesu
% ARIMA(p,d,q):
  pol=poly([ones(d,1)+eps;phir]); np = length(pol);
% a z nej AR parametry procesu ARIMA(p,d,q)
  PHID = [-fliplr(pol(1:np-1)./pol(np))].';
  [x,phir,thetar]=armagen(PHID,THETA,2*n,sigma);
  phir.', abs(phir).'
% Dva koreny polynomu PHI(z) lezi vne jednotkoveho kruhu, a zbyvajici
% je d-nasobny koren 1 porusujici kauzalitu (=stacionaritu) procesu
  thetar.', abs(thetar).'
% Koreny polynomu THETA(z) lezi vne jednotkoveho kruhu,
% takze proces ARIMA(p,d,q) je invertibilni
pause;
%-----------------------------------------------------------------------------
  xv = x(n+1:2*n);  % 2. usek dat pouzijeme pro overeni modelu
  xm = x(1:n);      % 1. usek dat pouzijeme pro odhad modelu
% Proverime nahodnost
  xrnd = trand(x-mean(x))
% Alespon jeden z testu by mel zamitnout nahodnost (vraci 0).
% Pokud tomu tak neni, je x bily sum a nema smysl pokracovat.
% V takovem pripade spustte davku znovu.
if xrnd
   disp('PREDCASNY KONEC: Vygenerovana data se jevi jako nahodny sum');
%==================================KONEC======================================
echo off
return;
end
pause;
%-----------------------------------------------------------------------------
% 2. IDENTIFIKACE PROCESU
%    --------------------
% a) ODHAD RADU DIFERENCOVANI, STREDNI HODNOTY A AUTOKOVARIANCNI FUNKCE
%    DIFERENCOVANEHO PROCESU (procedura acf)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Pro diferencovani pouzijeme proceduru 'NAGdiff' z knihovny NAG.
% Provedte 'help NAGdiff' (napoveda), resp. 'type NAGdiff' (algoritmus)
keyboard;
%-----------------------------------------------------------------------------
  [xro,xroE,mi,miE,xacf]=acf(x,figure); xacf(1) % acf + rozptyl
% Rada vykazuje (nahodne) promenny trend a jeji autokorelacni funkce
% vykazuje pomaly (priblizne linearni) pokles. Oboji prokazuje nestacionaritu.
%
% Budeme proto zvysovat rad diferencovani, dokud nebude diferencovana rada
% vykazovat stacionarni chovani.
% Pro optimalni rad diferencovani (v praxi nebyva vetsi nez 2) by melo byt
% dosazeno nejmensi hodnoty rozptylu xacf(1):
pause;
%-----------------------------------------------------------------------------
  [xdro,xdroE,mid,midE,xdacf]=acf(NAGdiff(1,x),figure); xdacf(1)
% Rozptyl vyrazne poklesl, pricemz autokorelacni funkce jiz vykazuje
% exponencialni pokles. To jiz signalizuje stacionarni proces.
% Pro jistotu zkusime jeste pokracovat v diferencovani:
pause;
%-----------------------------------------------------------------------------
  [xdro,xdroE,mid,midE,xdacf]=acf(NAGdiff(2,x),figure); xdacf(1)
% Rozptyl jiz vzrostl a obcas i 'acf' signalizuje mozne poruseni stacionarity.
% Doslo tedy asi jiz k 'prediferencovani'.
% Zvolime proto v souladu s ocekavanim rad d=1:
pause;
%-----------------------------------------------------------------------------
  close all;
  d = 1;
  [xd,xr] = NAGdiff(d,x);     % Optimalne diferencovana rada
  nd = length(xd);
% Provereni konstantnosti rozptylu:
  [lambdaE,y] = powtr(xd+1.5*max(0,-min(xd)),[],figure);
pause;
%-----------------------------------------------------------------------------
% Konstantni rozptyl je potvrzen. Kdyby tomu tak nebylo, bylo by nutne
% provest 'powtr' s puvodnimi daty a misto x dale pracovat s odpovidajicim
% zpusobem transformovanou radou.
%
% Stredni hodnota diferencovane rady a jeji konfidencni intervaly
  [xdro,xdroE,mid,midE,xdacf]=acf(xd,figure);
  [mid-midE(1),mid+midE(1)]   % pouze za predpokladu stacionarity
  [mid-midE(2),mid+midE(2)]   % presnejsi za predpokladu stacionarity a normality
% Autokorelacni funkce by mela klesat exponencialne k nule, pricemz
% v konfidencnim pasu nulovosti by nemela mit nahodny charakter.
% Proto usuzujeme, ze diferencovana rada neni MA proces.
pause;
%-----------------------------------------------------------------------------
% A nyni alternativne s vyuzitim procedury 'NAGacf':
  nq = length(xdro)-1;
  [NAGxdro,NAGxdv,NAGmid,statd] = NAGacf(xd,nq);
% vidime, ze ziskane vysledky jsou shodne v ramci strojove presnosti:
  max(abs(xdro(2:nq+1)-NAGxdro))
  mid-NAGmid
pause;
%-----------------------------------------------------------------------------
% Navic statistika 'statd' ma priblizne rozdeleni 'chi2(nq)' a umoznuje testovat
% hypotezu o nulovosti NAGxdro (tj. xd je bily sum) s rizikem zamitnuti alpha:
  alpha = 0.05;
  kvantil=chi2inv(1-alpha,nq), statd
% hypotezu o nulovosti NAGxro bychom meli zamitnout s rizikem alpha:
  statd>kvantil
pause;
%-----------------------------------------------------------------------------
% b) ODHAD PARCIALNI AUTOKORELACNI FUNKCE A AUTOREGRESNICH KOEFICIENTU phi
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%    (procedura pacf)
%    ~~~~~~~~~~~~~~~~
  [xdpacf,phi]=pacf(xdacf,nd,figure);
% Obalka parcialni autokorelacni funkce by mela klesat exponencialne k nule,
% pricemz v konfidencnim pasu nulovosti by nemela mit nahodny charakter.
% Proto predbezne usuzujeme, ze diferencovana rada neni AR proces konecneho
% radu.
pause;
%-----------------------------------------------------------------------------
% Porovname vysledky jeste s vystupy procedury 'NAGpacf':
  [NAGxdpacf,NAGphi,vr]=NAGpacf(NAGxdro,nq);
  max(abs(NAGxdpacf-xdpacf(2:nq+1)))
  max(abs(NAGphi-phi))
pause;
%-----------------------------------------------------------------------------
close all
% c) POKUS O ODHAD RADU p PRO MODEL AR(p)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Odhad v(p) z useku xd a vv(p) z useku xdv pro kazdy rad p=1,...,pmax a
% nalezeni optimalnich radu pro nektera kriteria spocteme pomoci 'NAGpqest'.
% v(p)  ... odhad rozptylu spocteny z tehoz useku xd uziteho pro odhad modelu
% vv(p) ... odhad rozptylu spocteny z nezavisleho useku xdv
% v(p)=NaN a vv(p)=NaN, pokud nemohly byt ziskany spolehlive odhady.
% To zpravidla indikuje preparametrizovani (prilis velky rad p).
%
% Provedte help NAGpqest (navod k pouziti), resp.
%          type NAGpqest (algoritmus).
%
keyboard;
%-----------------------------------------------------------------------------
  pmax = 10;
  xdv = xd(n+1:length(xd));   % 2. usek pouzijeme pro overeni modelu
  xd  = xd(1:n);              % 1. usek pouzijeme pro odhad modelu
  [pqo,v,vv,ke,pq,rp] = NAGpqest(AR(pmax),xd,xdv,figure);
  pqo(1,:)   % optimalni rady p spoctene porade z FPE,AIC,MDL a vv
  pq(:,1:ke) % rady p,q, z nichz mohl byt uspesne proveden optimalni vyber
             % v(ke+1) = vv(ke+1) = NaN, pokud ke<pmax.
pause;
  rp  % pocty parametru pripadajici na jedno pozorovani:
      % umoznuje spocist jine kriterium, napr. v.*(1+c*rp) pro vhodne c.
  v10 = v.*(1+10*rp)  % pro c=10.
  figure; plot(v10);
  title('Silne modifikujici kriterium (c=10)');
  xlabel('rad p');
  [v10min,kmin] = min(v10)   % kmin=po ... optimalni rad p dle kriteria v10

pause;
%-----------------------------------------------------------------------------
% Body zlomu vetsiny ztratovych funkci, resp. minima kriterialnich funkci
% by mely byt nevyrazna a davat spise vyssi navzajem se lisici hodnoty.
% Proto se na zaklade prubehu 'acf' a 'pacf' ve 2a)b) priklonime
% k neadkvatnosti AR modelu a pokusime se najit spise optimalni ARMA model.
pause;
%-----------------------------------------------------------------------------
close all
% d) ODHAD RADU p,q DIFERENCOVANEHO PROCESU ARMA(p,q), p>=0, q>0
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Odhad v(q) z useku xd a vv(q) z useku xdv pro kazdy rad
% p=0,1,...pmax; q=1,...,qmax a nalezeni optimalnich radu pro nektera
% kriteria spocteme opet pomoci 'NAGpqest', k->(p,q).
% v(k)  ... odhad rozptylu spocteny z tehoz useku xd uziteho pro odhad modelu
% vv(k) ... odhad rozptylu spocteny z nezavisleho useku xdv
% v(k)=NaN a vv(k)=NaN, pokud nemohly byt ziskany spolehlive odhady (errc=7).
% To zpravidla indikuje preparametrizovani (prilis velke p+q).
  pmax=4; qmax=4;
  [pqo,v,vv,ke,pq,rp] = NAGpqest(ARMA(pmax,qmax),xd,xdv,figure);
  pqo        % optimalni rady p,q spoctene porade z FPE,AIC,MDL a vv
  pq(:,1:ke) % rady p,q, z nichz mohl byt uspesne proveden optimalni vyber
             % v(ke+1) = vv(ke+1) = NaN, pokud ke<(pmax+1)*qmax.
pause;
  rp  % pocty parametru pripadajici na jedno pozorovani:
      % umoznuje spocist jine kriterium, napr. v.*(1+c*rp) pro vhodne c.
  v10 = v.*(1+10*rp)  % pro c=10.
  figure; plot(v10);
  title('Silne modifikujici kriterium (c=10)');
  xlabel('rady p,q');
  [v10min,kmin] = min(v10)
  pq(:,kmin) % optimalni rady p,q dle kriteria v10
%
% Grafy i minima kriterialnich funkci by mely urcovat bod zlomu a tedy
% rady p=2,q=2.
pause;
%-----------------------------------------------------------------------------
close all
% 3. ODHAD PARAMETRU PROCESU IDENTIFIKOVANEHO JAKO ARIMA(2,1,2)
%    ----------------------------------------------------------
% Pro odhad parametru nabizi NAG funkce:
% NAGpest ... predbezny odhad parametru
% NAGest  ... finalni odhad parametru minimalizaci souctu ctvercu odchylek,
%             ktery je obecne nelinearni funkci neznamych parametru.
% Provedte help NAGpest a help NAGest.
keyboard;
%-----------------------------------------------------------------------------
% a) Predbezny odhad pomoci 'NAGpest':
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  p = 2; q = 2;
  mr = ARIMA(p,d,q) % vektor radu pro SARIMA model redukovany na ARIMA(p,d,q).
  [par0,v0,isf] = NAGpest(mr,NAGxdro,NAGxdv);
% Predbezne odhady
  phi0   = par0{1}  % parametru phi
  theta0 = par0{2}  % parametru theta
  v0                % rozptylu sigma2
  isf(1:2)          % hodnota 1 by mela indikovat uspesnost odhadu
pause;
%-----------------------------------------------------------------------------
% Autokorelacni funkci nemusime pocitat predem, pokud 'NAGpest' zavolame
% takto:
  [par0,v0,isf] = NAGpest(mr,x);
% Pri tomto volani zaridi vse procedura 'NAGpest' sama a obdrzene vysledky
% jsou opet stejne
  phi0   = par0{1}  % parametru phi
  theta0 = par0{2}  % parametru theta
  v0                % rozptylu sigma2
  isf(1:2)          % hodnota 1 by mela indikovat uspesnost odhadu
pause;
%-----------------------------------------------------------------------------
% b) Finalni odhad pomoci 'NAGest':
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [par1,par1S,v1,est1INFO] = NAGest(mr,xm,par0);
% Takto jsme spocetli finalni odhady pomoci metody nejmensich ctvercu
% z NAG procedury 'G13AFF':
  phi1    = par1{1}  % parametru phi
  theta1  = par1{2}  % parametru theta
  phi1S   = par1S{1} % jejich smerodatnych odchylek
  theta1S = par1S{2}
  v1                 % rozptylu sigma2
% Z phi1S a theta1S spocteme 95% konfidencni polosirky s vyuzitim
% normalni aproximace chyb odhadu parametru.
% Obecne: phi1E=norminv(1-alpha/2)*phi1S, resp.
%         theta1E=norminv(1-alpha/2)*theta1S, kde 'norminv' je kvantilova funkce
%         normalniho rozlozeni a alpha riziko.
% Odtud:
  phi1E   = 1.96*phi1S   % 1.96 = norminv(1-0.05/2) odpovida riziku 5%
  theta1E = 1.96*theta1S
pause;
%-----------------------------------------------------------------------------
% Prohledneme jeste obsah struktury estINFO, ktera obsahuje dodatecne
% informace:
  est1INFO
% kde je:
         % xdmi        ... odhad stredni hodnoty diferencovane rady
         % xdmiS       ... jeji smerodatna odchylka
  est1INFO.corrmx   %  ... korelacni matice odhadu vsech parametru
                    %      (celkem pet parametru vcetne 'xdmi')
  est1INFO.isf(1:2) %  ... hodnota 1 by mela indikovat uspesnost odhadu
         % stset       ... stavova mnozina pro dalsi zpracovani
         % z           ... rezidua (mela by reprezentovat pozorovani
         %                          generujiciho bileho sumu Z(t))
         % rss         ... rezidualni soucet ctvercu
         % ndf         ... pocet stupnu volnosti
         % itc         ... skutecny pocet iteraci provedenych pri minimalizaci
         %                 souctu ctvercu odchylek.
pause;
%-----------------------------------------------------------------------------
% Do NAGest lze pridat dalsi vstupni parametry, napriklad vlastni
% pocatecni odhad 'xdmi', resp. potlacit odhadovani 'xdmi' zadanim xdmiEST=0:
  xdmi0 = mean(xd);
  xdmiEST = 0;
  [par2,par2S,v2,est2INFO] = NAGest(mr,xm,par0,[],xdmi0,xdmiEST);
% Finalni odhady
  phi2    = par2{1}  % parametru phi
  theta2  = par2{2}  % parametru theta
  phi2S   = par2S{1} % jejich smerodatnych odchylek
  theta2S = par2S{2}
  v2                 % rozptylu sigma2
  est2INFO
pause;
%-----------------------------------------------------------------------------
% Definitivni odhad spocteme metodou maximalni verohodnosti z NAG procedury
% 'G13BEF':
  [par3,par3S,v3,est3INFO] = NAGest(mr,xm,par0,'MLE');
% Finalni odhady
  phi3    = par3{1}       % parametru phi
  theta3  = par3{2}       % parametru theta
  phi3E   = 1.96*par3S{1} % jejich konfidencnich polosirek
  theta3E = 1.96*par3S{2}
  v3                      % rozptylu sigma2
  est3INFO
  xdmi3 = est3INFO.xdmi   % odhad stredni hodnoty diferencovane rady
pause;
%-----------------------------------------------------------------------------
% 4. OVERENI MODELU ODHADNUTEHO JAKO ARIMA(2,1,2)
%    --------------------------------------------
% a) Nejprve zkontrolujeme, zda hodnoty koncovych parametru jsou
%    vyznamne (jejich konfidencni intervaly neobsahuji nulu):
  phi3 = phi3'; phi3E = phi3E';
  [phi3-phi3E,phi3+phi3E]
  theta3 = theta3'; theta3E = theta3E';
  [theta3-theta3E,theta3+theta3E]
% Pokud by napriklad [theta3(q)-theta3E(q),theta3(q)+theta3E(q)] obsahoval
% nulu, znamenalo by to, ze theta3(q) se nelisi vyznamne od nuly a ze tedy
% rad modelu je pravdepodobne nadhodnocen.
% Cely proces odhadu bychom pak zopakovali pro ARIMA(p,d,q-1).
% Podobne bychom presli na model ARIMA(p-1,d,q), pokud phi(p) by byl
% nevyznamny parametr.
pause;
%-----------------------------------------------------------------------------
% b) Proverime, zda odhadnuty model a jeho komponenty splnuji teoreticke
%    predpoklady, za nichz byl nalezen.
%    V nasem pripade staci spocist tzv. rezidualni sum, tj. zbytkovy
%    sum, ktery by mel korespondovat s generujicim bilym sumem.
%    Uzijeme pro verifikaci nejprve usek dat 'xd', z nehoz byl nalezen odhad
%    parametru modelu:
  [xds,zd,vzd] = armasim(xd-xdmi3,phi3,theta3);
  xds = xds + xdmi3;                % korekce o stredni hodnotu
% Porovname s estINFO.z:
  max(abs(zd-est3INFO.z))
  plot([zd,est3INFO.z]); zoom on;
% odhady se lisi nejvice na pocatku v dusledku rozdilnych pocatecnich
% podminek pro Z(t). Totiz 'armasim' dosazuje nulove, zatimco 'NAGest'
% spocte Z(t) z prislusne invertibilni reprezentace.
pause;
  zoom off;
%-----------------------------------------------------------------------------
% Nyni pro verifikaci uzijeme nezavisly usek 'xdv':
  [xxdvs,zzdv,vzzdv] = armasim(xdv-xdmi3,phi3,theta3);
  xxdvs = xxdvs+xdmi3;
% Jeste jednou presneji s doplnenim pocatecnich podminek pro 'xd' i 'zd':
  [xdvs,zdv,vzdv] = armasim(xdv-xdmi3,phi3,theta3,xd-xdmi3,zd);
  xdvs = xdvs + xdmi3;
% Vysledky se nepatrne zmenily:
  max(abs(zdv-zzdv)), vzdv-vzzdv
  plot([xdvs,xxdvs]); zoom on;
pause;
% Proverime vlastnost IID rezidulaniho sumu napr. pomoci 'trand'
  zd = est3INFO.z;
  xds = xd - zd;
  trand(zd)
  trand(zdv)
% V nasem pripade by zadny test nemel zamitnout nahodnost.
% Pokud by v realnem pripade k zamitnuti nekterymi testy doslo,
% znacilo by to, ze v rezidualnim sumu existuji korelace a tedy
% je nutno prejit pravdepodobne k modelu ARIMA(p,d,q) s vyssim radem p nebo q.
pause;
%-----------------------------------------------------------------------------
% Spocteme jeste autokorelacni funkci rezidualniho sumu 'zd', resp. 'zdv' pomoci
% 'acf' nebo 'NAGacf' a proverime, zda se jedna o bily sum:
  acf(zd ,figure);
  acf(zdv,figure);
% V obou pripadech by vsechny hodnoty pro h>0 mely byt v 95% konfidencnim
% pasu nulovosti, takze hypotezu o bilem sumu nezamitneme.
pause;
%-----------------------------------------------------------------------------
% Podobne uzitim 'NAGacf':
  nqz = round(sqrt(n));
  [zdro,ans,ans,statd] = NAGacf(zd,nqz);
  statd, kvantil=chi2inv(0.95,nqz-(p+q))
% tato relace by nemela platit, pokud nezamitame hypotezu o bilem sumu.
  statd>kvantil
pause;
close all;
%-----------------------------------------------------------------------------
% Pro verifikaci rezidualniho sumu z nabizi NAG funkci 'NAGcheck'
% ktera navic spocte autokorelacni funkci pro 'zd' vcetne presnejsi
% portmanteau statistiky pro testovani, zda 'zd' je bily sum.
% Provedte help 'NAGcheck'.
keyboard;
%-----------------------------------------------------------------------------
  [zdro3,zdro3cm,statd3,idf3,siglevd3] = NAGcheck(mr,par3,zd);
  max(abs(zdro3-zdro))
  idf3, siglevd3
  statd3, kvantil=chi2inv(0.95,idf3)
% tato relace by nemela platit, pokud nezamitame hypotezu o bilem sumu.
  statd3>kvantil
% nebo alternativne:
% tato relace by opet nemela platit, pokud nezamitame hypotezu o bilem sumu.
  0.05>siglevd3
% spocteme jeste polosirky intervalu spolehlivosti pro 'zdro3' pri normalni
% aproximaci rozlozeni odhadu:
  zdro3E = 1.96*diag(zdro3cm);
  errorbar(1:length(zdro3),zdro3,zdro3E,zdro3E,'o');
  title('Rezidualni autokorelacni funkce');
pause;
%-----------------------------------------------------------------------------
% c) na zaver overime vizualne z grafu soulad pozorovanych hodnot
%    s hodnotami predikovanymi z p hodnot pozorovanych o 1 krok zpet
%    pri vyuziti odhadnuteho modelu.
%    Tuto predikci jiz drive spocetla funkce 'armasim' v 'xds' a v 'xdvs'.
%-----------------------------------------------------------------------------
% Porovname 1-krokovou predikci 'xds' s diferencovanymi daty 'xd':
   plot(xd,'y'); hold on; plot(xds,'m'); hold off;
   title('Diferencovana data xd (zluta), predikce (fialova)');
   zoom on
% Detaily lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
pause;
%-----------------------------------------------------------------------------
% Porovname 1-krokovou predikci 'xdvs' se skutecnym pozorovnim nasledujiciho
% nezavisleho useku xdv puvodnich dat:
   plot(xdv,'y'); hold on; plot(xdvs,'m'); hold off;
   title('Nezavisly usek diferencovanych dat (zluta), predikce (fialova)');
   zoom on
% Detaily opet lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
pause;
%-----------------------------------------------------------------------------
% Jeste jednou totez s tim, ze 1-krokovou predikci provedeme v puvodni casove
% rade. K tomu je nezbytne pretransformovat puvodni ARIMA model na
% ekvivalentni ARMA model pomoci funkce 'ARMApar':
% Provedte 'help ARMApar' (napoveda), resp. 'type ARMApar' (algoritmus).
% V NAG predikce pro puvodni nediferencovanou radu pocita funkce 'NAGpred'.
% Provedte 'help NAGpred'.
keyboard;
  [phi3x,theta3x,xc3]=ARMApar(mr,par3,xdmi3);
  [xp,z,vz] = armasim(xm,phi3x,theta3x); xp=xp+xc3;
% Porovname 1-krokovou predikci 'xp' s pozorovanymi daty 'xm':
   plot(xm,'y'); hold on; plot(xp,'m'); hold off;
   title('Pozorovana data x (zluta), predikce (fialova)');
   zoom on
% Detaily lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
pause;
%-----------------------------------------------------------------------------
% Porovname 1-krokovou predikci 'xpv' se skutecnym pozorovnim nasledujiciho
% nezavisleho useku xv puvodnich dat:
   [xpv,zv,vzv] = armasim(xv,phi3x,theta3x,xm); xpv=xpv+xc3;
   plot(xv,'y'); hold on; plot(xpv,'m'); hold off;
   title('Nezavisly pozorovany usek (zluta), predikce (fialova)');
   zoom on
% Detaily opet lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
pause;
%-----------------------------------------------------------------------------
% Spocteme na zaver predikci o 10 kroku pomoci 'NAGpred' a porovname ji
% s xv(1:10):
  npv = 10;   % delka predikce
  [xp3,xp3S] = NAGpred(mr,par3,npv,est3INFO);
% polosirka intervalu spolehlivosti pri normalni aproximaci chyb:
  xp3E = 1.96*xp3S;
  plot(xv(1:npv),'y'); hold on;
  plot(xp3,'m'); plot([xp3-xp3E,xp3+xp3E],'m--'); hold off;
  title('Skutecny usek pozorovanych dat (zluta), jeho predikce (fialova)');
pause;
%-----------------------------------------------------------------------------
% Kvalita odhadu modelu se nejvice pozna pri predikci na vetsi vzdalenost.
% Porovname idealni predikci s odhadnutou:
  estINFO = est3INFO; estINFO.xdmi = 0;
  xp = NAGpred(mr,{PHI,THETA,[],[]},npv,estINFO);  % IDEALNI PREDIKCE
% Do grafu zarazujeme i p predchozich hodnot, pomoci nichz se
% predikce spocetla
  t=(n-p+1:n+npv).';
  figure;plot(t,[xm(n-p+1:n);xp],'w',t,[xm(n-p+1:n);xp3],'y--')
  title('Predikce: idealni (-), odhadnuta (--)');
%==================================KONEC======================================
echo off
