Why non-linear response to random values is always positive?

I'm creating a non-linear response to a series of random values from {-1, +1} using a simple Volterra kernel:

With a zero mean for a(k) values I would expect r(k) to have a zero mean as well for arbitrary w values. However, I get r(k) with an always positive mean value, while a mean for a(k) behaves as expected: is close to zero and changes sign from run to run.

Why don't I get a similar behavior for r(k)? Is it because a(k) are pseudo-random and two different values from a are not actually independent?

Here is a code that I use:

import numpy as np
import matplotlib.pyplot as plt
import itertools

# array of random values {-1, 1}
A = np.random.randint(2, size=10000)
A = [x*2 - 1  for x in A]

# array of random weights
M = 3
w = np.random.rand(int(M*(M+1)/2))

# non-linear response to random values
R = []
for i in range(M, len(A)):
    vals = np.asarray([np.prod(x) for x in itertools.combinations_with_replacement(A[i-M:i], 2)])
    R.append(np.dot(vals, w))

print(np.mean(A), np.var(A))
print(np.mean(R), np.var(R))