show this function $f(m,n)=f(n,m)$.

First, let us find a simple bound: We have $|f(n,m)| \le m! n!$

In fact, for $m=0$ or $n=0$ this is true. If then $m=1$, we have $$f(1,n) = f(0,n) - n f(0,n-1) = 1- n.$$ Thus, we can suppose that $m \ge 2$. Then by induction, we get that \begin{align} |f(m,n)| &\le |f(m-1,n)| + n |f(m-1,n-1)| \\ & \le (m-1)! n! + n (m-1)! (n-1)! = 2[n! (m-1)!] \le n! m!. \end{align}

Next we determine the exponential generating function. (We will see that this function is symmetric.)

Note that we don't know if the generating function is convergent in an appropriate region. Therefore, we need to consider the exponential generating function.

Define $$G(z,w) = \sum_{m,n=0}^\infty \frac{f(m,n)}{m!n!} z^m w^n$$ and note that $G(z,w)$ is analytic on $|z|,|w| < 1$. Using the recurrence relation, we find \begin{align*} G(z,w) &= (e^z +e^w-1) + \sum_{m,n=1}^\infty \frac{f(m,n)}{m!n!} z^m w^n \\ &= (e^z +e^w-1) + \sum_{m,n=1}^\infty \Big(\frac{f(m-1,n)}{m!n!} - \frac{f(m-1,n-1)}{m! (n-1)!} \Big) z^m w^n \\ &= e^w + \sum_{m,n=0}^\infty\frac{f(m,n)}{(m+1)!n!} z^{m+1} w^n - \sum_{m,n=0}^\infty \frac{f(m,n)}{(m+1)! n!} z^{m+1} w^{n+1} \\ &= e^w+ \int_0^z G(t,w) \, dt - \int_0^z w G(t,w) \, dt \end{align*} This leads to the ordinary differential equation $$G_z(z,w) = (1-w) G(z,w).$$ Thus $G(z,w) = C(w) \exp ((1-w) z)$. Since $e^w = G(0,w) = C(w)$ we find that $$G(z,w) = \exp(z+w-zw).$$ The last identity implies that $G(z,w) =G(w,z)$ and thus the taylor coefficients must be symmetric, too.

Moreover, we can use this result in order determine an explicit formula for $f(n,m)$.

In fact, we have \begin{align*} \exp(z+w-zw) &= \sum_{n=0}^\infty \frac{1}{n!} (z+w-zw)^n \\ & = \sum_{m,n=0}^\infty z^m w^n \sum_{n= k_1+k_3, m= k_2+k_3} (-1)^{k_3} \frac{1}{k_1! k_2! k_3!} \end{align*}

and therefore $$f(n,m) = n! m! \sum_{n= k_1+k_3, m= k_2+k_3} (-1)^{k_3} \frac{1}{k_1! k_2! k_3!}.$$

Just an idea, first it's clear that, for $m\in\mathbb{N}^{+}$ $$f(m,0)=f(m-1,0)=f(m-2,0)=\ldots=f(m-(m-1),0)=f(1,0)=1$$ Therefore, $$f(m,1)=f(1,1)-f(1,0)-f(2,0)-f(3,0)-\ldots-f(m-1,0)=f(1,1)-(m-1)$$ $$=1-f(0,0)-(m-1)=2-f(0,0)-m\qquad (*)$$

For the other hand, $$f(1,m)=f(0,m)-mf(0,m-1)=1-m$$ From where it follows that, $f(1,1)=0$. Therefore, by $(*)$ we have $0=f(1,1)=1-f(0,0)$, i.e., $f(0,0)=1$. And from $(*)$ we have, $f(m,1)=1-m$. So, in this case, we have $f(1,m)=f(m,1)$.

Another realtion, that maybe help us is that $$f(2,m+1)=f(1,m+1)-(m+1)f(1,m)=m^2-m-1$$

From this, $f(2,2)=-1$ and for other hand $f(2,2)=f(1,2)-2f(1,1)=(1-2)-2*0=-1$. So this makes me think that the function is symmetrical.

With the MATHEMATICA script

Clear[f]
f[m_, n_] := If[(n > 0) && (m > 0), f[m - 1, n] - n f[m - 1, n - 1], 1]
Table[f[m, n], {m, 0, 10}, {n, 0, 10}] // MatrixForm   

we have the table for the first $10\times 10$ values

$$ \left[ \begin{array}{ccccccccccc} 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 0 & -1 & -2 & -3 & -4 & -5 & -6 & -7 & -8 & -9 \\ 1 & -1 & -1 & 1 & 5 & 11 & 19 & 29 & 41 & 55 & 71 \\ 1 & -2 & 1 & 4 & 1 & -14 & -47 & -104 & -191 & -314 & -479 \\ 1 & -3 & 5 & 1 & -15 & -19 & 37 & 225 & 641 & 1405 & 2661 \\ 1 & -4 & 11 & -14 & -19 & 56 & 151 & -34 & -1159 & -4364 & -11389 \\ 1 & -5 & 19 & -47 & 37 & 151 & -185 & -1091 & -887 & 6067 & 32251 \\ 1 & -6 & 29 & -104 & 225 & -34 & -1091 & 204 & 7841 & 14050 & -28419 \\ 1 & -7 & 41 & -191 & 641 & -1159 & -887 & 7841 & 6209 & -56519 & -168919 \\ 1 & -8 & 55 & -314 & 1405 & -4364 & 6067 & 14050 & -56519 & -112400 & 396271 \\ 1 & -9 & 71 & -479 & 2661 & -11389 & 32251 & -28419 & -168919 & 396271 & 1520271 \\ \end{array} \right] $$

as can we verify $f(m,n) = f(n,m)$