Otra solución podría ser la siguiente:
Es un recuento doble de las funciones suryentes del conjunto $\{1,\dots,j\}$ al conjunto $\{1,\dots, p\}$ .
Por un lado, es evidente que estas funciones son $p!$ si $j=p$ mientras que no existen tales funciones si $j<p$ .
Vamos a contar estas funciones de una manera diferente:
Para $1\leq i\leq p$ definamos los conjuntos $$A_i:=\{\text{the functions from }\{1,\dots,j\}\text{ to }\{1,\dots,p\}\text{ which doesn't contain }\; i\;\text{ in their image}\}.$$
Tenemos entonces
$$\begin{align} |A_i|&=(p-1)^j\\ |A_i\cap A_l|&=(p-2)^j\\ &\vdots&\\ \left|\bigcap_{i=1}^p A_i\right|&=0.\end{align}$$
Estamos interesados en calcular lo siguiente: $$|\overline{A_1}\cap\overline{A_2}\cap\dots\cap\overline{A_p}|=|\overline{A_1\cup A_2 \cup\dots\cup A_p}|= p^j-|A_1\cup A_2 \cup\dots\cup A_p|= $$ $$p^j-\sum_{i=1}^p\left((-1)^{i+1}\sum_{1\leq j_1<\dots<j_i\leq p}\left|\bigcap_{k=1}^i A_{j_k}\right|\right)=\sum_{k=0}^p(-1)^{p+k}\binom{p}{k}k^j.$$ De lo cual se desprende la tesis.
EDITAR Dado que hace un año dediqué mucho tiempo a este problema, me gustaría dar otros dos centavos sobre el tema, y daría otra solución, o al menos un esbozo de ella:
Entonces, consideremos el operador diferencial que realiza: $$D(f(x)) = \frac{\mathrm d}{\mathrm dx}(x\cdot f(x)) = f(x) + x\cdot\frac{\mathrm df}{\mathrm dx}.$$ Inmediatamente se ve que $$D(x^m)=m\cdot x^m,$$ por lo que $$\sum_{k=0}^{p}(-1)^k \binom{p}{k}\,k^j = \left.D^j\left((1-x)^p\right)\right|_{x=1}.$$ Así que, (casi) inmediatamente, de nuevo sigue la tesis.
Permítanme decir una última cosa. Esta fórmula está estrictamente relacionada con Números Stirling del segundo tipo pero fíjate que la segunda prueba que te di no dice mucho sobre lo que pasa si $j>p$ , mientras que el primero encaja perfectamente incluso en esta situación.