The anwer is **yes**, even if $p: \mathbb{R}^n \to \mathbb{R}$ is only a bounded continuous function.

**Preliminary observations:**

Let us first consider the case $p=0$.

Since $q$ is globally Lipschitz on $\mathbb{R}^n$, say with Lipschitz constant $L$, all solutions of the ODE $\dot x = q(x)$ on $\mathbb{R}^n$ exist globally and, if $(\varphi_t)_{t \in \mathbb{R}}$ denotes the induced flow in $\mathbb{R}^n$, each function $\varphi_t: \mathbb{R}^n \to \mathbb{R}^n$ is Lipschitz continuous with constant $e^{|t|L}$ (this follows from Grönwall's lemma).

Since we assumed $p=0$, the solution of your PDE is given by
$$
  f(t,x) = f_0(\varphi_t(x)) \tag{*}\label{1}
$$
for each $t$ and $x$. 

Hence, if $f_0$ decays faster than every polynomial, so does $f(t, \cdot)$ for each fixed time $t$.

Let us now phrase this in a more functional analytic language: For each integer $k \ge 0$, consinder the weight function $w_k: \mathbb{R}^n \to [1,\infty)$ given by $w_k(x) =1 + \lvert x \rvert^k$, and consider the Banach space $E_k$ of all continuous functions $f: \mathbb{R}^n \to \mathbb{R}$ which satisfies $w_k f \in C_0(\mathbb{R}^n; \mathbb{R})$; we endow the space $E_k$ with the weighted supremum norm $\|\cdot\|_{w_k}$ given by
$$
  \|f\|_{w_k} = \|w_k f\|_\infty.
$$

It follows from $\eqref{1}$ that the *Koopman group* $(T_k(t))_{t \in \mathbb{R}}$ given by
$$
  T_k(t)f = f \circ \varphi_t
$$
for each $f \in E_k$ is a well-defined one-parameter group of bounded linear operators on $E_k$. Moreover, the orbits of the group are continuous with respect to the weak topology on $E_k$ (use that (i) $E_k$ is isomorphic to $C_0(\mathbb{R}^n; \mathbb{R})$ via multiplication with $w_k$, that (ii) for each $x \in \mathbb{R}^n$, the mapping $t \mapsto \varphi_t(x)$ is continuous, and (iii) the Riesz representation theorem for measures on $C_0(\mathbb{R}^n; \mathbb{R})$).
Hence, by a standard result from $C_0$-semigroup theory, the the group $(T_k(t))_{t \in \mathbb{R}}$ is even continuous with respect to the strong operator topology, i.e., it is a $C_0$-group on $E_k$.

It is trivial, but important to observe, that the groups $(T_k(t))_{t \in \mathbb{R}}$ act consistently on the spaces $E_k$.


Using these preliminary observations, we can show:

**Theorem.** Let $p$ be continuous and bounded.

(i) For each $f_0 \in C_0(\mathbb{R}^n; \mathbb{R})$ there exists a unique mild solution $f: \mathbb{R} \to C_0(\mathbb{R}^n; \mathbb{R})$ of the Cauchy problem
$$
  \dot f = \langle q, \nabla \rangle f(t) + p, \qquad f(0) = f_0.
$$

(ii) If $f_0$ decays faster than every polynomial and $f: \mathbb{R} \to C_0(\mathbb{R}^n; \mathbb{R})$ denotes the mild solution from (i), then for each fixed time $t \in \mathbb{R}$, the function $f(t)$ also decays faster (in space) than every polynomial.

*Proof.* 
We note that multiplication with $p$ is a bounded linear operator on $E_k$ for each $k$. Hence, by standard perturbation theory for $C_0$-semigroups, the differential operator on the right hand side of the PDE generates a $C_0$-semigroup on $E_k$. Applying this to $k=0$ yields assertion (i).

Moreover, these $C_0$-semigroups are consistent on the spaces $E_k$ (this follows, for instance, from the Dyson-Phillips series and the fact that each space $E_{k+1}$ embedds continuously into $E_k$).

If $f_0$ decay faster then every polynomial, then we have $f_0 \in E_k$ for each $k$. Hence, for each time $t$ and each $k$, we have $T_0(t)f_0 = T_k(t)f \in E_k$. Thus, for each time $t$, the function $T_0(t)f_0$ decays faster in space than any polynomial. $\square$