crcv3
  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;
help powtr
        COMPUTING ESTIMATE OF PARAMETER lambda FOR THE POWER TRANSFORM:
                         y = x.^lambda   if lambda~=0
                           = log(x)      if lambda =0  (natural logarithm)
  =============================================================================
  function [lambdaE,y] = powtr(x,seglen,fig)
  INPUT:
  x      ... is an n x 1 or 1 x n data vector of POSITIVE observed data
             or vector of smoothed values (determnistic component) provided
             that the residual error component r describing the decomposition
             x = y + r is supplied instead of seglen
  seglen ... segment length used for computation of "std. deviation" growth
             (recommended: 4<=seglen<=12)
             or seglen=r is an n x 1 or 1 x n vector of the residual errors
             provided that the deterministic component y was supplied instead
             of x. In such a case the "std. deviation" growth will be evaluated
             more precisely from the dependence of r on y.
             optional: []=default segment length=8
  fig ... figure handle for the plot of estimated "std. deviation" versus mean
          optional: []=default=no plot
  OUTPUT:
  lambdaE ... =[lambda,E] is and 1 x 2 vector, where
               lambda is estimate of parameter lambda
               [lambda-E,lambda+E] is its confidence interval of 95%
  y  ... =data transformed according to parameter lambda
          vector of the same size as x

return
%==============================================================
% 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);
>>> CONSTANT VARIANCE (data not transformed): y=x <<<
  lambdaE
lambdaE =
   0.88897701949754   0.32871049790628
pause;
  [lambdaE,y] = powtr(x,10,figure);
>>> CONSTANT VARIANCE (data not transformed): y=x <<<
  lambdaE
lambdaE =
   1.05862009398911   0.17303456278115
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);
>>> POWER GROWTH OF VARIANCE (power transform): y=x.^1.4103 <<<
  lambdaE
lambdaE =
   1.41034862321665   0.28683644243842
pause;
  [lambdaE,y] = powtr(x,10,f+1);
>>> POWER GROWTH OF VARIANCE (power transform): y=x.^1.4047 <<<
  lambdaE
lambdaE =
   1.40472228654784   0.22612800643518
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);
>>> POWER GROWTH OF VARIANCE (power transform): y=x.^0.73024 <<<
  lambdaE
lambdaE =
   0.73023826579617   0.26195035389579
pause;
  [lambdaE,y] = powtr(x,10,f+1);
>>> POWER GROWTH OF VARIANCE (power transform): y=x.^0.7 <<<
  lambdaE
lambdaE =
   0.70000303963336   0.21220594451983
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);
>>> LINEAR GROWTH OF VARIANCE (logarithmic transform): y=log(x) <<<
  lambdaE
lambdaE =
   0.00842914620341   0.28577525324501
pause;
  [lambdaE,y] = powtr(x,10,f+1);
>>> LINEAR GROWTH OF VARIANCE (logarithmic transform): y=log(x) <<<
  lambdaE
lambdaE =
  -0.01452527226474   0.24824302489081
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);
>>> POWER GROWTH OF VARIANCE (power transform): y=x.^-0.571 <<<
  lambdaE
lambdaE =
  -0.57100296400894   0.30531977705222
pause;
  [lambdaE,y] = powtr(x,10,f+1);
>>> POWER GROWTH OF VARIANCE (power transform): y=x.^-0.42106 <<<
  lambdaE
lambdaE =
  -0.42106258782783   0.25225767415873
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);
>>> LINEAR GROWTH OF VARIANCE (logarithmic transform): y=log(x) <<<
  lambdaE
lambdaE =
  -0.23350784709744   0.42713487996569
pause;
  [lambdaE,y] = powtr(x,4,figure);
>>> LINEAR GROWTH OF VARIANCE (logarithmic transform): y=log(x) <<<
  xlabel(xinf(2,:));
  lambdaE
lambdaE =
  -0.10253855869531   0.21769530893707
pause;
%---------------------------------------------------------------
  figure; plot(y,'g');
  title(['Mocninna transformace dat s lambda=',num2str(lambdaE(1))]);
  xlabel(xinf(2,:));

%================================== KONEC ==================================
echo off
close all
help tperiod
         COMPUTING PERIODOGRAM AND DETECTING SIGNIFICANT PERIODICITIES
                       USING FISHER'S AND SIEGEL'S TEST
         Additive model: x(t)=y(t)+e(t) is assumed where
         y(t) ... deterministic component with periodicities to be detected
         e(t) ... error component, e: IID(0,sigma^2)
  =============================================================================
  function [I,kF,kS]=tperiod(x,fig)
  INPUT:
  x   ... is an n x 1 or 1 x n data vector of observations
  fig ... figure handle for the peridogram with emphasized significant
          harmonics
          optional: []=default=no plot
  OUTPUT:
  I ... periodogram: n2 x 1 or 1 x n2 vector sized according to x
                     n2=fix((n-1)/2)
  kF ... indices of significant harmonics using Fisher's test
  kS ... indices of significant harmonics using Siegel's test

diary off