

function [acf, lags, nse, rne, ineffactor] = ConvDiagnostics(x,nlags)

% Output
%--------------------------------------------------------------------------
% acf        : autocorrelation function
% lags       : lags for autocorrelation function
% nse        : numerical standard error using batch means method
% rne        : relative numerical efficiency, calculated using batch means
% method. Small values of rne denotes substantial autocorrelation
% ineffactor : reciprocal of rne. Large values denote the cost of working 
% with correlated samples 
% 
% Input
% x          : column vector
% nlags      : specify the number of lags, optional.
%--------------------------------------------------------------------------

[acf,lags] = autocorr(x,NumLags=400);

if nargin < 2 
   % nlags = 20;
   acfmin = min(abs(acf));
   if acfmin < 0.05
      nlags = find(abs(acf)<0.05,1,'first');
      disp('acf < 0.05')
   elseif acfmin >= 0.05 && acfmin < 0.15
      nlags = find(abs(acf)<0.15,1, 'first');
      disp('0.05 < acf < 0.15')
   elseif acfmin >= 0.15 &&  acfmin <0.30
      nlags = find(abs(acf)<0.30,1, 'first');
      disp('0.15 < acf < 0.30')
   else    
   nlags = size(x,1)/2;            % want nlags > 0, will make two batch
   disp('minimum acf > 0.30')
   end
end
 

%--------------------------------------------------------------------------
% Calculations for numerical standard error and relative numerical
% efficiency
n         = size(x,1);
nbatch    = floor(n/nlags);     % number of batches should be an integer
nuse      = nbatch*nlags;       % number of obs in series x, that we use 

% check
if rem(nuse,nlags) ~=0
   error('number of batches must be an integer')
end

x            = x(1:nuse,1);                % redefining x
xbatch       = reshape(x,nlags,nbatch);
mxbatch      = mean(xbatch);             % mean of xbatch

varxbatch    = sum((mxbatch' - mean(x)).^2)/(nbatch-1);
nse          = sqrt(varxbatch/(nbatch));

rne          = (std(x,1)/sqrt(nuse))/nse;
ineffactor   = 1/rne;

disp('       ntotal      nuse        nbatch      nlags')
disp([n      nuse    nbatch     nlags])

% n          : Length of the series x.
% nuse       : equals nbatch*nlags
% nbatch     : number of batches
% nlags      : Lag at which autocorrelation<0.05. Is only displayed not
% kept in output.

end


























