I am trying to find a solution to the recurrence relation
$$g_{k}(n) = n g_{k-1}(n) - n g_{k-1} (n-1)$$
such that $g_k(n)$ is written as a summation over $g_0(n)$ for various $n$.
Attempt
I started by reviewing a simpler recurrence relation: $f_{k}(n) = f_{k-1}(n) - f_{k-1}(n-1)$. This is the standard finite difference formula for the $k$th order backward difference for $f_k(n)$ as a function of $n$. By defining the operator $\hat{\gamma}$ that acts as $$ \hat{\gamma} f_k(n) = f_k(n-1),$$ we can show \begin{align} f_k(n) & = f_{k-1}(n) - f_{k-1}(n-1) \\ & = (1- \hat{\gamma}) f_{k-1}(n)\\ & = (1-\hat{\gamma})^kf_0(n) \\ &= \sum_{j=0}^{k}\binom{k}{j}(-1)^{j}f_0(n-j). \end{align} which is the standard higher order finite difference formula. Then, I returned to the main case. For this case, I defined a second operator $\hat{\delta}$ that acts as $$ \hat{\delta} g_k(n) = n g_{k}(n).$$ The recurrence relation then becomes $$ g_k(n) = n g_k(n) - n g_k(n-1) = (\hat{\delta} - \hat{\gamma}\hat{\delta})g_{k-1}(n) = (\hat{\delta} - \hat{\gamma}\hat{\delta})^k g_{0}(n).$$ (Note that $\hat{\delta}$ and $\hat{\gamma}$ only act on $g_k(n)$ or $f_k(n)$ and don't act on $n$, i.e., $\hat{\delta} n = n \hat{\delta}$ and $\hat{\gamma} n = n \hat{\gamma}$) It would be great to apply the binomial theorem again, but the operators $\hat{\gamma}$ and $\hat{\delta}$ are non-commuting so another approach seems necessary.
Notes
It's likely this operator-based methods isn't the best way to approach the problem. I explored some generating function approaches also to no avail. Also, it might be worth trying to solve the somewhat simpler recurrence relation $h_{k}(n) = n h_{k-1}(n) - (n-1) h_{k-1}(n-1)$ though I ran into similar problems
It turns out generating functions were the way to go. Defining $$ G(t, n) = \sum_{k=0}^{\infty} t^k g_{k}(n),$$ with boundary conditions $G(0, n) = g_{0}(n)$ and $G(t, 0) = \sum_{k=0}^{\infty} t^{k} \delta_{k,0} = 1$, we can show that $$g_{k+1}(n) = n g_{k}(n) - n g_{k}(n-1),$$ has the generating function form $$ (1-nt) G(t, n) = g_0(n) - nt G(t, n-1).$$ Solving for $G(t, n)$ and applying the solution recursively, we end up with $$ G(t, n) = \sum_{\ell=0}^{n}(-1)^{\ell} g_{0}(n-\ell) t^{\ell} (n)_{\ell}\prod_{j=0}^{\ell}\frac{1}{1-(n-j)t},$$ where $(n)_{\ell} = n(n-1)\cdots (n-\ell+1)$ is the falling factorial. Expanding the product in $t$ and selecting the term proportional to $t^k$ gives us a possible final result for $g_{k}(n)$: $$\boxed{g_{k}(n) = \sum_{\ell=0}^k (-1)^{\ell} g_{0}(n-\ell) (n)_{\ell} h_{k-\ell}(n, n-1, \ldots, n-\ell),}$$ where $h_{m}(X_1, X_2,\ldots, X_N)$ is the complete homogeneous symmetric polynomial.
Alternatively, we can express the answer in terms of Stirling Numbers. In particular, we can show \begin{align} t^{\ell}\prod_{j=0}^{\ell}\frac{1}{1-(n-j)t}& = \frac{t^{\ell}}{(1-n t)^{\ell+1}} \sum_{m=0}^{\infty}(-1)^{m} \frac{t^{m}}{(1-nt)^{m}} h_{m}(1, 2,\ldots, \ell)\\ & = \sum_{m=0}^{\infty}(-1)^{m} \frac{t^{m+\ell}}{(1-nt)^{m+\ell+1}} S(m+\ell, \ell), \end{align} where in the final line we used the relationship between complete homogeneous symmetric polynomials and the Stirling numbers of the second kind. By repeatedly differentiating $1/(1-u)$ we can show $$ \frac{t^{q}}{(1-nt)^{q+1}} = \sum_{k=0}^{\infty} \binom{k}{q}n^{k-q} t^{k},$$ which makes the generating function become $$ G(t, n) = \sum_{k=0}^{\infty} t^{k} \sum_{\ell=0}^{k}(-1)^{\ell} g_{0}(n-\ell)\, (n)_{\ell} \sum_{m=0}^{\infty} (-1)^{m} n^{k-(m+\ell)} \binom{k}{m+\ell} S(m+\ell, \ell).$$ Taking the term proportional to $t^{k}$ and making the change of variables $r = m +\ell$, we ultimately find $$\boxed{g_{k}(n) = \sum_{\ell=0}^k g_{0}(n-\ell) (n)_{\ell}\sum_{r=0}^{\infty} (-1)^{r} n ^{k-r} \binom{k}{r} S(r, \ell) .}$$