%========================================================================
%         BOX-JENKINSOVY ARMA MODELY PRO STACIONARNI CASOVE RADY
% Systematicky postup identifikace, odhadu parametru a overeni modelu AR
%                   (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 modelu AR
%                   (vyzaduje NAG FOUNDATION TOOLBOX)
%========================================================================
% V dalsim jiz budeme take vyuzivat funkce z NAG FOUNDATION TOOLBOXu (NAG).
% Zjisteni pritomnosti NAG FOUNDATION TOOLBOXu:
if ~exist('nagg')==3
   disp('CHYBA: NAG FOUNDATION TOOLBOX NENALEZEN');
   return;
end
% Libovolny stabilni kauzalni gaussovsky proces AR(p) vygenerujeme pomoci
% procedury 'armagen'.
% Na rozdil od pouze orientacniho pristupu v 'crcv11' provedeme nyni
% exaktnejsi analyzu s vyuzitim nekterych funkci, ktere nabizi
% NAG FOUNDATION TOOLBOX.
% 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 AR(2)
%=========================================================
% 1. GENEROVANI PROCESU
%    ------------------
  p=2; n=300;
% Rozptyl, resp. smerodatna odchylka generujiciho gaussovskeho bileho sumu
  sigma2=2.25; sigma=sqrt(sigma2);
% Generovani AR(2) procesu s parametry:
  PHI = [0.5;0.2];
  [x,PHIR,THETAR]=armagen(PHI,[],2*n,sigma);
  PHIR, abs(PHIR)
% Koreny polynomu PHI(z) lezi vne jednotkoveho kruhu,
% takze proces je kauzalni
  xv = x(n+1:2*n);  % 2. usek dat pouzijeme pro overeni modelu
  x  = 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 STREDNI HODNOTY A AUTOKOVARIANCNI FUNKCE PROCESU (procedura acf)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [xro,xroE,mi,miE,xacf]=acf(x,figure);
  mi                      % Stredni hodnota a jeji konfidencni intervaly
  [mi-miE(1),mi+miE(1)]   % pouze za predpokladu stacionarity
  [mi-miE(2),mi+miE(2)]   % presnejsi za predpokladu stacionarity a normality
% - Nula by mela lezet v intervalu spolehlivosti s polosirkou malou relativne
%   k velikosti dat, kdy nezamitame hypotezu o nulove stredni hodnote.
% - Autokorelacni funkce by mela klesat exponencialne k nule, pricemz
%   v konfidencnim pasu nulovosti by nemela mit nahodny charakter.
%   Proto usuzujeme, ze se nejedna o MA proces.
pause;
%-----------------------------------------------------------------------------
% A nyni alternativne s vyuzitim procedury 'NAGacf':
  nq = length(xro)-1;
  [NAGxro,NAGxv,NAGmi,stat] = NAGacf(x,nq);
% vidime, ze ziskane vysledky jsou shodne v ramci strojove presnosti:
  max(abs(xro(2:nq+1)-NAGxro))
  mi-NAGmi
pause;
%-----------------------------------------------------------------------------
% Navic statistika 'stat' ma priblizne rozdeleni 'chi2(nq)' a umoznuje testovat
% hypotezu o nulovosti NAGxro (tj. x je bily sum) s rizikem zamitnuti alpha:
  alpha = 0.05;
  kvantil=chi2inv(1-alpha,nq), stat
% hypotezu o nulovosti NAGxro bychom meli zamitnout s rizikem alpha:
  stat>kvantil
pause;
%-----------------------------------------------------------------------------
% b) ODHAD PARCIALNI AUTOKORELACNI FUNKCE A AUTOREGRESNICH KOEFICIENTU phi
%    (procedura pacf)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [xpacf,phi]=pacf(xacf,n,figure);
% Pro h>2 by pacf(h) mela lezet v nulovem konfidencnim pasu a vykazovat
% zde nahodne chovani.
% Podle ocekavani tedy usuzujeme na AR proces radu p=2.
pause;
%-----------------------------------------------------------------------------
% POZNAMKA: Pokud nepozadujeme vystupni parametr v, lze do pacf dosadit
%           misto autokovariancni funkce xacf autokorelacni funkci xro.
  [xpacf0,phi0,ans,v]=pacf(xro,n);
  max(abs(xpacf-xpacf0))
  max(abs(phi-phi0))
% Obdrzene vysledky jsou stejne.
pause;
%-----------------------------------------------------------------------------
% Porovname vysledky jeste s vystupy procedury 'NAGpacf':
  [NAGxpacf,NAGphi,vr]=NAGpacf(NAGxro,nq);
  max(abs(NAGxpacf-xpacf0(2:nq+1)))
  max(abs(NAGphi-phi0))
  vr(nq)-v
% Tedy posledni hodnota vr(nq)=(E(Xest(nq+1)-X(nq+1))^2)/var(X(1)) odpovida v,
% jestlize na vstup 'pacf' dame 'xro'.
pause;
%-----------------------------------------------------------------------------
close all
% c) ODHAD RADU p PROCESU AR(p)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~
% Spocteme odhady rozptylu sigma2 generujiciho bileho sumu pro sirsi rozsah
% radu, napr. pro p=1:10:
% v(p)  ... odhad rozptylu spocteny z tehoz useku x uziteho pro odhad modelu
% vv(p) ... odhad rozptylu spocteny z nezavisleho useku xv
  pmax = 10;
  xpacf=[]; phi=[]; vp=[];
  for p=1:pmax
     [xpacf,phi,phiE,v(p)] = pacf(xacf(1:p+1),n,[],xpacf,phi,vp);
     [xvs,vz,vv(p)] = armasim(xv,phi,[],x);
     vp = v(p);
echo off
  end
echo on
pause;
%-----------------------------------------------------------------------------
% Totez lze spocist i uzitim NAGpacf:
  [ans,ans,vr]=NAGpacf(NAGxro,pmax);
  max(abs(vr*xacf(1)-v'))
% a vykreslime prubeh v(p) ... tzv. ztratove funkce (loss function)
% Pozn.: misto v lze uzit i vr, nebot na multiplikativni konstante nezalezi
  figure; plot(v,'w');
  title('Ztratova funkce');
  xlabel('rad p');
pause;
%-----------------------------------------------------------------------------
% Z teorie vyplyva, ze bychom meli hledat p, pro nez je v(p) minimalni,
% nebot se v L_2 provadi ortogonalni projekce nahodne veliciny X(t) do
% podprostoru generovaneho velicinami X(t-1),....,X(t-p) a sigma2 odpovida
% kvadratu vzdalenosti X(t) od tohoto podprostoru v L_2 norme.
% Zde vsak v(p) je pouze odhad sigma2 zjisteny z jedineho
% pozorovani x(t), t=1,...,n.
% Se zvetsujicim se poctem parametru p se ovsem v takovem pripade model
% snazi co nejlepe popsat pouze toto dane pozorovani ve smyslu linearni
% regrese pro vektory, tj. zahrnuje falesne do modelu i prislusnou realizaci
% generujiciho bileho sumu z(t) (tzv. prefitovani).
% Proto krivka v(p) ma vzdy klesajici charakter a spise nez minimum
% se v ni snazime hledat bod zlomu, v nemz prudky pokles se stava
% pozvolnejsim.
pause;
% Bod zlomu se pak snaze urci jako minimum vhodne kriterialni funkce
% vcf(p), ktera zohledni falesny pokles za bodem zlomu tim, ze
% ho umele navysi v zavislosti na velicine rp=celkovy pocet parametru
% v modelu pripadajici na jeden pozorovany bod (v nasem pripade rp=p/n).
% Doporucovane kriterialni funkce:
% vcf(p) = v(p)*(1+rp)/(1-rp) ... tzv. Final Prediction Error (FPE) kriterium
% vcf(p) = v(p)*(1+c*rp) pro vhodnou vahu c>=0.
% Nejcasteji pouzivane vahy:
% c = 2      ... tzv. AIC kriterium (za predpokladu normality Z(t))
%                (Akaiake's Information Theoretic Criterion)
%                Pozn.: rp male => AIC dobre aproximuje FPE.
% c = log(n) ... tzv. MDL kriterium (zde log znaci prirozeny logaritmus)
%                (Rissanen's Minimum Description Length Criterion)
pause;
%-----------------------------------------------------------------------------
% Z prubehu v(p) by se jako nejpravdepodobnejsi bod zlomu mel jevit p=2.
% Overeni provedeme prikreslenim prubehu FPE, AIC a MDL:
  rp = (1:pmax)./n; AIC = 2; MDL = log(n);
  vfpe = v.*(1+rp)./(1-rp);
  vaic = v.*(1+ AIC*rp);
  vmdl = v.*(1+ MDL*rp);
  hold on
  plot(vfpe,'r--'); plot(vaic,'g-.'); plot(vmdl,'g:');
  title('Ztratova funkce (-), FPE (r--), AIC (g-.) a MDL (g:)');
  hold off;
  figure(gcf);
% Vidime, ze zejmena MDL kriterium dobre indikuje minimum pro p=2.
pause;
%-----------------------------------------------------------------------------
% Pokud pouzijeme odhad vv(p) spocteny z jineho nezavisleho pozorovani
% tehoz procesu, volime c=0, coz odpovida pripadu vcf(p)=vv(p).
  figure; plot(vv,'w');
  title('Ztratova funkce z nezavisleho pozorovani');
  xlabel('rad p');
% Meli bychom videt, ze skutecne vv(p) nabyva na rozdil od v(p) minima
% a to opet pro p=2, takze celkove budeme akceptovat model AR(2).
pause;
%-----------------------------------------------------------------------------
close all
% 3. ODHAD PARAMETRU PROCESU IDENTIFIKOVANEHO JAKO AR(2)
%    ---------------------------------------------------
% Zde budeme jiz prevazne vyuzivat funkce z NAG FOUNDATION TOOLBOXu,
% ktere umozni ziskat presnejsi odhad parametru, nez jaky dostavame z
% 'pacf', resp. z 'NAGpacf'.
%
% 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) Orientacni odhad pomoci 'pacf', resp. 'NAGpacf':
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  p = 2;
  [xpacf,phi,phiE,v]=pacf(xacf(1:p+1),n,figure);
  [phi,phiE]  % ORIENTACNI ODHAD PARAMETRU a jejich 95% konfidencni polosirky
  v*n/(n-1)   % ORIENTACNI ODHAD ROZPTYLU sigma2
% Vidime, ze ziskany odhad phi se blizi idealnim hodnotam PHI, ktere jsou
% zahrnuty v konfidencnim intervalu.
% Podobne v je dobrym odhadem sigma2.
  PHI, sigma2
pause;
%-----------------------------------------------------------------------------
% Obdobne vysledky (ale bez konfidencnich polosirek) dostaneme
% pomoci 'NAGpacf':
  [NAGxpacf,NAGphi,vr]=NAGpacf(NAGxro,p);
  NAGphi
  vr(p) - v/xacf(1)
pause;
%-----------------------------------------------------------------------------
% b) Predbezny odhad pomoci 'NAGpest':
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  mr = AR(p)     % vektor radu pro SARIMA model redukovany na AR(p).
  [par0,v0,isf] = NAGpest(mr,NAGxro,NAGxv);
% Predbezne odhady
  phi0 = par0{1}  % parametru phi
  v0              % rozptylu sigma2
  isf(1)          % hodnota 1 by mela indikovat uspesnost odhadu
% Vidime, ze v pripade AR modelu se predbezne odhady nelisi od orientacnich,
% takze v tomto pripade je pouziti 'pest' zbytecne (je totiz ekvivalentni
% s resenim Yule-Walkerova systemu linearnich rovnic tak, jak jej provadi
% pacf' nebo 'NAGpacf').
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}
  v0
  isf(1)
pause;
%-----------------------------------------------------------------------------
% c) Finalni odhad pomoci 'NAGest':
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [par1,par1S,v1,est1INFO] = NAGest(mr,x,par0);
% Takto jsme spocetli finalni odhady pomoci metody nejmensich ctvercu
% z NAG procedury 'G13AFF':
  phi1  = par1{1}  % parametru phi
  phi1S = par1S{1} % jejich smerodatnych odchylek
  v1               % rozptylu sigma2
% Z phi1S spocteme 95% konfidencni polosirky s vyuzitim
% normalni aproximace chyb odhadu parametru.
% Obecne: phi1E=norminv(1-alpha/2)*phi1S, 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%
% Srovnejme s polosirkami z orientacniho odhadu:
  phiE.'
pause;
%-----------------------------------------------------------------------------
% Prohledneme jeste obsah struktury estINFO, ktera obsahuje dodatecne
% informace:
  est1INFO
% kde je:
         % xdmi        ... odhad stredni hodnoty
         % xdmiS       ... jeji smerodatna odchylka
  est1INFO.corrmx   %  ... korelacni matice odhadu vsech parametru
                    %      (celkem tri parametry vcetne 'xdmi')
  est1INFO.isf(1)   %  ... 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 'xmi', resp. potlacit odhadovani 'xmi' zadanim xdmiEST = 0:
  xmi0 = mean(x);
  xmiEST = 0;
  [par2,par2S,v2,est2INFO] = NAGest(mr,x,par0,[],xmi0,xmiEST);
% Finalni odhady
  phi2  = par2{1}       % parametru phi
  phi2E = 1.96*par2S{1} % jejich konfidencnich polosirek
  v2                    % rozptylu sigma2
  est2INFO
pause;
%-----------------------------------------------------------------------------
% Definitivni odhad spocteme metodou maximalni verohodnosti z NAG procedury
% 'G13BEF':
  [par3,par3S,v3,est3INFO] = NAGest(mr,x,par0,'MLE');
% Finalni odhady
  phi3  = par3{1}       % parametru phi
  phi3E = 1.96*par3S{1} % jejich konfidencnich polosirek
  v3                    % rozptylu sigma2
  est3INFO
  xmi3 = est3INFO.xdmi; % odhad stredni hodnoty
pause;
%-----------------------------------------------------------------------------
%
% 4. OVERENI MODELU ODHADNUTEHO JAKO AR(2)
%    -------------------------------------
% a) Nejprve zkontrolujeme, zda hodnoty koncovych parametru jsou
%    vyznamne (jejich konfidencni intervaly neobsahuji nulu):
  phi3 = phi3'; phi3E = phi3E';
  [phi3-phi3E,phi3+phi3E]
% Pokud by napriklad [phi3(p)-phi3E(p),phi3(p)+phi3E(p)] obsahoval nulu,
% znamenalo by to, ze phi3(p) se nelisi vyznamne od nuly a ze tedy
% rad modelu je pravdepodobne nadhodnocen.
% Cely proces odhadu bychom pak zopakovali pro rad p-1 a znovu overili.
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 'x', z nehoz byl nalezen odhad
%    parametru modelu:
  [xs,z,vz] = armasim(x-xmi3,phi3);
  xs = xs + xmi3;                    % korekce o stredni hodnotu
% Porovname s estINFO.z:
  max(abs(z-est3INFO.z))
% Nebo spocteme 'z' jinak uzitim standardni funkce MATLABu 'filter':
  zz = filter([1;-phi3],1,x-xmi3); xxs = x-zz; vzz = (norm(z).^2)./n;
  max(abs(z-zz)), max(abs(xs-xxs)), vz-vzz
% Provedte help filter.
keyboard;
%-----------------------------------------------------------------------------
% Nyni pro verifikaci uzijeme nezavisly usek 'xv':
  [xvs,zv,vzv] = armasim(xv-xmi3,phi3,[],x-xmi3);
  xvs = xvs + xmi3;
  zzv = filter([1;-phi3],1,[x(n+1-p:n);xv]-xmi3); zzv = zzv(p+1:p+n);
  xxvs = xv-zzv; vzzv = (norm(zzv).^2)./n;
  max(abs(zv-zzv)), max(abs(xvs-xxvs)), vzv-vzzv
pause;
% Proverime, vlastnost IID rezidulaniho sumu napr. pomoci 'trand'
  trand(z)
  trand(zv)
% 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 ARMA(p,q) s vyssim radem p nrbo q.
% Protoze AR(p) = ARMA(p,0), museli bychom v nasem pripade uvazovat
% o modelu napr. AR(p+1) nebo o smisenem modelu pro p,q>0.
pause;
%-----------------------------------------------------------------------------
% Spocteme jeste autokorelacni funkci rezidualniho sumu z, resp. zv pomoci
% 'acf' nebo 'NAGacf' a proverime, zda se jedna o bily sum:
  acf(z ,figure);
  acf(zv,figure);
% V obou pripadech by skoro vsechny hodnoty pro h>0 mely byt v 95%
% konfidencnim pasu nulovosti, takze hypotezu o bilem sumu nezamitame.
pause;
%-----------------------------------------------------------------------------
% Podobne uzitim 'NAGacf':
  nqz = round(sqrt(n));
  [zro,ans,ans,stat] = NAGacf(z,nqz);
  stat, kvantil=chi2inv(0.95,nqz-p)
% tato relace by nemela platit, pokud nezamitame hypotezu o bilem sumu.
  stat>kvantil
pause;
close all;
%-----------------------------------------------------------------------------
% Pro vypocet rezidualniho sumu 'z' nabizi NAG funkci 'NAGcheck'
% ktera navic spocte autokorelacni funkci pro 'z' vcetne presnejsi
% portmanteau statistiky pro testovani, zda z je bily sum.
% Provedte help 'NAGcheck'.
keyboard;
%-----------------------------------------------------------------------------
  [zro3,zrocm3,stat3,idf3,siglev3] = NAGcheck(mr,par3,z);
  max(abs(zro3-zro))
  idf3, siglev3
  stat3, kvantil=chi2inv(0.95,idf3)
% tato relace by nemela platit, pokud nezamitame hypotezu o bilem sumu.
  stat3>kvantil
% nebo alternativne:
% tato relace by opet nemela platit, pokud nezamitame hypotezu o bilem sumu.
  0.05>siglev3
% spocteme jeste polosirky intervalu spolehlivosti pro 'zro3' pri normalni
% aproximaci rozlozeni odhadu:
  zro3E = 1.96*diag(zrocm3);
  errorbar(1:length(zro3),zro3,zro3E,zro3E,'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 'xs' a v 'xvs'.
%    V NAG predikce pocita funkce 'NAGpred'.
%    Provedte help 'NAGpred'.
keyboard;
%-----------------------------------------------------------------------------
% Porovname 1-krokovou predikci 'xs' se skutecnym pozorovnim puvodnich dat x:
   plot(x,'y'); hold on; plot(xs,'m'); hold off;
   title('Pozorovana data x (zluta), predikce (fialova)');
   zoom on
% Detaily lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
% Vsimnete si, ze graf predikce je jakoby o 1 krok posunuty vpravo, ale
% jinak kopiruje tvarove puvodni prubeh - to indikuje dobre odhadnuty
% model. Posunuti o 1 krok je dano tim, ze xs(t) ma nejvice informaci o
% datech do casu t-1. Jinak ovsem v praxi xs(t) pouzivame jako predikci
% x(t), i kdyz neni tak presna (x(t) totiz jeste nezname).
pause;
%-----------------------------------------------------------------------------
% Porovname 1-krokovou predikci 'xvs' se skutecnym pozorovnim nasledujiciho
% nezavisleho useku xv puvodnich dat:
   plot(xv,'y'); hold on; plot(xvs,'m'); hold off;
   title('Nezavisly usek pozorovanych dat (zluta), predikce (fialova)');
   zoom on
% Detaily opet lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
pause;
%-----------------------------------------------------------------------------
  zoom off
% Spocteme nyni predikci o 10 kroku pomoci 'NAGpred' a porovname ji s xv(1:10):
  npv = 10;   % delka predikce
  [xvs3,xvs3S] = NAGpred(mr,par3,npv,est3INFO);
% polosirka intervalu spolehlivosti pri normalni aproximaci chyb:
  xvs3E = 1.96*xvs3S;
  plot(xv(1:npv),'y'); hold on;
  plot(xvs3,'m'); plot([xvs3-xvs3E,xvs3+xvs3E],'m--'); hold off;
  title('Skutecny usek pozorovanych dat (zluta), jeho predikce (fialova)');
pause;
%-----------------------------------------------------------------------------
% Predikci (bez konfidencnich polosirek) muzeme spocist tez zadanim
% nulove smerodatne odchylky bileho sumu do 'armagen' :
  [xp3,phi3r] = armagen(phi3,[],npv,0,x-xmi3);
  xp3 = xp3 + xmi3;     % korekce o odhadnutou stredni hodnotu
% Porovnani s xvs3 dava shodny vysledek:
  max(abs(xvs3-xp3))
% Porovnejme jeste koreny idealniho polynomu PHI(z) a jeho odhadu phi(z):
  PHIR, phi3r
% Rovnez odhadnute koreny lezi vne jednotkoveho kruhu, coz garantuje
% kauzalitu procesu.
pause;
%-----------------------------------------------------------------------------
% Kvalita odhadu modelu se nejvice pozna pri predikci na vetsi vzdalenost.
% Porovname idealni predikci s odhadnutou:
  xp = armagen(PHI,[],npv,0,x); % IDEALNI PREDIKCE
% Do grafu zarazujeme i p predchozich hodnot, pomoci nichz se
% predikce spocetla
  t=(n-p+1:n+npv).';
  figure;plot(t,[x(n-p+1:n);xp],'w',t,[x(n-p+1:n);xp3],'y--')
  title('Predikce: idealni (-), odhadnuta (--)');
%==================================KONEC======================================
echo off
