%=================================================================
%               ODHAD STUPNE POLYNOMIALNIHO TRENDU
%          NEPARAMETRICKE METODY PRO ODHAD SEZONNI SLOZKY
%=================================================================

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

% ANALYZA SOUBORU deaths.dat:
% Mesicni pocty umrti pri automobilovych nehodach v USA
% x=absolutni pocet umrti
% t=1,2,...,72: 1973 az 1978
%==============================================================
echo on
format compact
format long
us_deaths=[crpath,'/brockwel.dat/deaths.dat'];
[x,xinf]=getdata(us_deaths);   % nacteni a vykresleni
n=length(x)
t=(1:n).';
pause;
%---------------------------------------------------------------
% Jiz znamym postupem identifikujeme periodicke komponenty napr.
% pomoci Fisherova testu:
  xTr=polfs(x,2,[],t);
  [pdg,kF]=tperiod(x-xTr(:,1),figure);  % periodogram pdg a jeho graf
  kF        % indexy vyznamnych harmonickych slozek uzitim Fisherova testu
% Nejvyraznejsi slozka kF(1)=kS(1)=6 (6-ta harmonicka) indikuje periodu
% 12 mesicu sezonni slozky.
% pripadnou zakladni periodu kF(#)=1 budeme ignorovat:
% kF(#)=[]
pause;
%--------------------------------------------------------------
% Pro dekompozici pouzijeme proceduru 'polfs'.
% Predpokladejme nejprve stupen 0 (konstantni trend):
  [yE,trE,pE,r,s,bE]=polfs(x,0,n./kF,[],[],[],figure);
  bE,       % hodnota konstatniho trendu + 95% polosirka intervalu spol.
% Rezidualni chybova slozka nevykazuje nahodne oscilace kolem nuly;
% stupen je tedy nutno zvetsit - predpokladejme linearni trend:
pause;
%---------------------------------------------------------------
  close all
  [yE,trE,pE,r,s,bE]=polfs(x,1,n./kF,[],[],[],figure);
  bE,       % koeficienty linearniho trendu + 95% polosirka intervalu spol.
% 95% interval spolehlivosti pro koeficient u prvni mocniny je:
  [bE(1,1)-bE(1,2),bE(1,1)+bE(1,2)]
% V tomto intervalu nelezi hodnota 0, takze s rizikem 5% zamitame
% hypotezu o nulovosti tohoto koeficientu a tedy predpoklad konstantniho
% trendu;
% To potvrzuje i prubeh rezidualni chybove slozky, ktera stale nevykazuje
% nahodne oscilace kolem nuly.
% Predpokladejme proto parabolicky trend:
pause;
%---------------------------------------------------------------
  close all
  [yE,trE,pE,r,s,bE]=polfs(x,2,n./kF,[],[],[],figure);
  bE,       % koeficienty linearniho trendu + 95% polosirka intervalu spol.
% 95% interval spolehlivosti pro koeficient u druhe mocniny je:
  [bE(1,1)-bE(1,2),bE(1,1)+bE(1,2)]
% V tomto intervalu nelezi hodnota 0, takze s rizikem 5% zamitame
% hypotezu o nulovosti tohoto koeficientu a tedy take predpoklad
% linearniho trendu;
% Protoze rezidualni chybova slozka jiz nahodne osciluje kolem nuly,
% lze se duvodne domnivat, ze trend je parabolicky.
% Pro kontrolu zkusime kubicky trend:
pause;
%---------------------------------------------------------------
  close all
  [yE,trE,pE,r,s,bE,cE,sE]=polfs(x,3,n./kF,[],[],[],figure);
  bE,       % koeficienty linearniho trendu + 95% polosirka intervalu spol.
% 95% interval spolehlivosti pro koeficient u treti mocniny je:
  [bE(1,1)-bE(1,2),bE(1,1)+bE(1,2)]
% V tomto intervalu lezi hodnota 0, takze nezamitame
% hypotezu o nulovosti tohoto koeficientu a definitivne prijimame
% predpoklad parabolickeho trendu.
pause;
close all;
%---------------------------------------------------------------
% Jina ponekud mene exaktni metoda odhadu stupne polynomialniho
% trendu je zalozena na diferencovani.
% Je pouzitelna pouze, pokud data obsahuji nejvyse jednu periodickou slozku
% o zname zakladni periode (zpravidla sezonni slozku).
% 1.krok:
% Jestlize data obsahuji periodickou slozku, odstranime ji diferencovanim
% na vzdalenost rovnajici se jeji zakladni periode d;
% pro nas soubor je
  d=12;
  xd12 = x(d+1:n)-x(1:n-d);
  figure; plot(xd12);
  title('Prubeh po odstraneni sezonni slozky diferencovanim');
% Timto diferencovanim jsme vsak snizili o 1 i stupen polynomialniho trendu.
% Chybovy prubeh neosciluje kolem konstantni hodnoty.
% Trend tedy nemohl byt ani konstantni (oscilace by musely byt kolem nuly),
% ani linearni (oscilace by musely byt kolem nenulove konstantni hodnoty).
pause;
%---------------------------------------------------------------
% 2.krok:
% Snizime stupen dalsim diferencovanim, tentokrat jiz standardne na
% vzdalenost 1:
  nd12=length(xd12);
  xd12d1=x(2:nd12)-x(1:nd12-1);
  figure; plot(xd12d1);
  title('Prubeh po dvojim diferencovani');
% Chybovy prubeh jiz osciluje kolem konstantni hodnoty.
% Trendovou komponentu proto muzeme pokladat za parabolickou.
% Kdyby tomu tak nebylo, zopakovali bychom krok 2 atd.
% V pripade, ze data neobsahuji sezonni slozku, pak krok 1 samozrejme
% vypustime a zacneme provadet rovnou krok 2.
pause;
close all;
%===========================================================================
% Pro separaci sezonni slozky pouzijeme nyni metodu maleho trendu
% implementovanou v procedure 'szsmtr':
  AM = 0; d=12;                      % 0=aditivni model; 12=perioda
  [sz,tr] = szsmtr(x,12,AM,figure);  % Metoda maleho trendu
pause;
%---------------------------------------------------------------
% Alternativne lze pouzit presnejsi metodu klouzaveho prumeru
% implementovanou v procedure 'szma':
  [sz,tr,msz] = szma(x,12,AM,figure); % Metoda klouzaveho prumeru
  msz,       % Aditivni korekce trendu je zanedbatelna
pause;
close all;
%===============================================================
% Predchozi metody lze uzit i primo pro pripad multiplikativniho modelu.
% Pro tento ucel uzijeme
% SOUBOR: airpass.dat
% Cestujici v mezinarodni letecke doprave
% x=mesicni souhrny v tisicich cestujicich
% t=1,2,...,144: leden 1949 az prosinec 1960
%===============================================================
  air_pass=[crpath,'/brockwel.dat/airpass.dat'];
  [x,xinf]=getdata(air_pass);   % nacteni a vykresleni
  n=length(x);
  [l,y]=powtr(x,4);
pause(2);
  figure; plot(y);
  title('Po logaritmicke transformaci');
  xlabel(xinf(2,:));
% Zavislost rozptylu na trendu je linearni, takze lze predpokladat,
% ze se data ridi multiplikativnim modelem.
pause;
%---------------------------------------------------------------
% Zjistime periodu sezonni slozky:
  [pdg,kF]=tperiod(y,figure);  % Model musi byt aditivni, tj. uzijeme y
                               % po logaritmicke transformaci.
  pause(2)
% Presneji odectenim parabolickeho trendu:
  yTr=polfs(y,2,[],(1:n).');
  [pdg,kF]=tperiod(y-yTr(:,1),figure);  
  kF
% kF(1)=12 je harmonicka slozka odpovidajici rocni sezonni periode,
% tj. za periodu zvolime opet 12 mesicu.
pause;
%---------------------------------------------------------------
  MM = 1; d=12;                      % 1=multiplikativni model; 12=perioda
  [sz,tr] = szsmtr(x,12,MM,figure);  % Metoda maleho trendu
pause;
%---------------------------------------------------------------
  [sz,tr,msz] = szma(x,12,MM,figure);    % Metoda klouzaveho prumeru
  msz,       % Multiplikativni korekce trendu je zanedbatelna
%---------------------------------------------------------------
% PROSTUDUJTE ALGORITMUS PROCEDUR 'szsmtr' a 'szma' pomoci prikazu 'type'.
%================================== KONEC ==================================
echo off
