%==============================================================
% ANALYZA 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-97
%==============================================================
echo on
format compact
  air_pass=[crpath,'/brockwel.dat/airpass.dat'];
  [x,xinf]=getdata(air_pass);   % nacteni a vykresleni
  n=length(x)
  t=(1:n).';
% Graf ukazuje na pritomnost trendove i peridockych slozek.
% Amplituda kmitu roste s trendem, takze bude pravdepodobne porusena
% konstantnost rozptylu chybove slozky (viz dale 1)
pause;
%---------------------------------------------------------------
% Pro blizsi posouzeni spocteme autokorelacni funkci a jeji graf:
  acf(x,figure);
  xlabel(xinf(2,:));
% Pomaly pokles a oscilacni charakter autokorelacni funkce
% potvrzuje rovnez pritomnost trendove i peridockych slozek
pause;
%---------------------------------------------------------------
% 1. Odhad parametru lambda pro mocninnou transformaci
%---------------------------------------------------------------
  [lambdaE,y] = powtr(x,4,figure);
  xlabel(xinf(2,:));
  title('Odhad parametru lambda pro mocninnou transformaci');
% Byl tedy potvrzen linearni rust rozptylu, nebot
% 0 lezi v konfid. intervalu
  [lambdaE(1)-lambdaE(2),lambdaE(1)+lambdaE(2)]
% To znamena, ze puvodni data x se ridi multiplikativnim modelem.
  pause;
%---------------------------------------------------------------------------
% Logaritmicky transformovana data xlog=y se pak ridi aditivnim modelem:
  xlog = y;
  figure; plot(t,xlog);
  xlabel(xinf(2,:));
  title('Logaritmicka transformace');
pause;
%---------------------------------------------------------------------------
% 2. Detekce periodicke slozky
%---------------------------------------------------------------------------
% 2.1. Detekce pomoci Fisherova a Siegelova testu periodicity:
  [pdg,kF,kS] = tperiod(xlog,figure); % periodogram pdg a jeho graf
  title('PERIODOGRAM');
  xlabel(xinf(2,:));
  kF        % indexy vyznamnych harmonickych slozek uzitim Fisherova testu
  kS        % indexy vyznamnych harmonickych slozek uzitim Siegelova testu
% o periodach:
  if length(kS)>length(kF)
     n./kS
  else
     n./kF
  end
% Mezi vyznamne slozky patri kF(4)=kS(4)=12 (12-ta harmonicka) indikujici
% podle ocekavani periodu 144/12 = 12 mesicu (rocni perioda).
  d = n./kF(4)
pause;
%---------------------------------------------------------------------------
% 2. Separace periodicke slozky
%---------------------------------------------------------------------------
  close all;
% Vyzkousime vsechny dostupne metody a porovname vysledky:
% a) metoda maleho trendu:
%    ---------------------
  MM=1; [sz1,tr1] = szsmtr(x,d,MM);             % Multiplikativni model
  AM=0; [szlog1,trlog1] = szsmtr(xlog,d,AM);    % Aditivni model
% Vykreslime soucasne oba odhady
% - periodicke slozky:
  figure; plot(t,[sz1,exp(szlog1)]);  % szlog1 z AM je nutno odlogaritmovat
  title('Sezonni slozka metodou maleho trendu v MM a AM');
  xlabel(xinf(2,:));
pause;
% - trendu
  figure; plot(t,[tr1,exp(trlog1)]);  % trlog1 z AM je nutno odlogaritmovat
  title('Trendova slozka metodou maleho trendu v MM a AM');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% b) metoda klouzaveho prumeru:
%    -------------------------
  MM=1; [sz2,tr2] = szma(x,d,MM);               % Multiplikativni model
  AM=0; [szlog2,trlog2] = szma(xlog,d,AM);      % Aditivni model
% Vykreslime soucasne oba odhady
% - periodicke slozky:
  figure; plot(t,[sz2,exp(szlog2)]);   % szlog2 z AM je nutno odlogaritmovat
  title('Sezonni slozka metodou klouzaveho prumeru v MM a AM');
  xlabel(xinf(2,:));
pause;
% - trendu
  figure; plot(t,[tr2,exp(trlog2)]);  % trlog2 z AM je nutno odlogaritmovat
  title('Trendova slozka metodou klouzaveho prumeru v MM a AM');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% c) regresni metoda pomoci Fourierovy rady (procedura 'polfs')
%    ----------------------------------------------------------
% Nejdrive je treba odhadnout stupen polynomu pro trend.
  xlog_sz2 = xlog - szlog2;  % = Tr+e po odecteni sezonni slozky odhadnute v b)
  figure; plot(t,xlog_sz2);
% Zkusime nejdrive parabolicky trend:
pause;
%---------------------------------------------------------------------------
  format long
  close;
  [y,E,r,s,bE]=polyfilt(2,n,xlog_sz2,t,[],[],figure);
  [bE(:,1)-bE(:,2),bE(:,1)+bE(:,2)] % Interval spolehlivosti pro koef. polynomu
% Koeficient u nejvyssi mocniny (1. radek) je vyznamny, nebot jeho interval
% nezahrnuje nulu.
pause;
%---------------------------------------------------------------------------
% Vyzkousime kubicky trend:
  [y,E,r,s,bE]=polyfilt(3,n,xlog_sz2,t,[],[],figure);
  [bE(:,1)-bE(:,2),bE(:,1)+bE(:,2)] % Interval spolehlivosti pro koef. polynomu
% Tentokrat je nula v intervalu pro koeficient u nejvyssi mocniny (1. radek),
% koeficient lze tedy povazovat za nulovy vzhledem k male sirce intervalu.
pause;
%---------------------------------------------------------------------------
% Pro regresi uzitim 'polfs' budeme tedy predpokladat kvadraticky trend
  close;
  kh = kF; kh(1) = []    % 1. harmonickou vyloucime
  [yE3,trE3,pE3] = polfs(xlog,2,kh,t,[],[],figure);
  szlog3 = pE3(:,1);
  trlog3 = trE3(:,1);
pause;
%---------------------------------------------------------------------------
% d) regresni metoda pro primy odhad sezonni slozky (procedura 'polper')
%    ------------------------------------------------------------------------
  [yE4,trE4,pE4] = polper(xlog,2,d,t,[],[],figure);
  szlog4 = pE4(:,1);
  trlog4 = trE4(:,1);
%---------------------------------------------------------------------------
pause;
% Vykreslime soucasne oba odhady c) a d) po odlogaritmovani
% - periodicke slozky:
  figure; plot(t,exp([szlog3,szlog4]));
  title('Sezonni slozka regresni metodou v AM');
  xlabel(xinf(2,:));
pause;
% - trendu
  figure; plot(t,exp([trlog3,trlog4]));
  title('Trendova slozka regresni metodou v AM');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% Vykreslime nyni vsechny odhady a) az d) v AM po odlogaritmovani
% - periodicke slozky:
  figure; plot(t,exp([szlog1,szlog2,szlog3,szlog4]));
  title('Sezonni slozka vsemi metodami v AM');
  xlabel(xinf(2,:));
pause;
% - trendu
  figure; plot(t,exp([trlog1,trlog2,trlog3,trlog4]));
  title('Trendova slozka vsemi metodami v AM');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% 3. Zaverecny odhad trendove slozky po separaci sezonni slozky
%---------------------------------------------------------------------------
% Provedeme pro kazdy odhad sezonni slozky dle a) az d), kterou odstranime
  close all
  xlog_sz1 = xlog - szlog1;
  xlog_sz2 = xlog - szlog2;
  xlog_sz3 = xlog - szlog3;
  xlog_sz4 = xlog - szlog4;
% A nyni nalezneme prislusne vysledne odhady trendu metodou klouzavych
% prumeru s delkou okna d=12, coz zajisti odstraneni pripadnych zbytkovych
% komponent s touto periodou.
  [trlog1,trlogE1]=polyfilt(2,d,xlog_sz1,t,[],[],figure);
pause;
  [trlog2,trlogE2]=polyfilt(2,d,xlog_sz2,t,[],[],figure);
pause;
  [trlog3,trlogE3]=polyfilt(2,d,xlog_sz3,t,[],[],figure);
pause;
  [trlog4,trlogE4]=polyfilt(2,d,xlog_sz4,t,[],[],figure);
pause;
%---------------------------------------------------------------------------
% Prubeh vsech odhadu trendu:
  figure; plot(t,trlog1,'r',t,trlog2,'g',t,trlog3,'b',t,trlog4,'w');
pause;
%---------------------------------------------------------------------------
% 3. Analyza rezidualni chybove slozky
%---------------------------------------------------------------------------
% Rezidualni chybove komponenty pro vsechny metody:
  close all
  elog1 = xlog - szlog1 - trlog1;
  trand(elog1)                        % Test nahodnosti
  elog2 = xlog - szlog2 - trlog2;
  trand(elog2)                        % Test nahodnosti
  elog3 = xlog - szlog3 - trlog3;
  trand(elog3)                        % Test nahodnosti
  elog4 = xlog - szlog4 - trlog4;
  trand(elog4)                        % Test nahodnosti
  figure; plot(t,elog1,'r',t,elog2,'g',t,elog3,'b',t,elog4,'w');
% Prubehy trendu i rezidualni chybove slozky dle metody b) (zelena) a
% d) (bila) vykazuji vynikajici shodu.
% Prubeh chyby metody c) (modra) vykazuje evidentne periodicke chovani,
% coz potvrzuje i nejhorsi vysledek testu nahodnosti (2 x odmitnuto)
pause;
%---------------------------------------------------------------------------
% Autokorelacni funkce:
  elog1ro=acf(elog1,figure);
% Neklesa dostatecne rychle k nule => metoda a) neodstranila dokonale trend
pause;
  elog2ro=acf(elog2,figure);
% Metoda b) vyhovujici
pause;
  elog3ro=acf(elog3,figure);
% Pomaly pokles a oscilacni charakter + velka neurcitost dana rozsirujicim
% se pasem spolehlivosti => metoda c) zcela nevyhovujici
pause;
  elog4ro=acf(elog4,figure);
% Metoda d) vyhovujici
pause;
%---------------------------------------------------------------------------
% Parcialni autokorelacni funkce:
  elog1pacf=pacf(elog1ro,n,figure);
pause;
  elog2pacf=pacf(elog2ro,n,figure);
pause;
  elog3pacf=pacf(elog3ro,n,figure);
pause;
  elog4pacf=pacf(elog4ro,n,figure);
pause;
%---------------------------------------------------------------------------
% 4. ZAVER
%---------------------------------------------------------------------------
% Metody b) a d) davaji verohodne a prakticky shodne vysledky:
  close all
  figure; plot(t,[szlog2,szlog4]);
pause;
  figure; plot(t,[trlog2,trlog4]);
pause;
  figure; plot(t,[elog2,elog4]);
pause;
% Pro zaverecnou prezentaci vysledku pouzijeme metodu d), ktera na
% rozdil od b) poskytuje i konfidencni pasy a navic vsechny testy
% potvrzuji nahodne chovani.
% Zpetne transformujeme odlogaritmovanim
  szlog  = szlog4;     sz  = exp(szlog);
  szlogE = pE4(:,2:3); szE = exp(szlogE);
  trlog  = trlog4;     tr  = exp(trlog);
  trlogE = trlogE4;    trE = exp(trlogE);
  elog   = elog4;      e   = exp(elog);
% a vykreslime
% tr s konfidencnimi pasy spolu s puvodnimi daty:
% --------------------------------------------------
  figure;
  plot(t,x,'go',t,tr,'w-',...
       t,[tr.*trE(:,1),tr./trE(:,1)],'y:',...
       t,[tr.*trE(:,2),tr./trE(:,2)],'r:');
  title('Data x (o), tr (-)');
  xlabel(xinf(2,:));
  trand(e)   % Opet plne potvrzeno nahodne chovani i po odlogaritmovani
pause;
%---------------------------------------------------------------------------
% sz + e s konfidencnimi pasy:
  figure;
  plot(t,sz,'w-',...
       t,[sz.*szE(:,1),sz./szE(:,1)],'y:',...
       t,[sz.*szE(:,2),sz./szE(:,2)],'r:');
  title('sz (-)');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% Rezidualni slozka e v MM:
  figure;
  plot(t,e,'y-');
  title('Rezidualni slozka e v MM');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% Rezidualni slozka elog v transformovanem AM:
  figure;
  plot(t,elog,'y-');
  title('Rezidualni slozka elog v transformovanem AM');
  xlabel(xinf(2,:));
pause;
%---------------------------------------------------------------------------
% AUTOKORELACNI FUNKCE REZIDUALNI CHYBOVE SLOZKY:
  [ero,eroE,mi,miE,eacf] = acf(elog,figure);
  xlabel(xinf(2,:));
  [mi-miE(1),mi+miE(1)]    % 95% interval spolehlivosti pro stredni hodnotu
  [mi-miE(2),mi+miE(2)]    % totez presneji za predpokladu normality chyb
% Oba intervaly jsou uzke a zahrnuji nulu, tj. nezamitame hypotezu o nulove
% stredni hodnote.
% Varovne hlaseni vsak upozornuje na moznou nestacionaritu.
% V grafu acf lze skutecne stale pozorovat slabe oscilatoricke chovani, ktere
% upozornuje na neuplne odstraneni periodicke slozky.
pause;
%---------------------------------------------------------------------------
% Pro kontrolu odhadneme jeste mi a miE linearni regresi stupne 0:
  [y,E,r,s,bE] = polyfilt(0,n,elog,[],[],[],figure);
  title('Regrese konstantou');
  xlabel(xinf(2,:));
  bE(1)                      % stredni hodnota
  [bE(1)-bE(2),bE(1)+bE(2)]  % konfidencni interval
% Je ponekud sirsi nez odhad poskytovany pomoci acf, ktery je za predpokladu
% stacionarity presnejsi. Nula je vsak rovnez zahrnuta.
pause;
%---------------------------------------------------------------------------
% PARCIALNI AUTOKORELACNI FUNKCE REZIDUALNI CHYBOVE SLOZKY:
  [epacf,phi,phiE,v] = pacf(eacf,n,figure);
  xlabel(xinf(2,:));
  v                          % Odhad rozptylu generujiciho bileho sumu
% Oba posledni grafy naznacuji moznost modelovat elog jako AR proces radu 11,
% nebot hodnoty pacf(h) pro h>11 se nachazi v nulovem konfidencnim pasu.
% jejich chovani v nulovem konfidencnim pasu se vsak nezda byt nahodne.
% Vzhledem k vysokemu radu AR pujde asi spise o ARMA.
% To spolu s vyse vypsanym varovanim naznacuje nutnost dalsiho podrobnejsiho
% zkoumani a hledani adekvatniho modelu pro chybovou slozku elog.
% Pro tento ucel je vsak nezbytny specialni software jako napr.
% SYSTEM IDENTIFICATION TOOLBOX pro MATLAB.
%
% Vychozi data a vysledek dekompozice ulozime do souboru airpass.mat:
  save airpass xinf x szlog szlogE sz szE trlog trlogE tr trE elog e

%================================== KONEC ==================================
echo off
