%==============================================================
% 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
us_deaths=[crpath,'/brockwel.dat/deaths.dat'];
[x,xinf]=getdata(us_deaths);   % nacteni a vykresleni
n=length(x);
pause;
%---------------------------------------------------------------
xacf=acf(x);
plotdata(xacf); % vypocet a vykresleni autokorelacni funkce
pause;

%---------------------------------------------------------------
% Separace sezonni slozky s periodou 12 mesicu
%---------------------------------------------------------------
% Odhad sezonni slozky pomoci metody maleho trendu:
[xsz,xtr]=szsmtr(x,12);
plotdata(xsz,xinf); ax=axis;
title('Metoda maleho trendu: sezonni slozka');
pause;
%---------------------------------------------------------------
% Graf: Puvodni data + data ocistena od sezonni slozky + trend
plotdata([x,x-xsz,xtr],xinf);
title('Metoda maleho trendu: x, x-sz, tr');
pause;
%---------------------------------------------------------------
% Chybova komponenta:
xe=x-xsz-xtr;
plotdata(xe,xinf); axis(ax);
title('Metoda maleho trendu: x-sz-tr');
pause;
%---------------------------------------------------------------
% Odhad sezonni slozky pomoci metody klouzaveho prumeru:
[xsz,xtr]=szma(x,12);
plotdata(xsz,xinf); ax=axis;
title('Metoda klouzaveho prumeru: sezonni slozka');
pause;
%---------------------------------------------------------------
% Graf: Puvodni data + data ocistena od sezonni slozky + trend
plotdata([x,x-xsz,xtr],xinf);
title('Metoda klouzaveho prumeru: x, x-sz, tr');
pause;
%---------------------------------------------------------------
% Chybova komponenta:
xe=x-xsz-xtr;
plotdata(xe,xinf); axis(ax);
title('Metoda klouzaveho prumeru: x-sz-tr');
pause;
%---------------------------------------------------------------
% Autokorelacni funkce chybove komponenty:
plotdata(acf(xe),xinf);
pause;
%---------------------------------------------------------------
% Odhad trendu a sezonni slozky pomoci regresniho modelu
% s Fourierovou radou
%---------------------------------------------------------------
% Nalezneme periodicke komponenty pomoci periodogramu:
  n2=fix(n./2);
  X=fft(x);          % Diskretni spektrum Fourierovych koeficientu
  mean(x)-X(1)./n    % X(1)/n = prumer
  Xs=X(2:n);         % Zbytek spektra=Harmonicke komponenty
  Xs(n2-5:n2+5)      % Jsou  komplexne sdruzene kolem stredu
pause;
%---------------------------------------------------------------
  figure;
  stem(2*abs(Xs))    % Graf amplitud: stredova symetrie je tez zretelna
pause;
%---------------------------------------------------------------
  figure;
  stem(180*atan2(imag(Xs),real(Xs))./pi) % Graf fazi: stredova antisymetrie
pause;
%---------------------------------------------------------------
  pdg=(2./n)*abs(Xs(1:n2)).^2;  % periodogram
  figure; stem(pdg)             % a jeho graf
  harm=find(pdg>0.15e7)         % Indexy vyznamnych harmonickych komponent:
% o periodach:
  n./harm
% 1.harmonicka = zakladni perioda dat=n => ignorujeme
% Nejvyraznejsi perioda 12 (6-ta harmonicka) indikuje periodu sezonni slozky.
pause;
%---------------------------------------------------------------
% do Fourierovy rady v linearnim modelu zahrneme proto jen
% harmonicke komponenty 6,12,18,30:
  w = 2*pi*harm(3:6)./n;   % uhlove rychlosti vybranych harm. komponent
  t=(1:n).';
  tt=linspace(1,n).'; N=length(tt);
  D=[t.^2 ,t,ones(n,1)];  % Cast matice planu pro parabolicky trend
  C=[tt.^2,tt,ones(N,1)]; % Cast matice param. funkci pro trend
  for wh=w.',
      D=[D,cos(wh.*t),sin(wh.*t) ] ; % Pridame do D cast pro wh
      C=[C,cos(wh.*tt),sin(wh.*tt)]; % Pridame do C cast pro wh
  end
  [bE,rE,sE,PbE,CbE]=linregr(x,D,C); % Provedeme regresi
  tr=1:3;         % indexy koeficientu paraboly trendu
  sz=4:size(C,2); % indexy Fourierovych koeficientu a,b,a,b,.. harmon. slozek
  bE(tr,:)        % a jejich odhady spolu s polosirkou intervalu spolehlivosti
  bE(sz,:)
pause;
%---------------------------------------------------------------
  figure;
  plot(tt,CbE(:,1),'y',tt,[CbE(:,1)-CbE(:,2),CbE(:,1)+CbE(:,2)],'r-');
  title('trend + sezonni slozka');
pause;
%---------------------------------------------------------------
  xtrf=C(:,tr)*bE(tr,1);   % Jen trend
  figure; plot(tt,xtrf);
  title('trend');
pause;
%---------------------------------------------------------------
  xszf=C(:,sz)*bE(sz,1);   % Jen sezonni slozka
  figure; plot(tt,xszf);
  title('sezonni slozka');
pause;
%---------------------------------------------------------------
% Porovnejte s odhadem klouzavymi prumery na obr. 3
  figure(3);
%================================== KONEC ==================================
echo off
