%==============================================================
% ARMA MODELOVANI 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-99
% Navazuje na zpracovani davkou airpass.m a vyzaduje
% NAG FOUNDATION TOOLBOX a vyuziva postupu demo davky 'crnag14'.
%========================================================================

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

format compact;
echo on
%==============================================================
% ARMA MODELOVANI 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-99
% Navazuje na zpracovani davkou airpass.m a vyzaduje
% NAG FOUNDATION TOOLBOX a vyuziva postupu demo davky 'crnag14'.
%========================================================================
%
% Zjisteni pritomnosti NAG FOUNDATION TOOLBOXu:
if ~exist('nagg')==3
   disp('CHYBA: NAG FOUNDATION TOOLBOX NENALEZEN');
   return;
end
pause;
%============================================================
% 1. NACTENI DAT PO PROVEDENI DEKOMPOZICE
%    ------------------------------------
  load airpass
  who
  plotdata(x,xinf);    title('Pozorovana data');
  plotdata(elog,xinf); title('Zlogaritmovana chybova komponenta');
  e = elog;
  n = length(e)
% Vzhledem k malemu poctu dat, neni vhodne pouzit cast dat pro
% overovani modelu.
pause;
%-----------------------------------------------------------------------------
% 2. IDENTIFIKACE PROCESU
%    --------------------
% Proverime nejprve nahodnost:
  ernd = trand(e-mean(e))
% Vsechny testy indikuji bily sum, coz svedci o dobre provedene dekompozici
% casove rady.
% Proverime presto moznost malezeni vhodneho ARMA modelu.
%
% a) ODHAD STREDNI HODNOTY A AUTOKOVARIANCNI FUNKCE PROCESU (procedura acf)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [ero,eroE,mi,miE,eacf]=acf(e,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 neindikuje bily sum, nebot zasahuje mimo konfidencni
%   pas nulovosti, pricemz nema nahodny charakter. Lze pozorovat dokonce
%   pravidelne se opakujici prubeh, coz svedci o zbytkovych korelacich
%   pochazejicich z periodickych komponent.
%   Proto usuzujeme, ze se pravdepodobne nebude jednat o MA proces.
pause;
%-----------------------------------------------------------------------------
% A nyni alternativne s vyuzitim procedury 'NAGacf':
  nq = sqrt(n)
  [NAGero,NAGev,NAGemi,state] = NAGacf(e,nq);
% Statistika 'stat' ma priblizne rozdeleni 'chi2(nq)' a umoznuje testovat
% hypotezu o nulovosti NAGero (tj. e je bily sum) s rizikem zamitnuti alpha:
  alpha = 0.05;
  kvantil=chi2inv(1-alpha,nq), state
% hypotezu o nulovosti NAGero zamitame s rizikem alpha:
  state>kvantil
pause;
%-----------------------------------------------------------------------------
% b) ODHAD PARCIALNI AUTOKORELACNI FUNKCE A AUTOREGRESNICH KOEFICIENTU phi
%    (procedura pacf)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [epacf,phi]=pacf(eacf,n,figure);
% Obalka parcialni autokorelacni funkce klesa exponencialne a velmi pomalu
% k nule, pricemz v konfidencnim pasu nulovosti nema nahodny charakter.
% Proto predbezne usuzujeme, ze se asi nejedna ani o AR proces konecneho radu.
pause;
%-----------------------------------------------------------------------------
close all
% c) POKUS O ODHAD RADU p PRO MODEL AR(p)
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Odhad v(p) z useku e pro kazdy rad p=1,...,pmax a
% nalezeni optimalnich radu pro nektera kriteria spocteme pomoci 'NAGpqest'.
% v(p)  ... odhad rozptylu spocteny z tehoz useku x uziteho pro odhad modelu
% v(p)=NaN, pokud nemohly byt ziskany spolehlive odhady.
% To zpravidla indikuje preparametrizovani (prilis velky rad p).
  pmax = 10;
  [pqo,v,vv,ke,pq,rp] = NAGpqest(pmax,0,e,[],figure);
  pqo(1,:)   % optimalni rady p spoctene porade z FPE,AIC a MDL
  pq(:,1:ke) % rady p,q, z nichz mohl byt uspesne proveden optimalni vyber
             % v(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;
%-----------------------------------------------------------------------------
% Kriterialni funkce FPE, AIC i MDL klesaji.
% Pouze silne modifikujici kriterium v10 nabizi moznost AR(4).
% Vzhledem k prubehum 'acf' a 'pacf' ve 2a)b) se ale priklanime spise
% k neadkvatnosti AR modelu a pokusime se najit smiseny ARMA model.
pause;
%-----------------------------------------------------------------------------
close all
% d) ODHAD RADU p,q PROCESU ARMA(p,q), p>=0, q>0
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Odhad v(q) z useku x 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 useku e uziteho pro odhad modelu
% v(k)=NaN, pokud nemohly byt ziskany spolehlive odhady.
% To zpravidla indikuje preparametrizovani (prilis velke p+q).
  pmax=5; qmax = 5;
  [pqo,v,vv,ke,pq,rp] = NAGpqest(pmax,qmax,e,[],figure);
  pqo        % optimalni rady p,q spoctene porade z FPE,AIC a MDL
  pq(:,1:ke) % rady p,q, z nichz mohl byt uspesne proveden optimalni vyber
             % v(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(1:ke))
  pq(:,kmin) % optimalni rady p,q dle kriteria v10
%
% Vsechny kriterialnich funkce indikuji rady p=2,q=1.
pause;
%-----------------------------------------------------------------------------
close all
% 3. ODHAD PARAMETRU PROCESU IDENTIFIKOVANEHO JAKO ARMA(2,1)
%    -------------------------------------------------------
% 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.
%
% a) Predbezny odhad pomoci 'NAGpest':
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  p = 2; q = 1;
  mr = ARMA(p,q)    % vektor radu pro SARIMA model redukovany na ARMA(p,q).
  [par0,v0,isf] = NAGpest(mr,NAGero,NAGev);
% 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
% Predbezny odhad parametru phi nebyl uspesny.
pause;
%-----------------------------------------------------------------------------
% b) Finalni odhad pomoci 'NAGest':
%    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  [par1,par1S,v1,est1INFO] = NAGest(mr,e,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}
  est1INFO.xdmi      % odhad stredni hodnoty
  v1                 % rozptylu sigma2
  est1INFO.isf(1:2)  % indikatory uspesnosti
  est1INFO.itc       % pocet iteraci
% Finalni odhad parametru theta pomoci metody nejmensich ctvercu nebyl
% uspesny.
pause;
%-----------------------------------------------------------------------------
% Zkusime to znovu, tentokrat bez odhadovani stredni hodnoty, kterou
% zvolime pevne (emiEST = 0) jako prumer:
  emi0 = mean(e);
  emiEST = 0;
  [par2,par2S,v2,est2INFO] = NAGest(mr,e,par0,[],emi0,emiEST);
% Finalni odhady
  phi2    = par2{1}  % parametru phi
  phi2S   = par2S{1} % jejich smerodatnych odchylek
  theta2  = par2{2}  % parametru theta
  theta2S = par2S{2} % jejich smerodatnych odchylek
  emi2  = est2INFO.xdmi       % stredni hodnota (pevna)
  v2                 % rozptylu sigma2
  est2INFO.isf(1:2)  % indikatory uspesnosti
  est2INFO.itc       % pocet iteraci
% Finalni odhad parametru theta opet nebyl zcela uspesny
pause;
%-----------------------------------------------------------------------------
% Definitivni odhad spocteme metodou maximalni verohodnosti z NAG procedury
% 'G13BEF':
  [par3,par3S,v3,est3INFO] = NAGest(mr,e,par0,'MLE');
% Finalni odhady
  phi3    = par3{1}       % parametru phi
  phi3E   = 1.96*par3S{1} % jejich konfidencnich polosirek
  theta3  = par3{2}       % parametru theta
  theta3E = 1.96*par3S{2} % jejich konfidencnich polosirek
  emi3  = est3INFO.xdmi       % odhad stredni hodnoty
  emi3E = 1.96*est3INFO.xdmiS % a jeji konfidencni polosirky
  v3                      % rozptylu sigma2
  est3INFO.isf(1:2)       % indikatory uspesnosti
  est3INFO.itc            % pocet iteraci
% Tentokrat jsou jiz odhady uspesne.
% Pocet iteraci dosahl implicitniho limitu 50. Proto jej zvetsime napr.
% na 1000:
pause;
%-----------------------------------------------------------------------------
  [par3,par3S,v3,est3INFO] = NAGest(mr,e,par0,'MLE',[],[],[],1000);
  est3INFO.itc
% Ani 1000 iteraci nestaci. Vznika podezreni, ze model neni dobry.
% Zkusime jej proverit:
pause;
%-----------------------------------------------------------------------------
% Proverime vlastnost IID rezidulaniho sumu pomoci 'trand':
  z = est3INFO.z;
  trand(z)
% a rezidua pomoci 'NAGcheck':
  [zro3,zro3cm,stat3,idf3,siglev3] = NAGcheck(mr,par3,z,sqrt(n));
  idf3, siglev3
  0.05>siglev3
% Relace vyrazne plati, coz vede k zamitnuti hypotezy o bilem sumu.
% Dokonce nektere testy z 'trand' popiraji nahodnost 'z', coz je horsi
% vysledek, nez u puvodnich dat 'e'.
% Model zrejme nevystihuje dobre data a pocet parametru je pravdepodobne
% nedostatecny. Zvolime vyssi rad q, jelikoz je jen 1 a byly problemy
% s odhadem theta;
pause;
close all;
%-----------------------------------------------------------------------------
  p = 2; q = 2;
  mr = ARMA(p,q);
% Pro pocatecni odhad pouzijeme jako vychozi odhad phi3 a theta3,
% pricemz v par0 polozime theta0(2)=0:
  par0 = par3; par0{2} = [par3{2},0];
  [par4,par4S,v4,est4INFO] = NAGest(mr,e,par0,'MLE');
  est4INFO.itc
% Zvetsime jeste pocet iteraci
pause;
%-----------------------------------------------------------------------------
  [par4,par4S,v4,est4INFO] = NAGest(mr,e,par0,'MLE',[],[],[],500);
  est4INFO.itc
% Dosazeny definitivni odhady po 112 iteracich, coz je vyhovujici.
% Provedme opet ihned provereni:
  z = est4INFO.z;
  trand(z)
% a rezidua pomoci 'NAGcheck':
  [zro4,zro4cm,stat4,idf4,siglev4] = NAGcheck(mr,par4,z,sqrt(n));
  siglev3, siglev4
  0.05>siglev4
% Opet je zamitan bily sum, i kdyz 'siglev4' je vetsi nez 'siglev3',
% coz je blize k nezamitnuti.
% Vyzkousime dalsi modely:
pause;
%-----------------------------------------------------------------------------
  p = 3; q = 1;
  mr = ARMA(p,q);
  par0 = par3; par0{1} = [par3{1},0];
  [parX,parXS,vX,estXINFO] = NAGest(mr,e,par0,'MLE',[],[],[],500);
  estXINFO.itc
% Provedme opet ihned provereni:
  z = estXINFO.z;
  trand(z)
% a rezidua pomoci 'NAGcheck':
  [zroX,zroXcm,statX,idfX,siglevX] = NAGcheck(mr,parX,z,sqrt(n));
  siglev4, siglevX
  0.05>siglevX
% Situace se zhorsila (siglevX se vyrazne zmensila oproti siglev4), 
% zkusime jiny model
pause;
%-----------------------------------------------------------------------------
  p = 2; q = 3;
  mr = ARMA(p,q);
  par0 = par3; par0{2} = [par3{2},0,0];
  [parX,parXS,vX,estXINFO] = NAGest(mr,e,par0,'MLE',[],[],[],500);
  estXINFO.itc
% Provedme opet ihned provereni:
  z = estXINFO.z;
  trand(z)
% a rezidua pomoci 'NAGcheck':
  [zroX,zroXcm,statX,idfX,siglevX] = NAGcheck(mr,parX,z,sqrt(n));
  siglev4, siglevX
  0.05>siglevX
% Je to lepsi: 'siglevX' vzrostlo oproti 'siglev4'
  siglevm = siglevX;
% Zkusime jeste zvetsit q:
pause;
%-----------------------------------------------------------------------------
  p = 2; q = 4;
  mr = ARMA(p,q);
  par0 = parX; par0{2} = [parX{2},0];
  [parX,parXS,vX,estXINFO] = NAGest(mr,e,par0,'MLE',[],[],[],500);
  estXINFO.itc
% Provedme opet ihned provereni:
  z = estXINFO.z;
  trand(z)
% a rezidua pomoci 'NAGcheck':
  [zroX,zroXcm,statX,idfX,siglevX] = NAGcheck(mr,parX,z,sqrt(n));
  siglevm, siglevX
  0.05>siglevX
% Situace se ponekud zhorsila, dalsi zvetsovani q nebude ta spravna cesta.
% Zkusime jeste AR(3,0):
pause;
%-----------------------------------------------------------------------------
  p = 3; q = 0;
  mr = ARMA(p,q);
  par0 = par3; par0{1} = [par3{1},0]; par0{2}=[];
  [parX,parXS,vX,estXINFO] = NAGest(mr,e,par0,'MLE',[],[],[],500);
  estXINFO.itc
% Provedme opet ihned provereni:
  z = estXINFO.z;
  trand(z)
% a rezidua pomoci 'NAGcheck':
  [zroX,zroXcm,statX,idfX,siglevX] = NAGcheck(mr,parX,z,sqrt(n));
  siglevm, siglevX
  0.05>siglevX
% Vyrazne zhorseni.
% Dalsi experimenty s rady modelu jiz zrejme nepovedou k vyraznemu zlepseni.
% Vratime se k ARMA(2,3) a podivame se na vyznamnost parametru theta(3):
pause;
%-----------------------------------------------------------------------------
  p = 2; q = 3;
  mr = ARMA(p,q);
  par0 = par3; par0{2} = [par3{2},0,0];
  [par5,par5S,v5,est5INFO] = NAGest(mr,e,par0,'MLE',[],[],[],500);
  theta5 = par5{2};
  theta5E = 1.96*par5S{2};
  [theta5(q)-theta5E(q),theta5(q)+theta5E(q)]
% Interval spolehlivosti zahrnuje nulu => theta5(q) je nevyznamny.
% Bude tedy rozumne se vratit k odhadum c.4 pro ARMA(2,2).
pause;
%-----------------------------------------------------------------------------
% Finalni odhady
  p = 2; q = 2;
  mr = ARMA(p,q);
  phi4    = par4{1}       % parametru phi
  phi4E   = 1.96*par4S{1} % jejich konfidencnich polosirek
  theta4  = par4{2}       % parametru theta
  theta4E = 1.96*par4S{2} % jejich konfidencnich polosirek
  emi4  = est4INFO.xdmi       % odhad stredni hodnoty
  emi4E = 1.96*est4INFO.xdmiS % a jeji konfidencni polosirky
  v4                      % rozptylu sigma2
pause;
%-----------------------------------------------------------------------------
% Prohledneme jeste obsah struktury estINFO, ktera obsahuje dodatecne
% informace:
  est4INFO
% kde je:
         % xdmi        ... odhad stredni hodnoty
         % xdmiS       ... jeji smerodatna odchylka
  est4INFO.corrmx   %  ... korelacni matice odhadu vsech parametru
                    %      (celkem pet parametru vcetne 'xdmi')
         % isf         ... indikatory uspesnosti
         % stset       ... stavova mnozina pro dalsi zpracovani
  z = est4INFO.z;
         % 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;
%-----------------------------------------------------------------------------
% 4. ZAVERECNE OVERENI MODELU ODHADNUTEHO JAKO ARMA(2,2)
%    ---------------------------------------------------
% a) Nejprve zkontrolujeme, zda hodnoty (zejmena koncovych) parametru jsou
%    vyznamne (jejich konfidencni intervaly neobsahuji nulu):
  phi4 = phi4'; phi4E = phi4E';
  [phi4-phi4E,phi4+phi4E]
  theta4 = theta4'; theta4E = theta4E';
  [theta4-theta4E,theta4+theta4E]
% Neobsahuji nulu, takze parametry phi4(p) a theta4(q) lze povazovat za
% vyznamne.
pause;
%-----------------------------------------------------------------------------
% b) Proverime, zda odhadnuty model a jeho komponenty splnuji teoreticke
%    predpoklady, za nichz byl nalezen.
  [ans,phi4r,theta4r] = armagen(phi4,theta4,1);
  format long
  abs(phi4r)   % koreny phi(z) jsou vne jednotkoveho kruhu => kauzalita
  abs(theta4r)
  format short
%    jeden koren theta(z) je jen tesne vne jednotkoveho kruhu =>
%    invertibilita je jen tesne dosazena;
%    blizkost korene k jednotkovemu kruhu byla pravdepodobne pricinou
%    predchozich obtizi pri odhadu parametru theta.
pause;
%-----------------------------------------------------------------------------
% Proverime vlastnost IID rezidulaniho sumu pomoci 'trand':
  trand(z)
% OK.
pause;
%-----------------------------------------------------------------------------
% Spocteme jeste autokorelacni funkci rezidualniho sumu z pomoci
% 'acf' a 'NAGacf' a proverime, zda se jedna o bily sum:
  nqz = round(sqrt(n));
  [zro4,zro4cm,stat4,idf4,siglev4] = NAGcheck(mr,par4,z,nqz);
  idf4, siglev4
  0.05>siglev4
  acf(z ,figure); xlabel('z rezidui');
% porovname s
  acf(e ,figure); xlabel('z puvodnich dat');
% Rezidua 'z' maji evidentne blize k bilemu sumu, nez 'e', takze
% model ARMA(2,2) byl prinosny, i kdyz bily sum byl zamitnut
% testem z 'NAGcheck' s rizikem 0.05.
pause;
%-----------------------------------------------------------------------------
% Podobne uzitim 'NAGacf':
  [zro,ans,ans,stat] = NAGacf(z,nqz);
  stat, kvantil=chi2inv(0.95,nqz-(p+q))
  stat>kvantil
% Tato relace je splnena, takze tento test rovnez zamita hypotezu o
% bilem sumu s rizikem alpha=0.05.
  state,stat
% Hodnota statistiky je vsak vyrazne mensi ve srovnani se 'state' spoctene
% z puvodnich dat, coz rovnez svedci ve prospech modelu ARMA(2,2).
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.
pause;
%-----------------------------------------------------------------------------
% Porovname 1-krokovou predikci 'es' s puvodnimi daty e:
   es = e - z;
   plot(e,'y'); hold on; plot(es,'m'); hold off;
   title('Pozorovana data e (zluta), predikce (fialova)');
   zoom on
% Detaily lepe uvidite s vyuzitim zoomovani pomoci mysi (viz help zoom).
% Vidime velmi dobrou shodu.
pause;
%-----------------------------------------------------------------------------
  zoom off
% Spocteme nyni orientacne predikci o 10 kroku pomoci 'NAGpred':
  npv = 10;   % delka predikce
  [ep4,ep4S] = NAGpred(mr,par4,npv,est4INFO);
% polosirka intervalu spolehlivosti pri normalni aproximaci chyb:
  ep4E = 1.96*ep4S;
  figure; plot(ep4,'m'); hold on; plot([ep4-ep4E,ep4+ep4E],'m--'); hold off;
  title('Predikce');
pause;
%-----------------------------------------------------------------------------
% Na zaver provedeme korekci trendu prictenim predikovane (vyhlazene)
% casti es (totiz elog = e = es + z):
% xlog = trlog + szlog + (es +  z) = (trlog + es) + szlog +  z:
  trlogs = trlog + es;
  figure; plot([trlog,trlogs]);
pause;
% Po odlogaritmovani:
  trs = exp(trlogs);
  figure; plot([tr,trs]);
pause;
%-----------------------------------------------------------------------------
% Tedy vysledny trend (bez sezonni slozky) a 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,:));
% Proverime nahodnost po odlogaritmovani:
  ez = exp(z);
  trand(ez)
% Vsechny testy ji potvrzuji
%==================================KONEC======================================
echo off
