Standard unbiased estimator Kurtosis wrong implementation?

Question

So I implemented the standard unbaised estimator for the Kurtosis; https://en.wikipedia.org/wiki/Kurtosis#Standard_unbiased_estimator

However, it gives values which I do not expect.

I tested it for the Laplace, Wigner Cricle and Normal distributions. I would expect that the Laplace would give a positive Kurtosis, the normal distribution a Kurtosis value around zero and the Wigner Cricle a negative Kurtosis. However, all Kurtosis values are at $-3$..., so I checked my code but I cannot found any problems...

close all; clear all; clc;
y1  = randl(1,1000);
y2  = normrnd(0,1,1,1000);
y3  = betarnd(1.400,1.400,1,1000);

figure;
subplot(131);
% plot(x,yy1,'linewidth',2);
histogram(y1);
xlabel('Value [-]');
ylabel('Frequency [-]');
title(['Kurtosis \kappa > ',num2str(kur(y1))]);
axis square; box on; grid on;
subplot(132);
histogram(y2);
xlabel('Value [-]');
title(['Kurtosis \kappa = ',num2str(kur(y2))]);
axis square; box on; grid on;
subplot(133);
histogram(y3);
xlabel('Value [-]');
title(['Kurtosis \kappa = ',num2str(kur(y3))]);
axis square; box on; grid on;
xlim([-0.1 1.1]);
set(gcf,'Position', [100,100,800,400]);

[a1,a2,a3,a4] = kur(y1);
[b1,b2,b3,b4] = kur(y2);
[c1,c2,c3,c4] = kur(y3);

[a1,a2,a3,a4;
 b1,b2,b3,b4;
 c1,c2,c3,c4]

With

function x = randl(m, n)
u1 = rand(m, n);
u2 = rand(m, n);
x = log(u1./u2);
x = bsxfun(@minus, x, mean(x));
x = bsxfun(@rdivide, x, std(x));
end

and

function [y1,y2,y3,y4] = kur(x)
    n = length(x);

    k1 = (n + 1)*n/((n - 1)*(n - 2)*(n - 3));
    k2 = moment(x,4)/moment(x,2)^2;
    k3 = 3*(n - 1)^2/((n - 2)*(n - 3));

    y1 = 1/n*k2;        % Baised kurtosis
    y2 = 1/n*k2 - 3;    % Baised excess kurtosis
    y3 = k1*k2;         % Unbaised kurtosis
    y4 = k1*k2 - k3;    % Unbaised excess kurtosis
end