Please feel welcome to contribute still.
Before moving onto 2., let's discuss 1.
Let us have $N$ fermionic modes. We represent $\xi_i \, \hat{=} \, c_i$ with $i=1...N$, $c_i$ is a matrix given by the Jordan-Wigner mapping (see question), and $\xi^{*}_i \, \hat{=} \, c_i$ with $i=N+1...2N$. As pointed out by Chiral Anomaly, the representation is not faithful in that the representation of $\xi^{*}_i$ is not the hermitian conjugate of the representation of $\xi_i$. I believe it is at least briefly mentioned in various sources that $\xi^{*}_i$ need to be handled as modes independent of $\xi_i$, which is consistent with this. One can formally define rules for complex conjugation which maps $c_i$, $i\leq N$ to $c_j$, $j = N+i$. I don't know if a rule like this creates mayhem though. This also strangely means that you need $2^{2N} \times 2^{2 N}$ matrices, instead of $2^{N} \times 2^{N}$, as with real-valued Grassmann variables, and with quantum fermions! I don't think I understand the reason for this, but I just want to underline that I am not claiming that this is a useful way to perform the integrals for large-scale systems.
2. Creutz's sketched idea (see question) seems to work. Let the function of Grassmann variables be of the form $f(\boldsymbol{\xi},\boldsymbol{\xi^{*}}) = e^{-\sum_i \xi^{*}_i \xi_i} g(\boldsymbol{\xi},\boldsymbol{\xi^{*}})$, $\boldsymbol{\xi} = (\xi_1,...,\xi_N)$, where $e^{-\sum_i \xi^{*}_i \xi_i}$ is the normalization factor related to Bargmann coherent states (which you can consider to be part of the measure), and $g$ is a general function of Grassmann variables.
Using the above conventions, I believe the integral can be represented with
\begin{align}
I &= \int \prod_i \mathrm{d} \xi^{*}_i \xi_i \, e^{- \sum_j \xi^{*}_j \xi_j} g(\boldsymbol{\xi},\boldsymbol{\xi^{*}})
\\
\ &\hat{=}\
\langle 0 | (-1)^{N(N-1)/2} e^{- \sum_j c_{N+j} c_j} g(c_1,...,c_{2N}) \big( c^{\dagger}_{2N} ... c^{\dagger}_1 | 0 \rangle \big), \label{eq:main_result}
\end{align}
Note that $e^{- \sum_j c_{N+j} c_j} = \prod_j e^{- c_{N+j} c_j} = \prod_j (1 - c_{N+j} c_j)$, which looks like a projector. Everything on the rhs can now be plugged onto a computer.
We should probably try to motivate the result. Let us motivate the phase factor $(-1)^{N(N-1)/2}$, which changes when the number of even modes changes, since it was glanced over by Creutz. To this end, let's look at the constant function $g=1$, and evaluate $I_{1} = \langle 0|\prod_{k=1}^{N}(1 - c_{N+k} c_k) c^{\dagger}_{2N} ... c^{\dagger}_{1} |0\rangle$. We know that the non-zero contributions come from terms where the lhs operators annihilate each mode exactly once. Thus, we can simplify
$$
\langle 0 | ... (1 - c_{N+k} c_k) c^{\dagger}_{2N} ... c^{\dagger}_{1} |0\rangle
\to
-\langle 0 | ... (c_{N+k} c_k) c^{\dagger}_{2N} ... c^{\dagger}_{1} |0\rangle,
$$
since the term proportional to $1$ in $1-c_{N+k} c_k$ will vanish, because nowhere in the corresponding amplitude are modes $N+k$ and $k$ annihilated, leaving a mismatch compared with the occupied ket.
For $k=N$, we swap $c_k$ over $N$ lowering operators, and employ $c_k c^{\dagger}_k = -c^{\dagger}_k c_k + 1$. Only the latter term survives, as can be seen by swapping the $c_k$ all the way to right to annihilate the vacuum. We then perform the same tricks for $c_{N+k} = c_{2N}$, which is already on the lhs of $c^{\dagger}_{2N}$. Thus far we have
$$
I_{1} = \langle 0| (-1)(-1)^{N} \prod_{k=1}^{N-1}(1 - c_{N+k} c_k) c^{\dagger}_{2N-1} ... c^{\dagger}_{N+1} c^{\dagger}_{N-1} ... c^{\dagger}_{1} |0\rangle.
$$
Notice that for $k=N-1$, $c_k$ only has to be swapped $N-1$ times. Plowing through the whole product, in total we gain the phase factor $(-1)^{2N}(-1)^{\sum_{l=0}^{N-1} (N - l)} = (-1)^{N^2 - N(N+1)/2} = (-1)^{N(N-1)/2}$. This is not a proof for integrating over an arbitrary function, but I believe the result remains unchanged. You can first reduce the contributions from the function, which gives you some phase $(-1)^{a(i,N)}$, where $a(i,N)$ is a function. You then reduce the Bargmann norm terms $k=N,...,N-(i-1)$, and $k=i-1...1$ as before, which gives you phases that cancel the $i$-dependent part in $(-1)^{a(i,N)}$ and leave just $(-1)^{2N}(-1)^{N(N-1)/2}$ as before. To be fair, I don't have a general proof, just tested some cases to see how it works.
Note that I have numerically tested the above expression for $I$ from $N = 1...8$, up to $5$th order terms, and got the same results as with pencil and paper.
