Wikipedia presents a recursive definition of the Pfaffian of a skewsymmetric matrix as $$ \operatorname{pf}(A)=\sum_{{j=1}\atop{j\neq i}}^{2n}(1)^{i+j+1+\theta(ij)}a_{ij}\operatorname{pf}(A_{\hat{\imath}\hat{\jmath}}),$$ where $\theta$ is the Heaviside step function. I have failed to find a proper reference for this result. Is it standard or does it require proper citation? What should one cite?

2$\begingroup$ Maybe a particular case of Theorem 1 in dx.doi.org/10.1006/aima.1995.1029 (no time to check, sorry). Should also be in Knuth's Overlapping Pfaffians paper. $\endgroup$ – darij grinberg Apr 12 '17 at 13:47

$\begingroup$ For this kind of identities or more complicated ones GrassmannBerezin integration is a powerful tool. See for example sciencedirect.com/science/article/pii/S0196885803001465 $\endgroup$ – Abdelmalek Abdesselam Apr 12 '17 at 13:57

$\begingroup$ @darijgrinberg: Yep, a prticular case of Theorem 1 of A Simple Proof of an Identity Concerning Pfaffians of Skew Symmetric Matrices by Dress, Wenzel. Thanks. $\endgroup$ – Vít Tuček Apr 12 '17 at 14:15

$\begingroup$ @AbdelmalekAbdesselam: Thanks. Could be useful some day. $\endgroup$ – Vít Tuček Apr 12 '17 at 14:16
Here is how you get the Pfaffian. $\newcommand{\bR}{\mathbb{R}}$ Consider the symplectic $2$form
$$ \omega = dx^1\wedge dx^2+\cdots +x^{2n1}\wedge x^{2n} \in\Omega^2(\bR^{2n}). $$
Note that
$$ \frac{1}{n!}\omega^{\wedge n}:=\frac{1}{n!}\underbrace{\omega\wedge \cdots \wedge \omega}_n=dx^1\wedge dx^2\wedge \cdots \wedge dx^{2n1}\wedge dx^{2n} $$
Suppose that $A=(a_{ij})$ is a skewsymmetric $2n\times 2n$ matrix. None that $\newcommand{\be}{\boldsymbol{e}}$
$$a_{ij}=(\be_i, A\be_j), $$
where $(\be_i)$ is the canonical basis of $\bR^{2n}$. Consider now the $2$form
$$ \omega_A=\sum_{i<j} aa_{ij}dx^i\wedge dx^j=\sum_{i<j} (A\be_i,\be_j) dx^i\wedge dx^j=\frac{1}{2}\sum_{i\neq j} a_{ij} dx^i\wedge dx^j. $$
The top exterior power of $\omega_A$ is a multiple $\mu$ of $\frac{1}{n!}\omega^{\wedge n}$
$$\omega_A^{\wedge n}=\frac{\mu}{n!} \omega^{\wedge n}. \tag{*}$$
The Pfaffian of $A$ is then the multiple $\mu$. If you write $\omega_A$ in the form
$$\omega_A =\frac{1}{2}\sum_j a_{ij} dx^i\wedge dx^j\frac{1}{2}\sum_{k,\ell\neq i} a_{k\ell} dx^k\wedge dx^\ell $$ and use the equality (*) then you will get your formula.

$\begingroup$ Thank you! This is pretty much the approach used in the paper that was linked by Darij Grinberg in comments where it is used to prove more general identity. $\endgroup$ – Vít Tuček Apr 12 '17 at 14:17

1$\begingroup$ This is essentially Berezin or fermionic integration. $\endgroup$ – Liviu Nicolaescu Apr 13 '17 at 1:10
Your formula can be found in
Fulton, William, Pragacz, Piotr: Schubert varieties and degeneracy loci, Lecture Notes in Mathematics 1689, Springer. xi, 148 p. (1998). ZBL0913.14016,
see in particular equation (D.1) p. 116.