Htykuut

Let $L$ be the collection of sequences indexed by $mathbb{N}$ taking values in $mathbb{R}$.
$$L = { (x_0, x_1, x_2, ldots ) : x_i in mathbb{R} }$$
$L$ will be a vector space with respect to componentwise addition and scalar multiplication.

Define a linear map $D$ on $L$ by
$$L ni x = (x_0, x_1, x_2, ldots ) quadmapstoquad Dx = (x_1, x_2, ldots ) in L$$
The result of applying $D - 1$ to a sequence $x$ is the finite difference of $x$.
$$left[ (D-1)x right]_n = x_{n+1} - x_n$$

If you apply $(D-1)^k$ to the sequence $u = (n^k)_{n=0}^infty$, you will find

$$left[ (D-1)^k u right]_n
= sum_{j=0}^{k}(-1)^{k-j}binom{k}{j} left[D^j uright]_n
= sum_{j=0}^{k}(-1)^{k-j}binom{k}{j} u_{n+j}
= sum_{j=0}^{k}(-1)^{k-j}binom{k}{j}(j+n)^k
$$
When $n = 0$, this is the LHS of the equality at hand.

To see LHS of your equality is indeed equal to the RHS, you can use following fact:

If $v$ is a sequence whose value is a polynomial in the index $n$, i.e.
$$v_n = a_p n^p + a_{p-1} n^{p-1} + cdots + a_0$$
for some constants $a_p, a_{p-1},ldots, a_0$, the
finite difference of $v$ will be a polynomial in $n$ with degree one less.
Furthermore, the leading coefficient is multiplied by a factor $p$. i.e.
$$left[(D-1)vright]_n = p a_p n^{p-1} + cdots$$

Apply these $k$ times to the sequence $u = (n^k)_{n=0}^infty$, you find
$$sum_{j=0}^{k}(-1)^{k-j}binom{k}{j}(j+n)^k = left[ (D-1)^k u right]_n = (D-1)^k n^k = k! n^0 = k!$$
which is a constant independent of $n$ and equal to RHS of the equality at hand.

A combinatorial answer.

The right-hand side is the number of permutations of ${0, 1, dots, k-1}$.

I'll let $k=5$ for concreteness.

We count the permutations in another way. Take a $5$-tuple drawn from ${0,1,dots,4}$ - that is, a member of $5^5$. This is probably not a permutations - for example, it might be $00000$.

Remove anything which does not contain a $4$. That is, remove anything which can be drawn from $4^5$.
This leaves all the permutations, but also leaves things like $44444$,
so in fact we should be removing anything which can be drawn from ${ 1, 2, 3, 4}$ and from ${ 0, 2, 3, 4}$ and so on.

That is, we remove $4^5$ things $binom{5}{1} = binom{5}{4}$ times.

OK, we've left all the permutations, but now we've removed some things more than once: we removed $00000$ because it doesn't contain a $4$ but also because it doesn't contain a $3$, for instance.
Generally, we've removed twice anything which has two things it doesn't contain, and we've removed three times anything which has three things it doesn't contain, and so on.
We need to add them back in enough times: we need to add once anything which has two things it doesn't contain, add twice anything which has three things it doesn't contain, and so on.

Add in $3^5$ things $binom{5}{2}$ times: namely, add in anything which does not contain a $0$ and does not contain a $1$, then add in anything which does not contain a $0$ and does not contain a $2$, then… then add in anything which does not contain a $3$ and does not contain a $4$.
We've re-added anything which has two things it doesn't contain, but we've re-added $binom{3}{2}$ times anything which has three things it doesn't contain, re-added $binom{4}{2}$ times anything which has four things it doesn't contain, and so on.

Repeat this procedure of adding and subtracting until we eventually finish.

As lulu said, this can be proved by inclusion-exclusion theorem. We have
begin{equation}
|cup_{i=1}^k A_i|=sum_{j=1}^k(-1)^{j+1}(sum_{1leq i_1<...<i_jleq k}|A_{i_1}cap...cap A_{i_j}|).
end{equation}
Now let $A_i$ be the subset of function from $[k]$ to $[k]$ such that $i$ is not an image. Then lhs of the equation is equal to $k^k-k!$, which means we subtrct the bijective functions from all functions. The $j$-th component of rhs is equal to $(-1)^{j+1}begin{pmatrix}k\jend{pmatrix}(k-j)^k$, which means that we choose $j$ elements such that they are not images, and the $k$ elements only have then $k-j$ possibilities to choose an image. Now let $tilde{j}=k-j$, and inserting this into the equation, you will find the solution of the given form.

More details: we obtain then
begin{align}
k^k-k!&=sum_{j=0}^{k-1}(-1)^{k-j+1}begin{pmatrix}k\k-jend{pmatrix}j^k\
&=-1sum_{j=0}^{k-1}(-1)^{k-j}begin{pmatrix}k\jend{pmatrix}j^k\
&=-1sum_{j=1}^{k}(-1)^{k-j}begin{pmatrix}k\jend{pmatrix}j^k+(-1)^{k-k}begin{pmatrix}k\kend{pmatrix}k^k\
&=-1sum_{j=1}^{k}(-1)^{k-j}begin{pmatrix}k\jend{pmatrix}j^k+k^k.
end{align}
Both sides subtrct $k^k$ and devided by $-1$, we obtain the result.

In this answer, it is shown that for any $l$,
$$
begin{align}
sum_{r=0}^nbinom{n}{r}(-1)^r(l-r)^n
&=sum_{r=0}^n(-1)^{n-r}binom{n}{r}(r-l)^n\
&=n!
end{align}
$$
The given identity is this identity using $l=0$.

Here is a variant based upon the coefficient of operator $[z^k]$ used to denote the coefficient of $z^k$ in a series. This way we can write e.g.
begin{align*}
[z^j](1+z)^k=binom{k}{j}qquadqquadtext{or}qquadqquad k![z^k]e^{jz}=j^k
end{align*}

We obtain
begin{align*}
sum_{j=1}^k(-1)^{k-j}j^kbinom{k}{j}
&=sum_{j=0}^infty(-1)^{k-j}k![z^k]e^{jz}[u^k](1+u)^ktag{1}\
&=(-1)^kk![z^k]sum_{j=0}^{infty}left(-e^{z}right)^j[u^k](1+u)^ktag{2}\
&=(-1)^kk^ktag{3}\
&=k![z^k]left(z+frac{z^2}{2!}+frac{z^3}{3!}+cdotsright)^ktag{4}\
&=k!
end{align*}
and the claim follows.

Comment:

• In (1) we apply the coefficient of operator twice and extend the lower and upper limit of the series without changing anything, since we add only zeros.




• In (2) we use the linearity of the coefficient of operator and do some rearrangements.




• In (3) we use the substitution rule with $u=-e^z$:

begin{align*}
A(z)=sum_{j=0}^infty a_jz^j=sum_{j=0}^infty z^j [u^k]A(u)
end{align*}

• In (4) we see that in each of the $k$ factors $(e^z-1)$ only $z^1$ provide a contribution to $z^k$.