help acf
  DIRECT COMPUTATION OF THE AUTOCORRELATION (AUTOCOVARIANCE) FUNCTION ESTIMATE
                  (suitable for short sequences)
  CALLING:
  function [xro,xroE,mi,miE,xacf,stat] = acf(x,fig,alpha,K,P)
  INPUT:
  x ... time series data vector of length n
  fig ... figure handle for the acf plot (optional: default=[]=no plot)
          at most 50 acf values plotted.
  alpha ... a percentual confidence risk (optional: default=[]=5%)
  K     ... length of leading part of xro to be used for the computation of
            Portmanteau statistics (see stat on output)
            (optional: default=[]=min(round(sqrt(n)),length(xro)))
  P     ... x residuals => P=number of parameters in the associated model
        ... x is observed time series => P=0 (default=[]=0)
  OUTPUT:
  xro ... estimate of the autocorrelation function of length n/2 (column vector)
          (leading significant part only)
  xroE ... xroE(q+1)=asymptotic (1-alpha)% confidence halfwidth for xro(j+1), j>q
           (q=0,1,...n/2-1) under the hypothesis that x is an MA(q) process.
           (column vector)
  mi   ... estimate of the mean
  miE  ... = [miE1,miE2] = (1-alpha)% asymptotic confidence halfwidths of the mean
           miE2 is under assumption of normality of the generating noise
           miE2<miE1.
  xacf ... full estimate of the autocovariance function of length n
           (column vector)
  stat ... vector 3 x 1 allowing to accomplish Portmanteau test:
           stat(1) = value of the Portmanteau statistics:
                   = n*sum(xro(2:K+1).^2)                  for P=0
                   = n*(n+2)*sum(xro(2:K+1).^2./(n-(1:K))) for P>0 (residuals)
           stat(2) = chi2inv(1-0.01*alpha,K-P) ... chi^2 quantile
                   = NaN if K<=P (too many parameters in the model)
           stat(3) = K-P ... degrees of freedom
           Under the hypothesis xro(h)=0 for h>1 (x=white noise)
           stat(1) is approximately distributed as chi^2(K-p) and thus
           stat(1) > stat(2) => we reject with risk alpha % the hypothesis that
                                x could be a white noise.

% Bily sum
n=1000; mi=10; sigma=3;
wn=mi+sigma*randn(n,1);  % ~WN(mi,sigma^2)
colordef none   % cerne pozadi obrazku
plotdata(wn,'Bily sum')
plotdata(acf(wn),'Autokorelacni funkce bileho sumu')
acf(wn,figure);
WARNING: The process behaves like non-stationary => the estimate miE may not be correct
mean(wn)
ans =
    9.8708
std(wn)
ans =
    2.8305

% Trend v casove rade zpomaluje pokles acf
t=(1:n).';
k=0; x=k*t+randn(n,1); acf(x,figure);
WARNING: The process behaves like non-stationary => the estimate miE may not be correct
k=0.001; x=k*t+randn(n,1); acf(x,figure);
k=0.002; x=k*t+randn(n,1); acf(x,figure);
k=0.003; x=k*t+randn(n,1); acf(x,figure);

% Periodicita v casove rade
x=sin(100*pi*t./n)+randn(n,1); acf(x,figure);
close all

% acf pro radu z vety 2.33
sigma=2; z=sigma*randn(n+1,1);
r=0.5; theta1=(1+sqrt(1-4*r^2))/(2*r);theta2=(1-sqrt(1-4*r^2))/(2*r);
x1=z(2:n+1)+theta1*z(1:n); x2=z(2:n+1)+theta2*z(1:n);
acf(x1,figure);
acf(x2,figure);
r=-0.5; theta1=(1+sqrt(1-4*r^2))/(2*r);theta2=(1-sqrt(1-4*r^2))/(2*r);
x1=z(2:n+1)+theta1*z(1:n); x2=z(2:n+1)+theta2*z(1:n);
acf(x1,figure);
acf(x2,figure);
% atd. pro ruzna -0.5<=r<=0.5
diary off