%=================================================================
%       FILTRACE POMOCI METODY EXPONENCIALNIHO VYROVNAVANI
%              (PROCEDURY EXP1FILT a EXP2FILT)
%=================================================================

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

format compact;
echo on
%
% Exponencialni vyrovnavani = vyhlazeni a predikce pomoci vazene regrese
% minulych hodnot polynomem zadaneho stupne.
% Pro tento ucel jsou k dispozici procedury:
% exp1filt ... jednoduche exponencialni vyrovnavani (s polynomem stupne 0)
% exp2filt ... dvojite exponencialni vyrovnavani    (s polynomem stupne 1)
%-----------------------------------------------------------------
% PROSTUDUJTE OVLADANI TECHTO PROCEDUR POMOCI PRIKAZU
% help exp1filt a help exp2filt
% (vypis algoritmu v tele procedury prikazem type misto help)
%-----------------------------------------------------------------
keyboard;
%=================================================================
% VYGENERUJEME SIMULACNI DATA JAKO V crcv7 a crcv8:
  n = 100;
  t = (0:n-1).';
  y = 2 + sin(2*pi*t./n);  % Sinusovy trend (1 perioda)
  sigma = 0.2;
  load sinus  % =totez jako: x = y + sigma*randn(n,1);
  figure; plot(t,x,'w'); hold on; plot(t,y,'r','LineWidth',2);
  title('SIMULACE SINUSOVKY KONTAMINOVANE BILYM SUMEM');
pause;
%---------------------------------------------------------------
% Vyzkousime nejprve jednoduche vyrovnavani s potlacenim automaticke
% korekce vyhlazovaciho parametru alfa (ACS=Inf):
  [y,E,r,S] = exp1filt(x,[],Inf,t);
pause;
%---------------------------------------------------------------
% Pocatecni odhad optimalni hodnoty alfa se dostal mimo interval
% [0.7,), je tedy nutno zvysit stupen polynomu o 1 a prejit na
% dvojite exponencialni vyrovnavani.
% Pouzijeme nyni standardni korekci alfa (implicitni hodnota ACS=5)
  [y,E,r,S,b] = exp2filt(x,[],[],t);
  S      % vypis pocatecnich (1.radek) a koncovych stavovych hodnot (2.radek):
         % 1.sloupec S(:,1) = hodnoty alfa
% Procedura hlasi korekci alfa. Pokud chceme detailne zachytit prubeh vsech
% stavovych promennych, zopakujeme vypocet po krocich v cyklu.
pause;
%---------------------------------------------------------------
  y1 = []; E1 = []; r1 = []; b1 = [];
  Sk = S(1,:);    % nastavime pocatecni stav
  for k=1:n
      [yk,Ek,rk,Sk,bk] = exp2filt(x(k),Sk,[],t(k));
      y1(k)=yk; E1(k)=Ek; r1(k,:)=rk; b1(:,k)=bk;
echo off
  end
echo on
  figure;  plot(t,[y,y1'],'-',t,[y-E,y+E],'y:',t,[y-E1',y+E1'],'g:');
pause;
%---------------------------------------------------------------
% Vypoctene hodnoty jsou stejne
  max(abs(y-y1.'))
  max(abs(E-E1.'))
  max(abs(r-r1))
  max(abs((b-b1).'))
pause;
%---------------------------------------------------------------
% Stavova matice S vsak jiz zachycuje stavove vektory celeho prubehu, nebot
% pri kazdem volani se do S pripisuje stav po vyhlazeni jedne hodnoty:
  size(Sk)
  Sk(1,1)                          % Pocatecni odhad alfa
  I=find(Sk(2:n,1)~=Sk(1:n-1,1));  % Indexy kroku, v nichz nastala zmena alfa
  [I,Sk(I,1)]                      % Vypis techto indexu a prislusnych alfa
pause;
%---------------------------------------------------------------
% Zopakujeme vypocet s vykreslenim prubehu vyhlazeni:
  [y,E,r,S,b] = exp2filt(x,[],[],[],[],[],figure);
  yp1 = x-r(:,2);  % Prubeh predikce y(k) z dat x(t) pro t<=k-1
                   % totiz r(:,2) = x-yp1 je chyba predikce
% Prikreslime predikce do existujiciho grafu tucnou carou:
  hold on; plot(t,yp1,'LineWidth',1.5); hold off
pause;
%---------------------------------------------------------------
% Postupne budeme nyni snizovat hodnotu ACS, cimz bude umoznena
% castejsi korekce alfa:
  [y,E,r,S,b] = exp2filt(x,[],4,[],[],[],gcf);
pause;
%---------------------------------------------------------------
  [y,E,r,S,b] = exp2filt(x,[],3,[],[],[],gcf);
pause;
%---------------------------------------------------------------
  [y,E,r,S,b] = exp2filt(x,[],2,[],[],[],gcf);
% Vidime, ze castejsi korekce alfa zvysuji riziko vybehnuti
% mimo povoleny interval [0.7,1). Doporucene hodnoty ACS=4 az 6.
pause;
close all
%===============================================================
% Dalsi demonstrace bude jiz s realnymi daty prevzatymi z Ciprovy
% ucebnice: Analyza dat s aplikacemi v ekonomii:
  [x,xinf]=getdata([crpath,'/cipra.dat/zakazky.dat']);
  n=length(x);
  x1=x(1:22);   % Prvy pozorovany usek
% zpracujeme jednododuchym exponencialnim vyrovnanim:
  [y1,E1,r1,S1] = exp1filt(x1,[],[],[],[],[],figure);
pause;
%---------------------------------------------------------------
% Provedeme nyni predikci o 3 kroky dopredu:
  tt = 23:25;
  [y1,E1,r1,S1] = exp1filt(x1,[],[],[],tt);
  y1
% Vsimnete si, ze predikce je konstantni, nebot delame regresi
% polynomem stupne 0.
pause;
%---------------------------------------------------------------
% Dokoncime vypocet, jakmile mame dalsi pozorovana data:
  x2=x(23:n);   % Druhy pozorovany usek
  [y2,E2,r2,S2] = exp1filt(x2,S1);
  y2(1:3)        % Toto jsou zpresnene hodnoty, predpovedene v y1
pause;
close all
%---------------------------------------------------------------
% Vezmeme nyni jiny soubor dat
  [x,xinf]=getdata([crpath,'/cipra.dat/odbyt.dat']);
  n=length(x)
  x1=x(1:24);   % Prvy pozorovany usek
% zpracujeme opet jednododuchym exponencialnim vyrovnanim:
  [y1,E1,r1,S1] = exp1filt(x1,[],[],[],[],[],figure);
% Pri pocatecnim odhadu (i dale pri vypoctu) se optimalni alfa
% dostava mimo povoleny rozsah. To signalizuje nutnost zvysit
% stupen polynomu.
% I z obrazku vidime, ze vyhlazeny prubeh nejde stredem dat
% Pouzijeme tedy dvojite exponencialni vyrovnani:
pause;
%---------------------------------------------------------------
  [y1,E1,r1,S1,b] = exp2filt(x1,[],[],[],[],[],figure);
pause;
%---------------------------------------------------------------
% Provedeme nyni predikci o 3 kroky dopredu:
  tt = 25:27;
  [y1,E1,r1,S1,b] = exp2filt(x1,[],[],[],tt);
  y1
% Vsimnete si, ze predikce je linearni, nebot delame regresi
% polynomem stupne 1.
pause;
%---------------------------------------------------------------
% Dokoncime vypocet, jakmile mame dalsi pozorovana data:
  x2=x(25:n);    % Druhy pozorovany usek
  [y2,E2,r2,S2,b] = exp2filt(x2,S1);
  y2(1:3)        % Toto jsou zpresnene hodnoty, predpovedene v y1
pause;
%---------------------------------------------------------------
% Zaverem ukazka predikce o 5 hodnot do minulosti i do budoucnosti
% z celeho pozorovani:
  tt = -4:n+5;
  [y,E,r,S,b] = exp2filt(x,[],[],[],tt,[],figure);
%================================== KONEC ==================================
echo off
