%=================================================================
% POUZITI LINEARNIHO REGRESNIHO MODELU PRI DEKOMPOZICI CASOVE RADY
%=================================================================

% (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
us_deaths=[crpath,'/brockwel.dat/deaths.dat'];
[x,xinf]=getdata(us_deaths);   % nacteni a vykresleni
n=length(x)
t=(1:n).';
pause;
%---------------------------------------------------------------
% Nalezneme vyznamne periodicke komponenty pomoci periodogramu a
% Fisherova-Siegelova testu:
  [pdg,kF,kS]=tperiod(x,figure);  % periodogram pdg a jeho graf
  kF        % indexy vyznamnych harmonickych slozek uzitim Fisherova testu
  kS        % indexy vyznamnych harmonickych slozek uzitim Siegelova testu
% Nejvyraznejsi slozka kF(1)=kS(1)=6 (6-ta harmonicka) indikuje periodu
% 12 mesicu sezonni slozky.
pause;
%---------------------------------------------------------------
% Zlepsime citlivost testovani odectenim parabolickeho trendu:
  D=[t.^2 ,t,ones(n,1)];  % Cast matice planu pro parabolicky trend
  [ans,ans,ans,ans,yTr]=linregr(x,D(:,1:3),D(:,1:3));
% Jednoduseji totez uzitim 'polfs' - viz dale
  xTr=polfs(x,2,[],t);
  max(abs(xTr(:,1)-yTr(:,1)))   % Kontrola shody obou vysledku
% Fisherova-Siegelova testu:
  [pdg,kF,kS]=tperiod(x-xTr(:,1),figure);  % periodogram pdg a jeho graf
  kF        % indexy vyznamnych harmonickych slozek uzitim Fisherova testu
  kS        % indexy vyznamnych harmonickych slozek uzitim Siegelova testu
% Nejvyraznejsi slozka kF(1)=kS(1)=6 (6-ta harmonicka) indikuje periodu
% 12 mesicu sezonni slozky.
pause;
%--------------------------------------------------------------
% do Fourierovy rady v linearnim modelu zahrneme nejprve jen
% nejvyznamnejsi harmonicke komponenty 6,12,18,30:
  w = 2*pi*[6,12,18,30]./n;   % uhlove rychlosti vybranych harm. komponent
  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
  D=[D,cos(t*w),sin(t*w)]  ; % Pridame do D cast pro Fourierovu radu
  C=[C,cos(tt*w),sin(tt*w)]; % Pridame do C cast pro Fourierovu radu
  [bE,rE,sE,PbE,CbE]=linregr(x,D,C); % Provedeme regresi
  tr=1:3;         % indexy koeficientu paraboly trendu
  sz=4:size(C,2); % indexy Fourier. koeficientu a,a,,..; b,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;
%---------------------------------------------------------------
% Vyse uvedeny postup je pohodlneji implementovan v procedure 'polfs',
% kterou nyni pouzijeme pro dekompozici s indexy kS, kterych je vice nez v kF
  dt = 1/12;
  t  = (1973:dt:1979-dt).';      % Casovou osu definujeme dle skutecnosti
  tt = linspace(1973,1979-dt).';
  [yE,trE,pE,r,s,bE,cE,sE]=polfs(x,2,n*dt./kS,t,tt,[],figure);
  bE,       % parametry parabolickeho trendu + 95% polosirka intervalu spol.
  cE,       % koeficienty u kosinu ve F.r. + 95% polosirka intervalu spol.
  sE,       % koeficienty u sinu   ve F.r. + 95% polosirka intervalu spol.
pause;
%---------------------------------------------------------------
% Alternativne muzeme pro dekompozici pouzit proceduru 'polper', ktera
% vsak umoznuje predikci na siti tt pouze pro trend.
  [yE,trE,pE,r,s,bE]=polper(x,2,n/kS(1),t,tt,[],figure);
  bE,       % parametry parabolickeho trendu + 95% polosirka intervalu spol.
pause;
%---------------------------------------------------------------
% Porovname graficky vysledky dekompozice z obou procedur na siti t:
  [yE1,trE1,pE1,r1,s1,bE1,cE1,sE1]=polfs(x,2,n*dt./kS,t);
  [yE2,trE2,pE2,r2,s2,bE2]=polper(x,2,n/kS(1),t);
  bE1, bE2,            % Porovnejte koeficienty parabolickeho trendu
pause;
  figure; plot(t,[yE1(:,1),yE2(:,1)]); title('Porovnani vyhlazeni');
pause;
  figure; plot(t,[trE1(:,1),trE2(:,1)]); title('Porovnani trendu');
pause;
  figure; plot(t,[pE1(:,1),pE2(:,1)]); title('Porovnani periodicke slozky');
pause;
  figure; plot(t,[r1(:,1),r2(:,1)]); title('Porovnani chybove slozky');
%---------------------------------------------------------------------------
% PROSTUDUJTE ALGORITMUS PROCEDUR 'polfs' a 'polper' pomoci prikazu 'type'.
%================================== KONEC ==================================
echo off
