draw k times from n numbers with replacement, each number appears at least m times The objective is to draw $k$ times from $n$ numbers with replacement and with remembering the order. What is the probability $P(k,n,m)$ that every number ($1\ldots n$) is drawn at least $m$ times?
For $m=0$ the probability is $P(k,n,0)=1$.
For $m=1$ (1498003) the probability is $P(k,n,1)=\mathcal{S}_k^{(n)}n!/n^k$ in terms of the Stirling number of the second kind.
Is there an efficient way of calculating $P(k,n,m)$ for higher values of $m$?
Here's an inefficient way to calculate these probabilities in Mathematica by enumerating all possible drawing results:
f[k_Integer /; k >= 0, n_Integer /; n >= 1] := f[k, n] = 
   Reverse@Accumulate@Reverse@BinCounts[
        Min[BinCounts[#, {1, n + 1, 1}]] & /@ Tuples[Range[n], k],
        {0, Floor[k/n] + 1, 1}]/n^k;
P[k_Integer /; k >= 0, n_Integer /; n >= 1, m_Integer /; m >= 0] := 
  With[{L = f[k, n]}, If[m >= Length[L], 0, L[[m + 1]]]]

Test of the known formula for $m=1$:
With[{k = 8, n = 3},
  P[k, n, 1] == StirlingS2[k, n]*n!/n^k]
(* True *)


update: faster code based on solution
@GCab's solution gives an efficient recipe for calculating these probabilities. Simplifying the multinomials and rising factorials a bit,
$$
P(k,n,m)
=\frac{1}{n^k} \frac{k!}{(k-nm)!} \sum_{l_j\ge0\\l_1+l_2+\ldots+l_n=k-n m}\binom{k-nm}{l_1,l_2,\ldots,l_n}\frac{1}{(1+l_1)^{\bar{m}}(1+l_2)^{\bar{m}}\ldots(1+l_n)^{\bar{m}}}\\
=\frac{1}{n^k} \sum_{l_j\ge0\\l_1+l_2+\ldots+l_n=k-n m}\binom{k}{l_1+m,l_2+m,\ldots,l_n+m}\\
=\frac{1}{n^k} \sum_{a_j\ge m\\a_1+a_2+\ldots+a_n=k}\binom{k}{a_1,a_2,\ldots,a_n}
$$
In Mathematica,
P[k_, n_, m_] := 
  1/n^k * Total[(Multinomial @@ #)*(Multinomial @@ Tally[#][[All,2]]) & /@ 
    IntegerPartitions[k, {n}, Range[m, k-(n-1)*m]]]

The extra factor of Multinomial @@ Tally[#][[All,2]] in this formula counts how many permutation of a given integer partition exist, because IntegerPartitions does not enumerate partitions that differ only by permutation.
recursion
Still much faster is @GCab's recursive solution: in Mathematica,
$RecursionLimit = 10^4;
R[k_, n_, m_] /; k < n*m = 0;
R[k_, n_, 0] = n^k;
R[k_, 1, m_] = 1;
R[k_, n_, m_] := R[k, n, m] = 
  n*(R[k-1, n, m] + Binomial[k-1, m-1]*R[k-m, n-1, m])
P[k_, n_, m_] := R[k, n, m]/n^k

With this recursion we can go to large values of $k$ easily:
P[1000, 100, 7] // N // AbsoluteTiming
(* {0.245745, 2.4066*10^-8} *)

 A: In terms of combinatorial classes we get
$$\def\textsc#1{\dosc#1\csod}
\def\dosc#1#2\csod{{\rm #1{\small #2}}}
\textsc{SEQ}_{=n}(\textsc{SET}_{\ge m}(\mathcal{Z})).$$
This gives the EGF
$$\left(\exp(z)-\sum_{q=0}^{m-1} \frac{z^q}{q!}\right)^n.$$
The desired quantity is then given by
$$\bbox[5px,border:2px solid #00A000]{
\frac{1}{n^k} k! [z^k]
\left(\exp(z)-\sum_{q=0}^{m-1} \frac{z^q}{q!}\right)^n.}$$
 Addendum. For  computational purposes  it  is best  to use  a
recurrence like the one shown below.

X :=
(k, n, m) ->
1/n^k*k!*coeftayl((exp(z)-add(z^q/q!, q=0..m-1))^n,
                  z=0, k);

R :=
proc(k, n, m)
option remember;

    if n=1 then
        return `if`(k >= m, 1/k!, 0);
    fi;

    add(1/q!*R(k-q, n-1, m), q=m..k);
end;

XR := (k, n, m) -> 1/n^k*k!*R(k, n, m);


 Addendum, II.  If we are required to use  integers only we can
use the  recurrence below  (e.g. for  implementation in  a programming
language that does not provide fractions as a data type).

R2 :=
proc(k, n, m)
option remember;

    if n=1 then
        return `if`(k >= m, 1, 0);
    fi;

    add(binomial(k,q)*R2(k-q, n-1, m), q=m..k);
end;

XR2 := (k, n, m) -> 1/n^k*R2(k, n, m);

A: This the same as 
laying $k$ (distinguished, since you consider the extraction order) balls into $n$ (distinguished) bins
but allow me to change $k$ into $s$ so to keep $k$ for use as index, thus
laying $s$ (distinguished) balls into $n$ (distinguished) bins, each bin having at least $m$ balls,
or equivalently
writing words of length $s$, from the alphabet $\{1,\cdots,n\}$, with each character appearing at least $m$ times.
Without the bound on $m$ (i.e. for $m=0$), there are $n^s$ ways of doing that, each represented by one term in the expansion of 
$$
\left( {x_{\,1}  + x_{\,2}  +  \cdots  + x_{\,n} } \right)^{\,s} 
$$
Since
$$
\left( {x_{\,1}  + x_{\,2}  +  \cdots  + x_{\,n} } \right)^{\,s}  = \sum\limits_{\left\{ {\begin{array}{*{20}c}
   {0\, \le \,j_{\,k} }  \\   {j_{\,1}  + \,j_{\,2} +\, \cdots +\,j_{\,n} \, = \,s}  \\
\end{array}} \right.\;} {\left( \begin{array}{c}
 s \\   j_{\,1} ,\,j_{\,2} ,\, \cdots ,\,j_{\,n}  \\ 
 \end{array} \right)\;x_{\,1} ^{j_{\,1} } \;x_{\,2} ^{j_{\,2} } \; \cdots \;x_{\,n} ^{j_{\,n} } } 
$$
then the multinomial coefficient counts the number of "occupancy" histograms with distinguished balls
corresponding to the same histogram when the balls are undistinguished.
(the number of histograms with undist. balls would be the number of weak compositions of $s$ into $n$ parts).
Now we are looking for
$$
\begin{array}{l}
 P(s,n,m) = \frac{{N(s,n,m)}}{{n^{\,s} }} = \frac{1}{{n^{\,s} }}\sum\limits_{\left\{ {\begin{array}{*{20}c}
   {m\, \le \,j_{\,k} }  \\    {j_{\,1}  + \,j_{\,2}  +  \cdots  + j_{\,n} \, = \,s}  \\
\end{array}} \right.\;} {\left( \begin{array}{c}
 s \\   j_{\,1} ,\,j_{\,2} ,\, \cdots ,\,j_{\,n}  \\ 
 \end{array} \right)}  =  \\ 
  = \frac{1}{{n^{\,s} }}\sum\limits_{\left\{ {\begin{array}{*{20}c}
   {0\, \le \,l_{\,k} }  \\    {l_{\,1}  + \,l_{\,2}  +  \cdots  + l_{\,n} \, = \,s - nm}  \\
\end{array}} \right.\;} {\left( \begin{array}{c}
 s \\   m + l_{\,1} ,\,m + l_{\,2} ,\, \cdots ,\,m + l_{\,n}  \\ 
 \end{array} \right)}  =  \\ 
  = \frac{1}{{n^{\,s} }}\frac{{s!}}{{\left( {s - nm} \right)!}}\sum\limits_{\left\{ {\begin{array}{*{20}c}
   {0\, \le \,l_{\,k} }  \\    {l_{\,1}  + \,l_{\,2}  +  \cdots  + l_{\,n} \, = \,s - nm}  \\
\end{array}} \right.\;} {\frac{{l_{\,1} !\,l_{\,2} ! \cdots l_{\,n} !}}{{\left( {m + l_{\,1} } \right)!\,\left( {m + l_{\,2} } \right)!\, \cdots \,\left( {m + l_{\,n} } \right)!}}
\left( \begin{array}{c}
 s - nm \\   l_{\,1} ,\,l_{\,2} ,\, \cdots ,\,l_{\,n}  \\ 
 \end{array} \right)}  =  \\ 
  = \frac{1}{{n^{\,s} }}\frac{{s!}}{{\left( {s - nm} \right)!}}\sum\limits_{\left\{ {\begin{array}{*{20}c}
   {0\, \le \,l_{\,k} }  \\    {l_{\,1}  + \,l_{\,2}  +  \cdots  + l_{\,n} \, = \,s - nm}  \\
\end{array}} \right.\;} {\left( \begin{array}{c}
 s - nm \\   l_{\,1} ,\,l_{\,2} ,\, \cdots ,\,l_{\,n}  \\ 
 \end{array} \right)\frac{1}{{\left( {1 + l_{\,1} } \right)^{\,\overline {\,m\,} } \,\left( {1 + l_{\,2} } \right)^{\,\overline {\,m\,} } \, 
\cdots \,\left( {1 + l_{\,n} } \right)^{\,\overline {\,m\,} } }}}  \\ 
 \end{array}
$$
where $x^{\,\overline {\,k\,} }$ represents  the Rising Factorial
Although other manipulations of the above formula are possible, I do not see it possible to reach to a closed
formula, unless (maybe) through the q-Stirling ? (*)
We can however establish a recursive relation in $s$ that might be more practical in computing.
Let's take the "words" scheme and indicate by $N^ *  (s,n,m)$ the 
Number of words of length $s$, from the alphabet $\{1,cdots,n\}$, with one character repeated exactly$m-1$ times and the other at least $m$ times.
That is
$$
\begin{array}{l}
 N^ *  (s,n,m) = \sum\limits_{\left\{ {\begin{array}{*{20}c}
   {m\, \le \,j_{\,k} }  \\    {j_{\,1}  + \,j_{\,2}  +  \cdots  + \,j_{\,n - 1}  + m - 1\, = \,s}  \\
\end{array}} \right.\;} {\left( \begin{array}{c}
 s \\   j_{\,1} ,\,j_{\,2} ,\, \cdots ,\,j_{\,n - 1} ,\,m - 1 \\ 
 \end{array} \right)}  =  \\ 
  = \sum\limits_{\left\{ {\begin{array}{*{20}c}
   {m\, \le \,j_{\,k} }  \\    {j_{\,1}  + \,j_{\,2}  +  \cdots  + \,j_{\,n - 1} \, = \,s - m + 1}  \\
\end{array}} \right.\;} {\left( \begin{array}{c}
 s \\   j_{\,1} ,\,j_{\,2} ,\, \cdots ,\,j_{\,n - 1} ,\,m - 1 \\ 
 \end{array} \right)}  =  \\ 
  = \frac{{s!}}{{\left( {s - m + 1} \right)!\left( {m - 1} \right)!}}\sum\limits_{\left\{ {\begin{array}{*{20}c}
   {m\, \le \,j_{\,k} }  \\    {j_{\,1}  + \,j_{\,2}  +  \cdots  + \,j_{\,n - 1} \, = \,s - m + 1}  \\
\end{array}} \right.\;} {\left( \begin{array}{c}
 s - m + 1 \\   j_{\,1} ,\,j_{\,2} ,\, \cdots ,\,j_{\,n - 1}  \\ 
 \end{array} \right)}  =  \\ 
  = \left( \begin{array}{c}
 s \\   m - 1 \\ 
 \end{array} \right)N(s - m + 1,n - 1,m) \\ 
 \end{array}
$$
Then a word of length $s+1$, and parameters $n,m$ can be composed
 - either by adding any character to a word $(s,n,m)$,
 - or by adding the missing character to a word counted by $N^* (s,n,m)$.
Therefore we have
$$
\begin{array}{l}
 N(s + 1,n,m) = n\left( {N(s,n,m) + N^ *  (s,n,m)} \right) =  \\ 
  = n\left( {N(s,n,m) + \left( \begin{array}{c}
 s \\  m - 1 \\ 
 \end{array} \right)N(s - m + 1,n - 1,m)} \right) \\ 
 \end{array}
$$
with a suitable tuning of the initial conditions / range of the parameters.
Note(*)
I stumbled now on learning that N(s,n,m) is related to the m-associated Stirling numbers of 2nd kind
$$
{{N(s,n,m)} \over {n!}} = S_{\,m} \left( {s,n} \right)
$$
re. for instance to these links 3, 4, 5.
