Estimating time-lagged mutual information for two signal samples

Question

This is an attempt to reproduce Moon et al. 1995, and the author's copy can be obtained through here.

As a benchmark, we estimate the time-lagged mutual information of a simple sine signal $\sin(0.02\pi t)$ using Gaussian kernels: $$K(y)=\dfrac{\exp\left(-\dfrac{(y-y_i)^TS^{-1}(y-y_i)}{2h^2}\right)}{(2\pi)^{d/2}h^d\det(S)^{1/2}},$$ where $h$ is chosen to be $$h=\left\{\dfrac{4}{(d+2)}\right\}^{1/d+4}n^{-1/(d+4)},$$ where $d$ is the dimension of the multivariate kernel, and $n$ is the number of data points.

I wrote a Matlab function (the code is at the end of the post).

Suppose $x_t$ is the original time series $\sin(0.02\pi t)$ with 400 data points $t=1,2,\cdots,400$, and $x_{t-\tau}$ is $\sin[0.02\pi (t-\tau)]$. The author plotted the estimated mutual information $I_{x_t,x_{t-\tau}}$ as a function of $\tau$:

However, what I got from my code is

While the qualitative features (like the diverging mutual information at lag 0, 50 and 100), but the magnitude is far off and the overall shape is not right.

My 1st suspicion is the fact that I used different window widths for joint pdf and marginal pdf, and they may not be consistent. So my question is: how to consistently construct joint and marginal pdf estimate in the context of mutual information estimation?

Second question is: is there anything wrong in the code besides the possibly inconsistent pdf KDEs?

The Matlab code:

function [MI] = get_MI(xt, xt_lag)
% xt is original time series, xt_lag is the lagged one, both are column
% vectors
x_pair = [xt' xt_lag'];
n=length(xt);

d=1;
h_d1=(4/(d+2))^(1/(d+4)) * n^(-1/(d+4));

d=2;
h_d2=(4/(d+2))^(1/(d+4)) * n^(-1/(d+4));

MI = 0.;

Sxy_pair = cov(x_pair); invS_pair = inv(Sxy_pair); detS_pair = det(Sxy_pair);
Sxy_xt = cov(xt'); invS_xt = inv(Sxy_xt); detS_xt = det(Sxy_xt);
Sxy_lag = cov(xt_lag'); invS_lag = inv(Sxy_lag); detS_lag = det(Sxy_lag);

for i=1:n
    Psq = p_mkde(x_pair(i,:)', x_pair', h_d2, invS_pair, detS_pair);
    Ps = p_mkde(x_pair(i,1), xt, h_d1, invS_xt, detS_xt);
    Pq = p_mkde(x_pair(i,2), xt_lag, h_d1, invS_lag, detS_lag);

    MI = MI + (log2(Psq) - log2(Ps) - log2(Pq));
end

MI = MI/n;

end

function [pxy]=p_mkde(x,X,h,invS,detS)

s1=size(X); d=s1(1); N=s1(2);

sum=0;
for ix=1:N
    p2=(x-X(:,ix))'*invS*(x-X(:,ix));
    sum=sum+1/sqrt((2*pi)^d*detS)*exp(-p2/(2*h^2));
end
pxy=1/(N*h^d)*sum;
end