Here's a different derivation. A derangement $D_n$ ($= D_{n,0}$) is a permutation of $n$ elements with no fixed points. We will prove an integral representation for $D_n$ that produces quick derivations of $D_n \approx n!/e$ and $D_n = n! \sum_{k=0}^n (-1)^k/k!.$
A combinatorial proof of the recurrence for $D_n$:
Here's a combinatorial proof of $D_n = (n-1)(D_{n-1} + D_{n-2})$ for $n \geq 2$ due to Euler.
For any derangement $(j_1, j_2, \ldots, j_n)$, we have $j_n \neq n$. Let $j_n = k$, where $k \in \{1, 2, \ldots, n-1\}$. We now break the derangements on $n$ elements into two cases.
Case 1: $j_k = n$ (so $k$ and $n$ map to each other). By removing elements $k$ and $n$ from the permutation we have a derangement on $n-2$ elements, and so, for fixed $k$, there are $D_{n-2}$ derangements in this case.
Case 2: $j_k \neq n$. Swap the values of $j_k$ and $j_n$, so that we have a new permutation with $j_k = k$ and $j_n \neq n$. By removing element $k$ we have a derangement on $n-1$ elements, and so, for fixed $k$, there are $D_{n-1}$ derangements in this case.
Thus, with $n-1$ choices for $k$, we have, for $n \geq 2$,
$$D_n = (n-1)(D_{n-1} + D_{n-2}). \tag{1}$$
An integral representation for $D_n$:
It turns out that $$D_n = \int_0^{\infty} (t-1)^n e^{-t} dt. \tag{2}$$
(The similarity of $(2)$ with the integral representation for the factorial as a gamma function is perhaps another reason it makes sense to call the derangements the "subfactorials.")
We will prove $(2)$ by showing that the integral satisfies $(1)$. Let $R_n = \int_0^{\infty} (t-1)^n e^{-t} dt$.
Applying integration by parts yields
$$R_n = (-1)^n + n R_{n-1} , \tag{3}$$
which is the recurrence $(8)$ in robjohn's answer. Applying the recurrence $(3)$ again with $R_{n-1}$ yields
$$
\begin{align}
R_n &= (-1)^n + (n-1)R_{n-1} + R_{n-1} \\
&= (-1)^n + (n-1)R_{n-1} + (-1)^{n-1} + (n-1)R_{n-2} \\
& = (n-1)(R_{n-1} + R_{n-2}).
\end{align}
$$
Since $R_0 = 1 = D_0$ and $R_1 = 0 = D_1$, Equation $(2)$ is established.
Equation $(2)$ is actually a special case of a more general result that says that the number of permutations with a specified set of fixed points can be represented by $\int_0^{\infty} R_{\tilde{G}}(t) e^{-t} dt$, where $\tilde{G}$ is the complement of $G$ in the complete bipartite graph on $n$ elements, and $R_G(t)$ is the associated rook polynomial for $G$. (See, for example, P. Mark Kayll, "Integrals Don't Have Anything to Do with Discrete Math, Do They?", Mathematics Magazine 84(2): 2011, 108-119.)
The approximation for $D_n$:
From $(2)$ we have
$$
\begin{align}
D_n &= \int_0^{\infty} (t-1)^n e^{-t} dt \\
&= \int_1^{\infty} (t-1)^n e^{-t} dt + \int_0^1 (t-1)^n e^{-t} dt \\
&= e^{-1} \int_0^{\infty} x^n e^{-x} dx + E_n \\
&= e^{-1} \Gamma(n+1) + E_n \\
&= \frac{n!}{e} + E_n.
\end{align}
$$
The quantity $E_n$ is small, too:
$$|E_n| < \left|\int_0^1 (t-1)^n dt \right| = \frac{1}{n+1},$$
so that $$D_n \approx \frac{n!}{e}.$$
The explicit formula for $D_n$:
Again from $(2)$ we have
$$
\begin{align}
D_n &= \int_0^{\infty} (t-1)^n e^{-t} dt \\
&= \int_0^{\infty} \sum_{k=0}^n \binom{n}{k} (-1)^k t^{n-k} e^{-t} dt \\
&= \sum_{k=0}^n \binom{n}{k} (-1)^k \int_0^{\infty} t^{n-k} e^{-t} dt \\
&= n! \sum_{k=0}^n \frac{(-1)^k}{k! (n-k)!} \Gamma(n-k+1) \\
&= n! \sum_{k=0}^n \frac{(-1)^k}{k!}.
\end{align}
$$
(I learned these arguments from the Kayll paper referenced above.)