I get $$f(a,b) = \frac{ (1-np)^b n!}{(n-a)!} \sum_{j=0}^b \binom{b}{j} \left\{ j \atop a \right\} \left(\frac{p}{1-np}\right)^j,$$ where $\left\{ j \atop a \right\}$ is a Stirling number of the second kind. I tested the formula in Mathematica against your recurrence for different values of $n$ and $p$, and they agree. Unfortunately I'm having trouble pasting the Mathematica input and output that verifies that agreement on this site without it turning into gobbledygook.

To get the formula I used Theorem 6 (which is actually due to Neuwirth) in my paper "On Solutions to a General Combinatorial Recurrence," *Journal of Integer Sequences* 14 (9): Article 11.9.7, 2011. The paper is about solution techniques for solving certain multivariate recurrences. Your recurrence happens to be in one of the forms for which the techniques work.

There might be a way to simplify the summation involving binomial coefficients and Stirling numbers, but I don't see it right now.

*Added*: The formula I used is the following.

**Theorem**. Suppose $R(n,k)$ satisfies the recurrence $$(\alpha(n-1) + \beta k + \gamma)R(n-1),k) + (\beta' + \gamma')R(n-1,k-1) + [n=k=0].$$ Then $$ R(n,k) = \left(\prod_{i=1}^k (\beta' i + \gamma') \right) \sum_{i=0}^n \sum_{j=0}^n \left[ n \atop i \right] \binom{i}{j} \left\{ j \atop k \right\} \alpha^{n-i} \beta^{j-k} \gamma^{i-j}.$$
(Here, $0^0$ is taken to be $1$.)

For the OP's recurrence, we have $\alpha = 0, \beta = p, \gamma = 1-np, \beta' = -p,$ and $\gamma' = (n+1)p$.