%=================================================================
%                       TESTOVANI NAHODNOSTI
%=================================================================

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

format compact;
echo on
%
% Pomoci testu nahodnosti overujeme, zda data maji nahodny charakter,
% presneji, zda jsou nezavisla a stejne rozlozena (IID).
% Je mozno je pouzit i pro testovani nahodnosti rezidualni chybove
% komponenty e_t po provedeni dekompozice, i kdyz v tomto pripade musime
% pocitat s tim, ze proces dekompozice zanesl do e_t
% korelacni zavislosti, tj. porusil predpoklad IID.
% Pokud testy nepotvrdi nahodnost, pak bud dekompozice nebyla dokonala
% (v datech jsou relikty trendu, pripadne periodickych a cyklickych
% slozek - proverime) a nebo zbytkovy sum muze jiz byt stacionarni
% avsak se silnymi korelacemi (autokorelacni funkce klesa pomalu k nule).
% V takovem pripade hledame pro e_t vhodny model (zpravidla typu ARMA)
% uzitim Box-Jenkinsovy metodologie (viz crcv11).
%
% Pro testovani nahodnosti je k dispozici procedura 'trand', v niz je
% implementovano 5 testu nahodnosti zpracovanych podle monografie
% T. Cipra: Analyza casovych rad s aplikacemi v ekonomii,
% SNTL/ALFA Praha 1986, kap. 10 (str.94-99)
% Testy c.1,3,4 maji byt citlivejsi na pritomnost (linearniho) trendu,
% zatimco zbyvajici dva na periodicke komponenty.
% Test c.1 (pocet znamenek diferenci) rozhodne nebude reagovat na
% periodicky se opakujici zmeny.
%-----------------------------------------------------------------
% PROSTUDUJTE OVLADANI, RESP. ALGORITMUS PROCEDURY trand POMOCI
% PRIKAZU: help trand, resp. type trand
%-----------------------------------------------------------------
keyboard;
%=================================================================
% PROVERIME FUNKCI PROCEDURY trand NEJPRVE NA TESTOVACICH DATECH
% (Cipra: tab.10.1 a 10.2)
%=================================================================
  x1 = getdata([crpath,'/cipra.dat/tab10_1.dat']); % Data z tab. 10.1
  [x1trslt,x1tstat,x1tvalue,x1tname]=trand(x1);
% nebo
  [x1trslt,x1tstat,x1tvalue,x1tname]=trand(x1,0);  % potlaci vypis vysledku
% Vysledky testu (xtrslt), hodnoty statistik (xtstat) a hlavnich pocitanych
% velicin (xtvalue) jednotlivych testu pak muzeme zobrazit ve forme
% tabulky (sloupce odpovidaji jednotlivym testum), napr. takto:
  [x1trslt;x1tstat;x1tvalue]
pause;
%---------------------------------------------------------------
% nebo podrobneji
  [x1tname,num2str([x1trslt;x1tstat;x1tvalue].')] % jen pro ver.5
% Vidime, ze vysledky vsech testu odpovidaji numerickym hodnotam
% z Ciprovy monografie
pause;
%---------------------------------------------------------------
% Podobne proverime testovaci soubor z tab. 10.2:
  x2 = getdata([crpath,'/cipra.dat/tab10_2.dat']);  % Data z tab. 10.2
  [x2trslt,x2tstat,x2tvalue,x2tname]=trand(x2,0);
  [x2tname,num2str([x2trslt;x2tstat;x2tvalue].')]   % jen pro ver.5
%  [x2trslt;x2tstat;x2tvalue]                        % alternativa pro ver.4
pause;
close all;
%=================================================================
% PROVERIME NYNI CITLIVOST TESTU NA PRITOMNOST PERIODICKE SLOZKY
%=================================================================
  n = 200;
  t = (0:n-1).';
  mi = 2;               % nenulova stredni uroven
  wn = randn(n,1);      % Bily sum WN(0,1)
  xp = sin(4*pi*t./n);  % 2 periody sinusovky
pause;
%---------------------------------------------------------------
% 1. Cisty bily sum:
  x = wn;
  figure; plot(t,x); title('Bily sum s nulovou stredni hodnotou');
  trand(x)
  x = mi + wn;
  plot(t,x); title('Bily sum WN(2,1) s nenulovou stredni hodnotou');
  trand(x)
% Zadny test nepopira podle ocekavani nahodnost bez ohledu na velikost
% stredni hodnoty
pause;
%---------------------------------------------------------------
% 2. Bily sum + sinusovka s amplitudou 0.1*sigma az 1*sigma:
  figure;
  for A=0.1:0.1:1
      x = mi + A*xp + wn;
      plot(t,x);
      title(['Bily sum WN(2,1) + sinusovka s amplitudou A=',num2str(A)]);
      A,trand(x)
pause;
echo off
  end
echo on
% Posledni 3.testy reaguji dobre.
pause;
%---------------------------------------------------------------
% 3. Bily sum + sinusovka s amplitudou 5*sigma az 10*sigma:
  figure;
  for A=5:10
      x = mi + A*xp + wn;
      plot(t,x);
      title(['Bily sum WN(2,1) + sinusovka s amplitudou A=',num2str(A)]);
      A,trand(x)
pause;
echo off
  end
echo on
% 1. a 2. test nereaguji vubec.
pause;
close all;
%=================================================================
% PROVERIME NYNI CITLIVOST TESTU NA PRITOMNOST LINEARNIHO TRENDU
%=================================================================
% 4. Bily sum + primka se zdvihem h=xtr(n)=0.1*sigma az 1*sigma:
  xtr = t-n./2;
  figure;
  for h=0.1:0.1:1
      x = mi + (h./xtr(n))*xtr + wn;
      plot(t,x);
      title(['Bily sum WN(2,1) + primka se zdvihem h=',num2str(h)]);
      h,trand(x)
pause;
echo off
  end
echo on
% 3. a 4.test na bazi Kendallova a Spearmanova koeficientu reaguji
% nejdrive, o neco pozdeji pak 5. test (medianovy).
pause;
%---------------------------------------------------------------
% 5. Bily sum + primka se zdvihem h=xtr(n)=40*sigma az 50*sigma:
  figure;
  for h=40:50
      x = mi + (h./xtr(n))*xtr + wn;
      plot(t,x);
      title(['Bily sum WN(2,1) + primka se zdvihem h=',num2str(h)]);
      h,trand(x)
pause;
echo off
  end
echo on
% Test c.2 (body zvratu) zabira nejpozdeji.
%
pause;
close all;
%=================================================================
% VYZKOUSIME CITLIVOST NA SMES PERIDICKE KOMPONENTY
% A LINEARNIHO TRENDU
%=================================================================
% 6. Bily sum + primka se zdvihem h=xtr(n)=0.1*sigma az 1*sigma:
  figure;
  for h=0.1:0.1:1
      x = mi + (h./xtr(n))*xtr + h*xp + wn;
      plot(t,x);
      title(['Bily sum WN(2,1) + primka + sinusovka s A=h=',num2str(h)]);
      h,trand(x)
pause;
echo off
  end
echo on
% Posledni tri testy se opet jevi jako nejspolehlivejsi
pause;
%=================================================================
% NAKONEC OVERIME CITLIVOST NA KORELACE V DATECH
%=================================================================
% x=stacionarni MA(q) proces s korelacemi na vzdalenost q delky n
  n = 200;
  figure;
  for q=1:10
% zvolime nahodne realne koreny vne jednotkoveho kruhu, aby proces MA
% byl invertibilni
      r = 10*rand(1,q)-5; r = sign(r) + r
      theta=poly(r); theta=fliplr(theta);
      theta=theta./theta(1); theta(1)=[]
      x = armagen([],theta,n);
      plot(t,x);
      title(...
      ['Stacionarni MA(q) proces s korelacemi na vzdalenost q=',num2str(q)]);
      q,trand(x)
pause;
echo off
  end
echo on
% Na korelace v datech reaguje nejlepe test c.2 (body zvratu)
% a to vetsinou jiz pri korelacich na vzdalenost 1.
%================================== KONEC ==================================
echo off
