%==============================================================
%    DETEKCE HETEROSKEDASTICITY (NEKONSTANTNOSTI ROZPTYLU)
%==============================================================

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

% Simulacni studie pro dekompozicni model:
% x(t) = d(t) + e(t),  t=1:n   ... !! MUSI BYT x(t)>0 !! - jinak posuneme
% d(t) = 2 + sin(2*pi*t/n)     ... sinusovy trend (1 perioda)
% e(t) ~ N(0,sigma*d(t)^theta) ... gaussovsky nekorelovany sum
%        s exponencialni zavislosti smerodatne odchylky
%        (a tedy i rozptylu) na stredni hodnote (trendu) d(t)
% Konstantniho rozptylu dosahneme mocninnou transformaci:
% y = x.^lambda   if lambda~=0
%   = log(x)      if lambda =0  (natural logarithm)
% kde lambda = 1-theta
% Pro odhad lambda uzijeme procedury 'powtr'.
  echo on
  format compact
  format long
%---------------------------------------------------------------
% Nastudujte nyni pomoci 'help powtr' pouzivani teto procedury,
% pripadne provedte jeji uplny vypis pomoci 'type powtr' a prostudujte
% podrobneji cely algoritmus.
keyboard;
%==============================================================
% GENEROVANI SIMULACNICH DAT
  n = 400;
  t = (1:n).';
  d = 2 + sin(2*pi*t./n);  % Sinusovy trend (1 perioda)
  sigma = 0.2;
pause;
%===============================================================
% 1) SIMULACE KONSTANTNIHO ROZPTYLU (theta=0 <=> lambda=1)
  theta = 0;
  x = d + sigma*d.^theta.*randn(n,1);
  figure; plot(t,x,'w'); hold on; plot(t,d,'r','LineWidth',2); hold off;
  f=gcf;
  title('SIMULACE KONSTANTNIHO ROZPTYLU');
pause;
%---------------------------------------------------------------
% Odhadneme lambda pomoci 'powtr':
% Pro porovnani nejprve presnejsi metodou vyuzivajici idealni rezidua x-d
  [lambdaE,y] = powtr(d,x-d,figure);
  lambdaE
pause;
  [lambdaE,y] = powtr(x,10,figure);
  lambdaE
pause;
%---------------------------------------------------------------
  figure; plot(t,x,'w',t,0.7+x.^lambdaE(1),'g');
  title(['Mocninna transformace dat s lambda=',num2str(lambdaE(1))]);
pause;
%===============================================================
% 2) SIMULACE POMALU (DLE PREVRACENE ODMOCNINY) KLESAJICIHO ROZPTYLU
%    (theta=-0.5 <=> lambda=1.5)
  theta = -0.5;
  x = d + sigma*d.^theta.*randn(n,1);
  figure(f); plot(t,x,'w'); hold on; plot(t,d,'r','LineWidth',2); hold off;
  title('SIMULACE POMALU (DLE 1/sqrt) KLESAJICIHO ROZPTYLU');
pause;
%---------------------------------------------------------------
% Odhadneme lambda pomoci 'powtr':
% Pro porovnani nejprve presnejsi metodou vyuzivajici idealni rezidua x-d
  [lambdaE,y] = powtr(d,x-d,f+1);
  lambdaE
pause;
  [lambdaE,y] = powtr(x,10,f+1);
  lambdaE
pause;
%---------------------------------------------------------------
  figure(f+2); plot(t,y,'g');
  title(['Mocninna transformace dat s lambda=',num2str(lambdaE(1))]);
pause;
%===============================================================
% 3) SIMULACE POMALU (DLE ODMOCNINY) ROSTOUCIHO ROZPTYLU
%    (theta=0.5 <=> lambda=0.5)
  theta = 0.5;
  x = d + sigma*d.^theta.*randn(n,1);
  figure(f); plot(t,x,'w'); hold on; plot(t,d,'r','LineWidth',2); hold off;
  title('SIMULACE POMALU (DLE ODMOCNINY) ROSTOUCIHO ROZPTYLU');
pause;
%---------------------------------------------------------------
% Odhadneme lambda pomoci 'powtr':
% Pro porovnani nejprve presnejsi metodou vyuzivajici idealni rezidua x-d
  [lambdaE,y] = powtr(d,x-d,f+1);
  lambdaE
pause;
  [lambdaE,y] = powtr(x,10,f+1);
  lambdaE
pause;
%---------------------------------------------------------------
  figure(f+2); plot(t,y,'g');
  title(['Mocninna transformace dat s lambda=',num2str(lambdaE(1))]);
pause;
%===============================================================
% 4) SIMULACE LINEARNE ROSTOUCIHO ROZPTYLU (theta=1 <=> lambda=0)
  theta = 1;
  x = d + sigma*d.^theta.*randn(n,1);
  figure(f); plot(t,x,'w'); hold on; plot(t,d,'r','LineWidth',2); hold off;
  title('SIMULACE LINEARNE ROSTOUCIHO ROZPTYLU');
pause;
%---------------------------------------------------------------
% Odhadneme lambda pomoci 'powtr':
% Pro porovnani nejprve presnejsi metodou vyuzivajici idealni rezidua x-d
  [lambdaE,y] = powtr(d,x-d,f+1);
  lambdaE
pause;
  [lambdaE,y] = powtr(x,10,f+1);
  lambdaE
pause;
%---------------------------------------------------------------
  figure(f+2); plot(t,y,'g');
  title(['Logaritmicka transformace dat pro lambda=',num2str(lambdaE(1))]);
pause;
%===============================================================
% 5) SIMULACE EXPONENCIALNE ROSTOUCIHO ROZPTYLU
%    (theta=1.5 <=> lambda=-0.5)
  theta = 1.5;
  x = d + sigma*d.^theta.*randn(n,1);
  figure(f); plot(t,x,'w'); hold on; plot(t,d,'r','LineWidth',2); hold off;
  title('SIMULACE EXPONENCIALNE ROSTOUCIHO ROZPTYLU');
pause;
%---------------------------------------------------------------
% Odhadneme lambda pomoci 'powtr':
% Pro porovnani nejprve presnejsi metodou vyuzivajici idealni rezidua x-d
  [lambdaE,y] = powtr(d,x-d,f+1);
  lambdaE
pause;
  [lambdaE,y] = powtr(x,10,f+1);
  lambdaE
pause;
%---------------------------------------------------------------
  figure(f+2); plot(t,y,'g');
  title(['Mocninna transformace dat s lambda=',num2str(lambdaE(1))]);
pause;
%===============================================================
% 6) VYZKOUSIME NYNI POSTUP NA REALNYCH DATECH:
% 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);
  t=(1:n).';
  pause;
%---------------------------------------------------------------
% Odhad parametru lambda pro mocninnou transformaci pomoci 'powtr':
% Pro porovnani nejprve presnejsi metodou vyuzivajici rezidua x-d,
% kde d nalezneme prolozenim paraboly zabudovanou funkci 'polyfit'.
  p = polyfit(t,x,2);
  d = polyval(p,t);
  figure(f); plot(t,x,'w'); hold on; plot(t,d,'r','LineWidth',2); hold off;
  title(xinf(1,:));
  xlabel(xinf(2,:));
pause;
  [lambdaE,y] = powtr(d,x-d,f+1);
  lambdaE
pause;
  [lambdaE,y] = powtr(x,4,figure);
  xlabel(xinf(2,:));
  lambdaE
pause;
%---------------------------------------------------------------
  figure; plot(y,'g');
  title(['Mocninna transformace dat s lambda=',num2str(lambdaE(1))]);
  xlabel(xinf(2,:));

%================================== KONEC ==================================
echo off
