%==============================================================
% ANALYZA NAHODNE SLOZKY SOUBORU airpass.dat:
% Cestujici v mezinarodni letecke doprave
% x=mesicni souhrny v tisicich cestujicich
% t=1,2,...,144: leden 1949 az prosinec 1960
%
% Zpracoval: V. Vesely, KAM PrF MU (C) 1996-98
% Navazuje na zpracovani davkou airpass.m a vyzaduje
% SYSTEM IDENTIFICATION TOOLBOX a vyuziva postupu demo davky 'crcv14'.
%========================================================================
format compact
echo on
%-----------------------------------------------------------------------------
% Zjisteni pritomnosti SYSTEM IDENTIFICATION TOOLBOXu:
if ~exist('ident')==2
   disp('CHYBA: SYSTEM IDENTIFICATION TOOLBOX NENALEZEN');
   return;
end
%============================================================
% 1. NACTENI DAT PO PROVEDENI DEKOMPOZICE
%    ------------------------------------
  load airpass
  who
  plotdata(x,xinf);    title('Pozorovana data');
  plotdata(elog,xinf); title('Zlogaritmovana chybova komponenta');
  n = length(x)
% Vzhledem k malemu poctu dat, neni vhodne pouzit cast dat pro
% overovani modelu.
pause;
%-----------------------------------------------------------------------------
% 2. IDENTIFIKACE PROCESU
%    --------------------
% a) ODHAD STREDNI HODNOTY A AUTOKOVARIANCNI FUNKCE PROCESU (procedura acf)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [elogro,elogroE,mi,miE,elogacf]=acf(elog,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)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [elogpacf,phi]=pacf(elogacf,n,figure);
% Obalka parcialni autokorelacni funkce klesa exponencialne k nule, pricemz
% v konfidencnim pasu nulovosti nema nahodny charakter.
% Proto usuzujeme, ze se nejedna o AR proces.
% Bude se tedy asi jednat o smiseny proces ARMA.
% Proverime jeste nahodnost
  elogrnd = trand(elog)
% Testy nahodnosti potvrzuji bily sum, coz je ponekud v rozporu
% s prubehem acf a pacf.
% Pokusime se proto o nalezeni vhodneho ARMA modelu.
pause;
%-----------------------------------------------------------------------------
close all
% c) POKUS O ODHAD RADU p PRO MODEL AR(p)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Spocteme odhady rozptylu generujiciho bileho sumu pro sirsi rozsah
% radu, napr. pro p=1:20:
  pmax = 20;
  V = arxstruc(elog,elog,(1:pmax).'); % Ztratova funkce jen z x
  figure; p = selstruc(V)       % graf ztratove funkce, v p nabyva minima
pause;
%-----------------------------------------------------------------------------
% Vykreslime grafy nekolika kriterialnich ztratovych funkci soucasne
  rp = (1:pmax).'./n; AIC = 2; MDL = log(n);
  v  = V(1,1:pmax).';
  vaic = v.*(1+ AIC*rp);
  vmdl = v.*(1+ MDL*rp);
  plot(v,'w-');
  hold on; plot(vaic,'g--'); plot(vmdl,'g:');
  title('v(p): Ztratova funkce (-), AIC (--) a MDL (:)');
  xlabel('rad p');
  hold off;
  figure(gcf);
pause;
%-----------------------------------------------------------------------------
% Prohledneme kriterialni funkce jeste pro c=3:8:
  c  = 3:8;
  vc = (v*ones(size(c))).*(1+rp*c);
  figure; plot(vc);
  title('Kriterialni ztratove funkce pro c=3,...,8');
  xlabel('rad p');
% a spocteme minima
  p = selstruc(V,'aic')  % p, kde AIC nabyva minima
  p = selstruc(V,'mdl')  % p, kde MDL nabyva minima
  for cc=c
      cc
      p = selstruc(V,cc) % p, kde kriterium s vahou cc nabyva minima
pause;
  end
%-----------------------------------------------------------------------------
% Body zlomu ztratovych funkci, resp. minima kriterialnich funkci jsou
% nevyrazna, davaji minimum pro 16.
% Tak vysoky rad je nepravdepodobny a vzhledem k malemu poctu dat n=144
% nevhodny.
% To potvrzuje nase ocekavani o neadkvatnosti AR modelu, ktere naznacil
% jiz prubeh acf.
% V dalsim se pokusime najit optimalni ARMA model.
pause;
%-----------------------------------------------------------------------------
close all
% d) ODHAD RADU p,q PROCESU ARMA(p,q), q>0
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Pro ARMA model IDENT bohuzel nenabizi (ponekud nelogicky) funkce,
% ktere by byly analogii funkci 'arxstruc' a 'selstruc' pro AR model.
% Pro odhad parametru modelu ARMA vsak muzeme pouzit funkci 'armax',
% kterou budeme opakovane volat pro vsechny kombinace p,q z vhodne
% zvoleneho rozsahu radu, napr. pro p=0:5, q:1:5.
% Pomoci 'th2par', resp. 'pe' pak spocteme odhady rozptylu sigma2
% generujiciho bileho sumu:
% v(p)  ... odhad rozptylu spocteny z tehoz useku x uziteho pro odhad modelu
  pmax = 8; qmax = 8;
% Vypocet dvojic [p,q] v uvedenem rozsahu:
  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); [pq;Npq]
  pq = pq.';
  kmax = size(pq,1)    % pocet dvojic
% Provedte help 'armax'
keyboard;
%-----------------------------------------------------------------------------
% Odhad v(k) z useku x pro kazdou dvojici pq(k,:):
  for k=1:kmax
      th = armax(elog,pq(k,:));
      [par,P,v(k)] = th2par(th);
%      z = pe(x ,th); v(k)  = std(z).^2;
echo off
  end
echo on
% Vykreslime prubeh v(k) a prislusnych kriterialnich funkci AIC a MDL:
  rp   = Npq.' ./n;
  vaic = v.*(1+ AIC*rp);
  vmdl = v.*(1+ MDL*rp);
  vc3  = v.*(1+   3*rp);
  plot(v,'w-');
  hold on; plot(vaic,'g--'); plot(vmdl,'g:');
  title('v(k): Ztratova funkce (-), AIC (--) a MDL (:)');
  xlabel('k');
  hold off;
pause;
% a spocteme minima
  [vm,I]    = min(v)   ; pq(I,:), vm     % p,q kde v   nabyva minima
  [vaicm,I] = min(vaic); pq(I,:), vaicm  % p,q kde AIC nabyva minima
  [vmdlm,I] = min(vmdl); pq(I,:), vmdlm  % p,q kde MDL nabyva minima
% Jeste pro c=3:8
  c  = 3:8;
  vc = (v*ones(size(c))).*(1+rp*c);
  figure; plot(vc);
  title('Kriterialni ztratove funkce pro c=3,...,8');
  xlabel('rad p');
% a spocteme minima
  for k=1:length(c)
      c(k)
      [vcmin,I] = min(vc(:,k)); I, pq(I,:), vcmin  % p,q kde vc nabyva minima
pause;
  end
% Vzhledem k prubehu grafu se jako nejvhodnejsi jevi p = 2, q = 2.
pause;
%-----------------------------------------------------------------------------
close all
% 3. ODHAD PARAMETRU PROCESU IDENTIFIKOVANEHO JAKO ARMA(2,2)
%    -------------------------------------------------------
  p = 2, q = 2
  th = armax(elog,[p,q]);
  present(th);
  [par,P,LAM]=th2par(th)
% kde
  phi   = - par(1:p).'
  theta = par(p+1:p+q).'
% je odhad parametru
% spolu s polosirkami 95%-ho konfidencniho intervalu
  phithetaE = 1.96*sqrt(diag(P));
  phiE   = phithetaE(1:p)
  thetaE = phithetaE(p+1:p+q)
% a
  v = LAM
% je odhad rozptylu generujiciho bileho sumu.
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'.
keyboard;
%-----------------------------------------------------------------------------
% Tedy koreny i poly jsou ulozeny v
  zepo = th2zp(th)
% pricemz prvy sloupec submatice
  zp = zepo(2:q+1,1:4)
% jsou skutecne (az na poradi) prave prevracene hodnoty korenu
% polynomu theta(z)
  thetar = roots(flipud([1;theta])) % koreny pro odhad theta
  zer = 1./thetar                   % koreny dle IDENT
% a navic druhy sloupec zp pro tyto koreny udava smerodatne odchylky
% a tedy polosirky 95% konfidencniho intervalu ziskame jako
  1.96*zp(:,2)
pause;
%-----------------------------------------------------------------------------
% Podobne treti sloupec zp jsou skutecne (az na poradi) prave
% prevracene hodnoty korenu polynomu phi(z)
  phir = roots(flipud([1;-phi])) % koreny pro odhad phi
  pol = 1./phir                  % poly dle IDENT
% a ctvrty sloupec zp pro tyto koreny udava smerodatne odchylky
% a tedy polosirky 95% konfidencniho intervalu ziskame jako
  1.96*zp(:,4)
% Rozmisteni korenu a polu vykresli funkce 'zpplot':
  figure; zpplot(zepo)
% Vidime, ze koreny (o) i poly (x) skutecne lezi uvnitr jednotkoveho kruhu.
pause;
%-----------------------------------------------------------------------------
close all
% 4. OVERENI MODELU ODHADNUTEHO JAKO ARMA(2,2)
%    -----------------------------------------
% a) Nejprve zkontrolujeme, zda hodnoty koncovych parametru jsou
%    vyznamne (jejich konfidencni intervaly neobsahuji nulu):
  [phi-phiE,phi+phiE]
  [theta-thetaE,theta+thetaE]
% Ani jeden z intervalu spolehlivosti neobsahuje nulu a tedy vsechny
% parametry se vyznamne lisi od nuly.
% Pokud by napriklad [theta(q)-thetaE(q),theta(q)+thetaEqp)] obsahoval nulu,
% znamenalo by to, ze theta(q) se nelisi vyznamne od nuly a ze tedy
% rad modelu je pravdepodobne nadhodnocen.
% Cely proces odhadu bychom pak zopakovali pro ARMA(p,q-1).
% Podobne bychom presli na model ARMA(p-1,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.
  [elogs,z,vz] = armasim(elog,phi,theta);  % uzijeme-li usek, z nehoz byl provaden
                                           % odhad modelu
% Proverime, vlastnost IID rezidulaniho sumu napr. pomoci 'trand'
% (existuji i specialni testy pouze na nekorelovanost, tj. jen na bily sum)
  trand(z)
% Az na jeden test je nahodnost potvrzena.
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 = resid(elog,th);    % vypocet residualniho sumu + graf autokor. funkce
pause;
  z1ro=acf(z1,figure);    % graf autokorelacni funkce pomoci 'acf'
pause;
  pacf(z1ro,n,figure);    % graf parcialni autokorelacni funkce
pause;
%-----------------------------------------------------------------------------
% Vidime, ze rovnez grafy acf a pacf potvrzuji jiz bily sum pro z1.
% Residualni sumy spoctene pomoci 'armasim' a 'resid' se shoduji:
%
  figure; plot([z,z1])
  max(abs(z-z1))
pause;
%-----------------------------------------------------------------------------
close all;
% c) na zaver overime vizualne z grafu soulad pozorovanych hodnot
%    s hodnotami predikovanymi z q hodnot pozorovanych o 1 krok zpet
%    pri vyuziti odhadnuteho modelu.
%    Tuto predikci jiz drive spocetla funkce 'armasim' v 'xs'.
%    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:
  elogs1  = predict(elog, th);
  max(abs(elogs-elogs1))
  figure; zoom on; plot([elogs,elogs1]);   % Pouzijte zoom (viz help zoom)
% Shoda je takrka uplna.
pause;
%-----------------------------------------------------------------------------
% Podobne postupujeme s 'compare', ktery navic zobrazi srovnavaci graf:
  figure; zoom on;
  elogs2  = compare(elog, th,1);
% Detaily lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
  max(abs(elogs1-elogs2))
pause;
%-----------------------------------------------------------------------------
% Kvalita modelu se nejvice pozna pri predikci na vetsi vzdalenost:
  figure; zoom on;
  compare(elog, th,2);
pause;
  figure; zoom on;
  compare(elog, th,3);
% I v techto pripadech je shoda tvaru prekvapive dobra, takze
% celkem muzeme nalezeny model ARMA(2,2) definitivne prijmout za verohodny.
pause;
%-----------------------------------------------------------------------------
% Na zaver provedeme korekci trendu prictenim predikovane casti elogs
% (totiz elog = elogs + z):
% xlog = trlog + szlog + (elogs +  z) = (trlog + elogs) + szlog +  z:
  trlogs = trlog + elogs;
  figure; plot([trlog,trlogs]);
pause;
  trs = exp(trlogs);
  figure; plot([tr,trs]);
pause;
%-----------------------------------------------------------------------------
% Tedy vysledny trend spolu s konfidencnimi pasy je po odlogaritmovani tvaru:
  t = 1:n;
  figure;
  plot(t,x,'go',t,trs,'w-',...
       t,[trs.*trE(:,1),trs./trE(:,1)],'y--',...
       t,[trs.*trE(:,2),trs./trE(:,2)],'r--');
  title('Data x (o), tr (-)');
  xlabel(xinf(2,:));
  ez = exp(z);
  trand(ez)
% Az na jeden test potvrzeno opet nahodne chovani i po odlogaritmovani.
%==================================KONEC======================================
echo off
