error caused by Wien-Hammerstein system (Bussgang?)

Hi,
I have the following setup:
x --->[ G(z) ]----->[ (.)^3 ]----->[ H(z) ]-----> y
where x is a Gaussian input with E(x)=0 and autocorrelation Rxx.
For simplicity, I can set Rxx=sigma^2*delta(tau) and G(z)=1, i.e. a pure
Hammerstein system driven with white Gaussian noise (I would like to
generalize afterwards if possible):
x --->[ (.)^3 ]----->[ H(z) ]-----> y
Now I would like to obtain the error in terms of its normalized mean
squared error caused by this system: error = E((y-x)^2)/E(x^2)
When H(z)=1, this is trivial: E((x^3)^2)/E(x^2) = 15*sigma^4.
With H(z) != 1, I think maybe Bussgang's theorem can help but I do not
see how. As it is written in Wikipedia, it works for Wiener systems but
not Hammerstein.
But I found in [1]: "For a Gaussian input signal, the iterative approach
can be motivated by Bussgang's theorem. This result implies that for a
Hammerstein system, the LTI model that minimizes the mean-square error
E((y(t) − G(q)u(t))^2)
will be equal to a scaled version of the LTI part of the system." This
implies to me that the best fit is obtained by C*H(z) where C is
obtained using the formula from Wikipedia [2]: C = 3*sigma^2.
However, a quick MATLAB check shows that something is wrong:
N = 2^16;
sigma = 1;
C = 3*sigma^2;
b = fir1(10, 0.1);
x = sigma*randn(N,1);
y1 = filter(b, 1, x.^3);
y2 = C*filter(b, 1, x);
plot([y1 , y2])
y1 and y2 are way too different, std(y1) and std(y2) completely different.
Is there any mistake. Or is there a better approach in the first place
to solve this problem?
Thanks,
Peter
[1] Martin Enquist, "Identification of Hammerstein Systems Using
Separable Random Multisines"
[2] https://en.wikipedia.org/wiki/Bussgang_theorem