EDIT
The method you want to use is ok, and gives a quick result. Here it is:
$$I=\int\prod d\theta^{*}d\theta\,\theta_{k}^{*}\theta_{l}exp(\theta^{*}B\theta+\eta^{*}\theta+\theta^{*}\eta)\\=\left(\frac{\partial}{\partial\eta_{k}^{*}}\right)\left(\frac{\partial}{\partial\eta_{l}}\right)\int\prod d\theta^{*}d\theta\, exp(\theta^{*}B\theta+\eta^{*}\theta+\theta^{*}\eta)$$
from where it follows that $I$ is equal to
$$\left.I=\mathrm{det}B\left(\frac{\partial}{\partial\eta_{k}^{*}}\right)\left(\frac{\partial}{\partial\eta_{l}}\right)\mathrm{exp}(\eta_{j}^{*}(B^{-1})_{ij}\eta_{i})\right\vert_{\eta^{*}=\eta=0}$$
$$I=\mathrm{det}B(B^{-1})_{kl}$$
Note: I use the shift $\theta_{i}\rightarrow\theta_{i} +(B^{-1})_{ij}\eta_{j}$ and $\theta_{i}^{*}\rightarrow\theta_{i}^{*} +\eta_{j}^{*}(B^{-1})_{ji}$. Here I used the fact that the second moments can be calculated as derivatives in the following way.
$$\left.\langle\eta_{l}\eta_{k}^{*}\rangle=\left(\frac{\partial}{\partial\eta_{k}^{*}}\right)\left(\frac{\partial}{\partial\eta_{l}}\right)\mathrm{exp}(\eta_{j}^{*}(B^{-1})_{ij}\eta_{i})\right\vert_{\eta^{*}=\eta=0}=(B^{-1})_{ij}$$
More or less, you can take this as a definition. But if you still want to prove this, do the following:
$$J=\left(\frac{\partial}{\partial\eta_{k}^{*}}\right)\left(\frac{\partial}{\partial\eta_{l}}\right)\mathrm{exp}(\eta_{j}^{*}(B^{-1})_{ij}\eta_{i})=\left(\frac{\partial}{\partial\eta_{k}^{*}}\right)\left(\frac{\partial}{\partial\eta_{l}}\right)\prod_{ij}(1-\eta_{j}^{*}(B^{-1})_{ij}\eta_{i})$$
After straight forward differentiation you get
$$J=\prod_{ij}\delta_{kj}\delta_{li}(B^{-1})_{ij}=(B^{-1})_{kl}$$
You can ignore whats below the line, that was my first answer. But I'll leave it because it may be instructive for others.
Since I cannot comment yet, i will sketch a proof for a simpler case
$$I=\int\prod_{i}d\theta_{i}^{*}d\theta_{i}exp(\theta_{i}^{*}B_{ij}\theta_{j})=det(B)$$
and hope that this will help you compute your integrals. After expanding the all the exponential we arrive at
$$I=\frac{1}{N!}\int d\theta_{1}^{*}d\theta_{1}\dots d\theta_{N}^{*}d\theta_{N}(\theta_{i_1}^{*}B_{i_{1}j_{1}}\theta_{j_{1}})(\theta_{i_2}^{*}B_{i_{2}j_{2}}\theta_{j_{2}})\dots (\theta_{i_N}^{*}B_{i_{N}j_{N}}\theta_{j_{N}})$$
At this point some explanations are in order. The factor $1/N!$ appears from the expansion of the exponentials. To see how the integrand was obtained, let's look at the case when $N=2$. We will have something like this
$$\int d\theta_{1}^{*}d\theta_{1}d\theta_{2}^{*}d\theta_{2}\left(1+\theta_{i}^{*}B_{ij}\theta_{j}+\frac{(\theta_{i}^{*}B_{ij}\theta_{j})^2}{2!}+\cdots\right)$$
It is obvious that only the quadratic term will contribute to the above integral, because only this term can saturate the number of grassmann variables in the integral measure.
$$\int d\theta_{1}^{*}d\theta_{1}d\theta_{2}^{*}d\theta_{2}\left(\frac{(\theta_{i}^{*}B_{ij}\theta_{j})^2}{2!}\right)=\frac{1}{2}\int d\theta_{1}^{*}d\theta_{1}d\theta_{2}^{*}d\theta_{2}\sum_{i_{1},i_{2},j_{1},j_{2}}^{2}(\theta_{i_1}^{*}B_{i_{1}j_{1}}\theta_{j_{1}})(\theta_{i_2}^{*}B_{i_{2}j_{2}}\theta_{j_{2}})$$
Expanding the sum and performing all four integrals, we get
$$\int d\theta_{1}^{*}d\theta_{1}d\theta_{2}^{*}d\theta_{2}\left(\frac{(\theta_{i}^{*}B_{ij}\theta_{j})^2}{2!}\right)=\frac{1}{2!}[2(B_{11}B_{22}-B_{12}B_{21})]=detB$$
Now, let us return to our original integral $I$. The next step before performing the integrals is to reorder the integrals and the grassmann numbers.
$$I=\frac{1}{N!}\int d\theta_{1}^{*}\dots d\theta_{N}^{*}\theta_{i_1}^{*}\dots\theta_{i_1}^{*}\int\theta_{1}\dots d\theta_{N}\theta_{j_1}\dots\theta_{j_1}B_{i_{1}j_1}\dots B_{i_{N}j_N}$$
From where we finally arrive at
$$I=\frac{1}{N!}\epsilon_{i_{1}\dots i_{N}}\epsilon_{j_{1}\dots j_{N}}B_{i_{1}j_1}\dots B_{i_{N}j_N}=\mathrm{det}B$$
(note: the ordering result $a_{1}b_{1}\dots a_{N}b_{N}=a_{1}\dots a_{N}b_{1}\dots b_{N}(-1)^{N(N-1)/2}$ has to be used twice for the integrals).
I found this method to be the most simple one of all when dealing with these sort of integrals. I hope this helps you prove those relations. The same method can be applied very easy in your case (just a strait forward extension). And with the method you described, it seams you are over complicating yourself. However, I'll work on it a bit and see what comes up.
Note that in OP's eq. (1) it is implicitly assumed that the matrix $M$ is antisymmetric
$$\tag{A} M^T~=~-M.$$
[A symmetric part in eq. (1) would not contribute to the integrand (1).] The term $\eta^TM^{-1}\eta$ in eq. (3) does in general not vanish
$$ 2S(\psi) ~:=~ (\psi - M^{-1}\eta)^T M(\psi - M^{-1}\eta)
~\stackrel{(A)}{=}~(\psi^T + \eta^T M^{-1}) M(\psi - M^{-1}\eta)$$
$$~=~\psi^T M\psi+\eta^T\psi -\psi^T\eta-\eta^TM^{-1}\eta . \tag{B}$$
One may check that if we vary the shifted action (B) wrt. the integration variable $\psi$, we unsurprisingly get the classical value
$$ \psi~\approx~M^{-1}\eta, \tag{C}$$
which reflect the shift (2) that OP performed in the first place.
As a check, note that the classical action vanishes
$$S\left(\psi=M^{-1}\eta\right)~\stackrel{(B)}{=}~0 ,\tag{D}$$
as it should.
Best Answer
The "real" integral evalautes to the Pfaffian of $A$ where for an $2n$-by$2n$ skew symmetric matrix $A$ $$ {\rm Pf}A= \frac 1{2^n n!} \epsilon_{i_1, \ldots, i_{2n}} A_{i_1,i_2}\ldots A_{2n-1,2n}. $$
The Pfaffian has the property that $({\rm Pf} A)^2= {\rm det}A$, and has applications many places in combinatorics. For Majorana fermion path integrals it replaces the one-loop Matthews-Salam determinant.
2: