function [kim_t, dkim_t] = compute_implied_rate_with_error(x_t, x_ddag, dt, tmax)
% Compute implied rate estimate (and standard error) as a function of observation time.
%
% [kim_t, dkim_t] = compute_implied_rate_with_error(x_t, x_ddag, tmax)
%
% ARGUMENTS
%   x_t - timeseries
%   x_ddag - dividing surface separating A and B states
%   dt - time interval between observations in x_t
%   tmax - maximum lag time to compute implied rate
%
% RETURNS
%   kim_t - implied rate as a function of lag time, up to time tmax
%   dkim_t - statistical error in implied rate; dkim_t(t) is the standard error in kim_t(t)

% Determine length of timeseries.
T = length(x_t);

% Quantize using dividing surface to define h_A and h_B functions.
h_A = (x_t <= x_ddag);
h_B = (x_t > x_ddag);

% Compute stationary probabilities.
PA = mean(h_A);
PB = mean(h_B);

% Compute second eigenvalue of 2x2 transition matrix.
mu2_t = zeros(1,tmax);
dh_A = h_A - PA;
denom = mean(dh_A.*dh_A);
for t = 1:tmax
  mu2_t(t) = mean(dh_A(1:T-t) .* dh_A(1+t:T)) / denom;
end
mu2_t = [1 mu2_t];

% Compute implied rate constant.
tval = [1:tmax] * dt;
kim_t = - tval.^(-1) .* log(mu2_t(2:end));

% Compute number of effectively uncorrelated snapshots from maximum estimate of relaxation time constant.
tau = max(kim_t.^(-1));
g = 1 + 2*tau;
Neff = T / g;

% Compute variance in estimate of eigenvalue of 2x2 transition matrix.
d2mu2_t = zeros(1,tmax);
for t = 1:tmax
  F_n = h_A(1:T-t);
  G_n = h_A(1:T-t) .* h_B((1+t):T);
  F = mean(F_n);
  G = mean(G_n);
  d2F = var(F_n) / Neff;
  d2G = var(G_n) / Neff;
  dFdG = cov(F_n, G_n) / Neff; dFdG = dFdG(1,2);
  dmu_dF = G*(1-2*F)/(F^2 * (1-F)^2);
  dmu_dG = 1.0 / (F*(1-F));
  d2mu2_t(t) = dmu_dF.^2 .* d2F + dmu_dG.^2 .* d2G + 2*dmu_dF*dmu_dG*dFdG;
end
d2mu2_t = [0 d2mu2_t];

% Compute standard error estimate.
dkim_t = sqrt(tval.^(-1) .* d2mu2_t(2:end) .* mu2_t(2:end).^(-2));

return
