function [m, m_var, m_lb, m_ub] = calc_subX_stats(fname, level, method, diagnosticPlot, opt1, opt2)

% mean, variance and CI for number of subcrossings of crossing of size 2^level
%
% fname: base file name
% level: >= 1, level of crossings used
% method: 'iid' for iid Z_i, 'srd' for correlated short-range dept Z_i,
%     'lrd for long-range dept Z_i
% diagnosticPlot: if 1 then diagnostic plots are produced (for methods srd or lrd only)
% opt1, opt2: option parameters
%     'iid' uses no options
%     'srd' uses opt1 in (0,1) for proportion of data used to calc delta, default 0.5
%     'lrd' uses opt1 and opt2, both passed to calc_alpha_exponent, defaults 10 and 1
%     opt1 is order of the wavelet basis, opt2 is the onset ot the scaling region
%
% requires function calc_alpha_exponent
%
% Y.Shen 12.2002, O.D.Jones 7.2003

if level < 10
    fnameXX = sprintf('%s_0%d', fname, level);
else
    fnameXX = sprintf('%s_%d', fname, level);
end
eval(['load ' fnameXX '.sub;']);
eval(['Z = ' fnameXX ';']);

if method == 'iid'
    
    m = mean(Z);
    S2 = var(Z);
    n = length(Z);
    m_var = S2/n;
    err = 1.96*sqrt(S2/n);
    m_lb = m - err;
    m_ub = m + err;
    
elseif method == 'srd'
    
    if nargin <= 4
        opt1 = 0.5;
    end
    [rho, m, S2] = autocorr(Z);
    n = length(Z);
    for N = 2:n*opt1
        x = 2*(1 - [1:N-1]/N).*rho(2:N);
        delta(N) = sum(x);
    end
    m_var = S2*(1 + delta(N))/n;
    err = 1.96*sqrt(S2*(1 + delta(N))/n);
    m_lb = m - err;
    m_ub = m + err;
    
    % diagnostic plot, if series is srd then delta(N) converges
    if diagnosticPlot == 1
        clf
        subplot(2,1,1);
        plot(1 + delta);
        title('variance correction factor')
        subplot(2,1,2);
        plot(rho);
        title('autocorrelation')
    end
    
elseif method == 'lrd'
    
    if nargin <= 4
        opt1 = 10;
    end
    if nargin <= 5
        opt2 = 1;
    end
    m = mean(Z);
    n = length(Z);
    
    [alpha, cfc, var_a, var_b, cov_ab] = calc_alpha_exponent(Z, diagnosticPlot, opt1, opt2);
    % a = log2(cfc);
    b  = 1 - alpha;
    if b > 1
        disp('WARNING: sequence of subcrossings apparently not LRD');
    end
    m_var = cfc*(2*pi*n)^(-b)*f(b)^2;
    [lb, ub] = ZW_CI(n, b, var_a, var_b, cov_ab, 0.95);
    m_lb = m - sqrt(m_var)*ub;
    m_ub = m - sqrt(m_var)*lb;
end


function [rho, m, S2] = autocorr(x)

% empirical autocorrelation function
% rho(r+1) is lag r autocorrelation of x
% m is mean of x
% S2 is variance of x

m = mean(x);
S2 = var(x);
x = x - m;

n = length(x);
for r = 0:n-1
    rho(r+1) = sum(x(1:n-r).*x(r+1:n))/(n-r)/S2;
end


function [lb, ub] = ZW_CI(n, b, var_a, var_b, cov_ab, sig_lvl);

% Monte-Carlo estimate of CI for Z/W

sample_size = 10000; % size of Monte-Carlo sample
delta_b = 0.001; % used to estimate f'
Z = randn(1,sample_size);
C = [var_a cov_ab;cov_ab var_b];
R = chol(C); 
L = R';
Z2 = L*randn(2,sample_size);
W = 2.^(Z2(1,:)/2).*(2*pi*n).^(Z2(2,:)/2).*(1 + Z2(2,:)*(f(b+delta_b) - f(b))/delta_b/f(b));
ZW = sort(Z./W);
lb = ZW(round((1-sig_lvl)/2*sample_size));
ub = ZW(round((1 - (1-sig_lvl)/2)*sample_size));


function y = f(x)

% used by ZW_CI

y2 = 2*pi*x/((1-x)*(2-x)*gamma(1-x)*sin(pi*x/2)*(2-2^(1-x)));
y = sqrt(y2);