%========================================================================
%         BOX-JENKINSOVY ARMA MODELY PRO STACIONARNI CASOVE RADY
% Systematicky postup identifikace, odhadu parametru a overeni modelu AR
%              (vyzaduje SYSTEM IDENTIFICATION TOOLBOX)
%========================================================================

% (C) 1996-98, V. Vesely, Masaryk University of Brno (Czech Republic)

format compact;
echo on
% 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
% SYSTEM IDENTIFICATION TOOLBOX.
% V pripade AR modelu pouzijeme take alternativne funkce 'pacf' a 'armasim'
% z knihovny 'cas_rady' a porovname vysledky.
% Provedte help armagen, help armasim, help acf, help pacf (navod k pouziti),
% resp.
%          type armagen, type armasim, type acf, type pacf (algoritmus)
keyboard;
%---------------------------------------------------------
% DEMONSTRACE ODHADU PARAMETRU SIMULOVANEHO PROCESU AR(2)
%=========================================================
% 1. GENEROVANI PROCESU
%    ------------------
  P=2; n=500;
% 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
  x  = x - mean(x); % Pro jistotu data vycentrujeme
  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)
% 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.
% - Obalka autokorelacni funkce klesa exponencialne k nule, pricemz
%   v konfidencnim pasu nulovosti nema nahodny charakter.
%   Proto usuzujeme, ze se nejedna o MA proces.
pause;
%-----------------------------------------------------------------------------
% b) ODHAD PARCIALNI AUTOKORELACNI FUNKCE A AUTOREGRESNICH KOEFICIENTU phi
%    (procedura pacf)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [xpacf,phi]=pacf(xacf,n,figure);
% Vidime, ze pro h>2 lezi pacf(h) v nulovem konfidencnim pasu a vykazuje
% 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]=pacf(xro,n,figure);
  max(abs(xpacf-xpacf0))
  max(abs(phi-phi0))
% Obdrzene vysledky jsou stejne.
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,z,vv(p)] = armasim(xv,phi,[],x(n+1-p:n));
     vp = v(p);
echo off
  end
echo on
% a vykreslime prubeh v(p) ... tzv. ztratove funkce (loss function)
  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 (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).
% Casto se pouziva nasledujici trida kriterialnich funkci:
% vcf(p) = v(p)*(1+c*rp) pro vhodnou vahu c>=0.
% Nejcasteji pouzivana kriteria:
% c = 2      ... tzv. AIC kriterium
%                (Akaiake's Information Theoretic Criterion)
% c = log(n) ... tzv. MDL kriterium (zde log znaci prirozeny logaritmus)
%                (Rissanen's Minimum Description Length Criterion)
pause;
%-----------------------------------------------------------------------------
% Z prubehu v(p) se jako nejpravdepodobnejsi bod zlomu jevi skutecne p=2.
% Overeni provedeme prikreslenim prubehu AIC a MDL:
  rp = (1:pmax)./n; AIC = 2; MDL = log(n);
  hold on
  vaic = v.*(1+ AIC*rp);
  vmdl = v.*(1+ MDL*rp);
  plot(vaic,'g--'); plot(vmdl,'g:');
  title('Ztratova funkce (-), AIC (--) a MDL (:)');
  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');
% Vidime, ze skutecne vv(p) nabyva na rozdil od v(p) minima a to opet pro p=2
% a muzeme tedy celkove akceptovat model AR(2).
pause;
%-----------------------------------------------------------------------------
% V dalsim jiz budeme prevazne vyuzivat funkce ze SYSTEM IDENTIFICATION
% TOOLBOXu (IDENT), ktere nabizi vice moznosti.
%
% Zjisteni pritomnosti SYSTEM IDENTIFICATION TOOLBOXu:
if ~exist('ident')==2
   disp('CHYBA: SYSTEM IDENTIFICATION TOOLBOX NENALEZEN');
   return;
end
% ODHAD RADU p PROCESU AR(p) pomoci IDENT:
% pro odhad radu p nabizi IDENT funkce:
% arxstruc ... vypocet ztratove funkce pro zadany rozsah radu
% selstruc ... nalezeni optimalniho radu modelu
% Provedte help arxstruc a help selstruc.
keyboard;
%-----------------------------------------------------------------------------
% Odhadneme rad uzitim techto funkci:
  V = arxstruc(x,x,(1:pmax).');    % Ztratova funkce jen z x
  V(1,1:pmax)
% koresponduje s
  v
  figure; plot(1:pmax,v,'w',1:pmax,V(1,1:pmax),'g');
pause;
  figure;
  p = selstruc(V) % graf ztratove funkce, v p nabyva minima
pause;
%-----------------------------------------------------------------------------
  [p,vmod] = selstruc(V,'aic')    % p, kde AIC nabyva minima
  p = selstruc(V,'mdl')  % p, kde MDL nabyva minima
  p = selstruc(V,3)      % p, kde kriterium s vahou c=3 nabyva minima
  vmod(1,:)
% koresponduje s
  log(vaic)
  figure; plot(1:pmax,log(vaic),'w',1:pmax,vmod(1,:),'g');
pause;
%-----------------------------------------------------------------------------
  VV = arxstruc(x,xv,(1:pmax).');  % Ztratova funkce uzitim xv
  VV(1,1:pmax)
% koresponduje s
  vv
  figure; plot(1:pmax,vv,'w',1:pmax,VV(1,1:pmax),'g');
pause;
  figure;
  p = selstruc(VV) % graf ztratove funkce, v p nabyva minima
pause;
%-----------------------------------------------------------------------------
  p = selstruc(VV,'aic') % p, kde AIC nabyva minima
  p = selstruc(VV,'mdl') % p, kde MDL nabyva minima
  p = selstruc(VV,3)     % p, kde kriterium s vahou c=3 nabyva minima
% Grafy i minima kriterialnich funkci dobre urcuji bod zlomu a tedy
% rad pro p=2
pause;
%-----------------------------------------------------------------------------
close all
% 3. ODHAD PARAMETRU PROCESU IDENTIFIKOVANEHO JAKO AR(2)
%    ---------------------------------------------------
  p = 2;
  [xpacf,phi,phiE,v]=pacf(xacf(1:p+1),n,gcf);
  [phi,phiE]  % ZAVERECNY ODHAD PARAMETRU a jejich konfidencni polosirky
  v           % ZAVERECNY 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 = 2.5.
PHI, sigma2
% IDENT nabizi pro odhad parametru AR procesu funkci 'ar'.
% Vysledky odhadu a dalsi informace o modelu jsou ulozeny ve specialnim
% formatu (tzv. theta format) do vystupni matice 'th' (viz help theta).
% Pro vetsi pohodli nabizi IDENT radu funkci pro manipulaci s theta
% formatem:
% - prezentace obsahu theta na displej: funkce present
% - extrakce informaci z  'th' (jmena funkci zacinaji th2...)
% - ukladani informaci do 'th' (jmena funkci konci ...th2)
% Provedte help ar, help present, help th2arx, th2par, th2poly
keyboard
%-----------------------------------------------------------------------------
% Treti parametr 'approach' funkce 'ar' dokonce dovoluje vyber
% z nekolika ruznych metod odhadu.
% Metoda 'yw' koresponduje s metodou pouzitou v 'pacf', ktera je zalozena
% na reseni Yule-Walkerova systemu rovnic Durbin-Levinsonovym algoritmem.
% Provedeme proto odhad nejprve touto metodou a porovname vysledky.
  th = ar(x,p,'yw');
  present(th);
  [A,B,dA,dB]=th2arx(th)
echo on
  phi1 = -A(2:p+1).'
% koresponduje s
  phi
% a v dA jsou smerodatne odchylky pro A, tj.
  phi1E=1.96*dA(2:p+1).'
% koresponduje s
  phiE
% Vidime pomerne dobrou shodu z obou odhadu.
pause;
%-----------------------------------------------------------------------------
% Misto 'th2arx' muzeme pro ziskani vysledku odhadu pouzit funkci
% 'th2par', ktera vraci v 'par' vsechny odhadovane parametry porade
% odpovidajici polynomum A,B,C,D,F a navic i jejich kovariancni matici P
% a rovnez odhad LAM rozptylu sigma2:
  [par,P,LAM]=th2par(th)
% kde
  phi1 = -par.'; phi1E = 1.96*sqrt(diag(P));
  [phi1,phi1E]
% koresponduje opet s
  [phi,phiE]
% a navic
  v1 = LAM
% koresponduje s
  v
pause;
%-----------------------------------------------------------------------------
% Podobne 'th2poly' vraci koeficienty vsech polynomu A,B,C,D,F
% (vcetne neodhadovaneho absolutniho clenu 1), LAM a navic i casovy
% krok diskretizace T, ktery jsme nespecifikovali, takze v nasem
% pripade je T=1:
  [A,B,C,D,F,LAM,T]=th2poly(th)
echo on
pause;
%-----------------------------------------------------------------------------
% Koreny (zeros) a poly ARMA modelu odpovidaji v IDENT prevracenym hodnotam
% korenu polynomu porade theta(z) a phi(z).
% Jejich odhady spolu se smerodatnymi odchylkami vraci
% ve specialnim formatu 'zepo' funkce 'th2zp'.
% 'zepo' = znaci ZEros & POles
% Pokud koreny phi(z) lezi vne jednotkoveho kruhu (tj. v IDENT uvnitr
% jednotkoveho kruhu), pak proces je kauzalni a stabilni, nebot
% 1/phi(z) je pak holomorfni v uzavrenem jednotkovem kruhu a da se zde
% rozvest do absolutne konvergentni Taylorovy rady.
% Podobne koreny theta(z) vne jednotkoveho kruhu (tj. v IDENT uvnitr
% jednotkoveho kruhu) zarucuji invertibilitu procesu.
% Rozmisteni korenu (o) a polu (x) a tedy invertibilitu, resp. kauzalnost
% a stabilitu systemu nejsnaze overime z grafu, ktery kresli funkce 'zpplot'.
% Pro AR proces se zobrazi pouze poly, nebot theta(z)=1 nema koreny.
% Provedte help th2zp a help zpplot.
keyboard;
%-----------------------------------------------------------------------------
% Tedy poly jsou ulozeny v
  zepo = th2zp(th)
% pricemz prvy sloupec submatice
  zp = zepo(2:p+1,3:4)
% jsou skutecne (az na poradi) prave prevracene hodnoty korenu
% polynomu phi1(z)
  phir  = roots(flipud([1;-PHI]))    % koreny pro presne parametry PHI
  phi1r = roots(flipud([1;-phi1]))   % koreny pro odhad phi1
  pol = 1./phi1r                     % poly dle IDENT
% a navic druhy sloupec zp pro tyto poly udava smerodatne odchylky
% a tedy polosirky 95% konfidencniho intervalu ziskame jako
  1.96*zp(:,2)
% Rozmisteni polu vykresli funkce 'zpplot':
  figure; zpplot(zepo)
% Vidime, ze poly skutecne lezi uvnitr jednotkoveho kruhu.
pause;
%-----------------------------------------------------------------------------
% Spocteme jeste odhady pomoci dalsich metod, ktere nabizi funkce 'ar':
% Metoda forward-backward: approach='fb'
  th1 = th;
  th = ar(x,p);   % implicitni metoda, tj. totez jako  th = ar(x,p,'fb');
  [par,P,LAM]=th2par(th);
  phi2 = -par.'; phi2E = 1.96*sqrt(diag(P)); v2 = LAM;
% Metoda least-squares: approach='ls'
  th = ar(x,p,'ls');
  [par,P,LAM]=th2par(th);
  phi3 = -par.'; phi3E = 1.96*sqrt(diag(P)); v3 = LAM;
% Metoda Burg's lattice: approach='burg'
  th = ar(x,p,'burg');
  [par,P,LAM]=th2par(th);
  phi4 = -par.'; phi4E = 1.96*sqrt(diag(P)); v4 = LAM;
% Metoda geometric lattice: approach='gl'
  th = ar(x,p,'gl');
  [par,P,LAM]=th2par(th);
  phi5 = -par.'; phi5E = 1.96*sqrt(diag(P)); v5 = LAM;
% Porovname vsechny odhady s presnymi hodnotami (1.sloupec):
% Parametry:
  [PHI,phi,phi1,phi2,phi3,phi4,phi5]
% Polosirky 95%-ho intervalu spolehlivosti:
  [zeros(p,1),phiE,phi1E,phi2E,phi3E,phi4E,phi5E]
% Rozptyl bileho sumu:
  [sigma2,v,v1,v2,v3,v4,v5]
% Vidime dobrou shodu vsech metod
pause;
%-----------------------------------------------------------------------------
close all
  th = th1; clear th1;
% 4. OVERENI MODELU ODHADNUTEHO JAKO AR(2)
%    -------------------------------------
% a) Nejprve zkontrolujeme, zda hodnoty koncovych parametru jsou
%    vyznamne (jejich konfidencni intervaly neobsahuji nulu):
  [phi1-phi1E,phi1+phi1E]
% Pokud by napriklad [phi1(p)-phi1E(p),phi1(p)+phi1E(p)] obsahoval nulu,
% znamenalo by to, ze phi1(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.
  [xs,z,vz] = armasim(x,phi1);  % uzijeme-li usek, z nehoz byl provaden
                                % odhad modelu
% Nebo jinak uzitim standardni funkce MATLABu 'filter':
  zz = filter([1;-phi1],1,x); xxs = x-zz; vzz = std(z).^2;
  max(abs(z-zz)), max(abs(xs-xxs)), vz-vzz
% Provedte help filter.
keyboard;
%-----------------------------------------------------------------------------
  [xvs,zv,vzv] = armasim(xv,phi1,[],x(n+1-p:n));
                                % uzijeme-li nezavisly usek
  zzv = filter([1;-phi1],1,[x(n+1-p:n);xv]); zzv = zzv(p+1:p+n);
  xxvs = xv-zzv; vzzv = std(zzv).^2;
  max(abs(zv-zzv)), max(abs(xvs-xxvs)), vzv-vzzv
% Proverime, vlastnost IID rezidulaniho sumu napr. pomoci 'trand'
% (existuji i specialni testy pouze na nekorelovanost, tj. jen na bily sum)
  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 q.
% Protoze AR(p) = ARMA(p,0), museli bychom v nasem pripade uvazovat
% o smisenem modelu pro vhodne q>0.
pause;
%-----------------------------------------------------------------------------
% Pro vypocet rezidualniho sumu z nabizi IDENT funkci 'pe', resp. 'resid',
% ktera navic zobrazi graf autokorelacni funkce pro z, jejiz prubeh by
% mel indikovat bily sum.
% Provedte help pe a help resid.
keyboard;
%-----------------------------------------------------------------------------
  z1 = pe(x,th);    % vypocet residualniho sumu
  acf(z1,figure);   % graf autokorelacni funkce pomoci 'acf'
pause;
%-----------------------------------------------------------------------------
% podobne, ale trochu jinak pro zv:
  ev1 = pe(xv,th);
  figure; zv1 = resid(xv,th);
  max(abs(ev1-zv1))
% Tedy skutecne 'pe' a 'resid' pocitaji totez.
pause;
  zvro1 = acf(zv1);
  figure; plot(zvro1(1:25));
% Porovnejte posledni dva obrazky - prubehy jsou stejne, i kdyz spocteny
% ruznymi funkcemi.
pause;
%-----------------------------------------------------------------------------
% Vidime, ze rovnez grafy autokorelacnich funkci potvrzuji bily sum pro
% z1 a zv1.
% Krome pocatecniho useku delky p (dano jinymi pocatecnimi podminkami pro x)
% se residualni sumy spoctene pomoci 'armasim' a 'resid' shoduji
%
  figure; zoom on; plot([z,z1]);    % Pouzijte zoom (viz help zoom)
pause;
  figure; zoom on; plot([zv,zv1]);    
pause;
  max(abs(z (p+2:n)-z1 (p+2:n)))
  max(abs(zv(p+2:n)-zv1(p+2:n)))
% Pro kontrolu jeste porovnejme vsechny odhady sigma2
  [v,v1,vz,vzv]
% s presnou hodnotou
  sigma2
pause;
%-----------------------------------------------------------------------------
close all;
% 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 IDENT predikce pocita funkce 'predict'.
%    Funkce 'compare' navic zobrazi do grafu ruznymi carami predikci
%    soucasne s pozorovanymi hodnotami.
%    Obe funkce predepsanim tretiho parametru umoznuji i predikci
%    z p hodnot o predepsany pocet kroku k>=1 zpet.
%    Provedte help predict a help compare.
keyboard;
%-----------------------------------------------------------------------------
% Spocteme nejprve predikce a porovname je s xs a xvs:
  xs1  = predict(x, th);
  xvs1 = predict(xv,th);
  max(abs(xs-xs1))
  max(abs(xvs(p+1:n)-xvs1(p+1:n)))
% Prvych p hodnot neporovnavame, nebot byly v 'armasim' generovany pomoci
% posledniho useku z x, coz 'predict' neumoznuje.
% Shoda je jinak uplna.
pause;
%-----------------------------------------------------------------------------
% Podobne postupujeme s 'compare', ktery navic zobrazi srovnavaci graf:
  figure; zoom on;
  xs2  = compare(x, th,1);
% Detaily lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
pause;
  figure; zoom on;
  xvs2 = compare(xv,th,1);
  max(abs(xs1-xs2))
  max(abs(xvs1-xvs2))
% 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;
%-----------------------------------------------------------------------------
% Kvalita modelu se nejvice pozna pri predikci na vetsi vzdalenost:
  figure; zoom on;
  compare(x, th,2);
pause;
  figure; zoom on;
  compare(x, th,3);
% I v techto pripadech je shoda tvaru do znacne miry zachovana, takze
% celkem muzeme nalezeny model AR(2) definitivne prijmout za verohodny.
pause;
%-----------------------------------------------------------------------------
%==================================KONEC======================================
echo off
