Gaussian process with polynomial covariance function integrating to 0 In Rasmussen and Williams (2006, p. 88), one can find the following piecewise polynomial covariance function with compact support:
$$ k_\text{pp}(t,t') = (1-|t-t'|)_+ $$
where $(x)_+ = x \times [x>0] $ is the ramp function.
From this I construct a zero-mean windowed Gaussian process $$X(t) \sim \mathcal{GP}(0,k)$$ with the following covariance function:
$$ k(t,t') = \Pi(t) \ k_\text{pp}(t,t') \ \Pi(t') $$
where $\Pi(x) = [|x| < \frac{1}{2}]$ is the boxcar function. This ensures $X(t)$ is almost surely integrable as it is identically zero outside the boxcar window $[-\frac{1}{2},\frac{1}{2}]$.
Now, according to the accepted answer to Distribution over functions that integrate to 0, it is then possible to select a subspace of the $X(t)$ functions which integrate to zero by considering the joint distribution of $(X,Y)$, with $$Y = \int_{-\infty}^\infty X(t) \ dt \sim \mathcal{N}(0,C^2)$$ and $$C^2 = \int_{-\infty}^\infty \int_{-\infty}^\infty k(t,t') \ dt \ dt' = 2/3.$$
This is done by constructing a new Gaussian process $Z(t) = X(t) \ | \ Y = 0 \sim \mathcal{GP}(0,k_0)$. According to my understanding, its covariance function is then determined by the conditioning formula for Gaussians and given by
$$ k_0(t,t') = k(t,t') - \frac{c(t) \ c(t')}{C^2}, $$
where $$c(t) = \int_{-\infty}^\infty k(t,t') \ dt' = \frac{1}{4}(3-4t^2) \Pi(t).$$
Thus
$$ k_0(t,t') = \Pi(t) \{ (1-|t-t'|)_+ - \frac{3}{32} (3-4t^2)(3-4t'^2) \} \Pi(t'), $$
which indeed induces negative correlations at certain points $t,t'$, as the zero integral of $Z(t)$ implies that "what goes up must go down".
However, on sampling from $Z(t)$, I do not find that the functions integrate to zero. Can someone please point out the error in my reasoning?

I've done the math again and I have not found anything wrong with it. Here is a simple python snippet that shows that samples actually are zero mean...
import numpy as np

def k(t1, t2):
    ct1 = 0.75 - t1**2
    ct2 = 0.75 - t2.T**2
    return 1 - np.abs(t1 -t2.T) - 1.5 * ct1 * ct2

T = np.linspace(-0.5, 0.5, 101)[:, None]

K = k(T, T)

Z = np.random.multivariate_normal(np.zeros(101), K)
np.mean(Z)

