Ruby, 101 98 bytes

->n,k{[w=[]].product(*(1..n).map{|z|(0..k).map{|x|[z]*x}}){|x|w|=[*x.flatten.permutation]};w.size}

Try it online!

R, 86 bytes

\(n,k)sum(apply(expand.grid(rep(list(0:k),n)),1,\(r)gamma(sum(r)+1)/prod(gamma(r+1))))

Attempt This Online!

Uses the formula from the question.

05AB1E, 10 bytes

LIи怜€`Ùg

Brute-force approach. Output includes the empty item.

Try it online or verify the smaller test cases.

Explanation:

L           # Push a list in the range [1, first (implicit) input]
 Iи         # Repeat each value the second input amount of times
   æ        # Take the powerset of this list (including empty list)
    €       # Map over each inner list:
     œ      #  Get all permutations
      €`    # Flatten the list of list of lists of integers one level down
        Ù   # Uniquify this list of lists
         g  # Pop and push its length
            # (which is output implicitly as result)

JavaScript (ES6), 59 bytes

Expects (n)(k). Not counting the empty string.

n=>g=(k,i=n)=>i--&&g(k,i)-((g[i]|=0)<k&&g[i]++-g(k)-g[i]--)

Try it online!

Commented

n =>            // outer function taking the alphabet size n
g = (           // inner recursive function taking:
  k,            //   the occurrence limit k
  i = n         //   a counter i initialized to n
) =>            //
i-- &&          // decrement i; if it was not 0:
  g(k, i)       //   do a recursive call to try the next character type
  - (           //   and subtract:
    (g[i] |= 0) //     if still undefined, set g[i] to 0
    < k &&      //     if g[i] is less than k:
      g[i]++    //       increment g[i]
      - g(k)    //       subtract g(k) (effectively adding this character)
      - g[i]--  //       subtract g[i] and decrement it
  )             //   NB: g[i]++ - g[i]-- is -1
                //       and -(g[i]++ - g(k) - g[i]--) is g(k) + 1

Haskell, 66 63 bytes

-3 bytes thanks to xnor

a?(b:c)=1+[]?(a++[b-1|b>1]++c)+(b:a)?c
_?_=0

((?)[].).replicate

Try it online!

Doesn't count the empty string.

This loosely translates my other answer.

The only reason it is worth posting is that it beats the current best Haskell. However, I'm not sure this is optimal, even for this strategy.

Vyxal, 7 bytes

ɾẋfÞxUL

Try it Online!

Inefficient method, as it generates all strings. Takes n then k. Outputs a(n, k) - 1

Explained

ɾẋfÞxUL­⁡​‎‎⁡⁠⁡‏‏​⁡⁠⁡‌⁢​‎‎⁡⁠⁢‏‏​⁡⁠⁡‌⁣​‎‎⁡⁠⁣‏‏​⁡⁠⁡‌⁤​‎‎⁡⁠⁤‏⁠‎⁡⁠⁢⁡‏‏​⁡⁠⁡‌⁢⁡​‎‎⁡⁠⁢⁢‏‏​⁡⁠⁡‌⁢⁢​‎‎⁡⁠⁢⁣‏‏​⁡⁠⁡‌­
ɾ        # ‎⁡Range [1, n]
 ẋ       # ‎⁢Repeated k times
  f      # ‎⁣And flattened. This gives us a list `nk` with each character in the alphabet `n` repeated `k` times.
   Þx    # ‎⁤Get all combinations without replacement of `nk`. This is `combinations(nk, size = 1) ++ combinations(nk, size = 2) ++ ... ++ combinations(nk, size = n * k)`
     U   # ‎⁢⁡Remove duplicates
      L  # ‎⁢⁢And get the length of the resulting list
💎

Created with the help of Luminespire.

Japt -x, 14 bytes

There has to be a better way. But, it's been a long weekend and this is the best I can manage.

Takes k as the first input. Excludes the empty string.

ÆVoÃc
£áYÄ â Ê

Try it

ÆVoÃc\n£áYÄ â Ê     :Implicit input of integers U=k & V=n
Æ                   :Map the range [0,U]
 Vo                 :  Range [0,V]
   Ã                :End map
    c               :Flatten
     \n             :Assign that to U
       £            :Map each 0-based index, Y
        á           :  Permutations of U of length
         YÄ         :    Y+1
            â       :  Deduplicate
              Ê     :  Length
                    :Implicit output of sum of resulting array

Uiua, 16 bytes

/+≡⌟(⧻◴⧅≠)+1°⊏⊚▽

outputs a(n,k)-1, and hits memory limits for n>3, k>3

Try it!

/+≡⌟(⧻◴⧅≠)+1°⊏⊚▽
               ⊚▽  # k copies of n characters
          +1°⊏     # 1..length of the string
  ≡⌟(    )         # for each value in the range:
     ⧻◴⧅≠          #   get the number of unique permutations with that length
/+                 # sum

Charcoal, 30 bytes

≔⊕ＮθＩΣＥＸθＮ÷Π…·¹↨↨ιθ¹ΠＥ↨ιθΠ…·¹λ

Attempt This Online! Link is to verbose version of code. Takes inputs in the order k, n. Explanation:

  Ｎ                             Input `k` as a number
 ⊕                              Incremented
≔  θ                            Store in variable
        θ                       Variable `k+1`
       Ｘ                        Raised to power
         Ｎ                      Input `n` as a number
      Ｅ                         Map over implicit range
                 ι              Current value
                ↨ θ             Converted to base `k+1`
               ↨   ¹            Sum of digits
           Π…·¹                 Factorial
          ÷                     Integer divided by
                       ι        Current value
                      ↨ θ       Converted to base `k+1`
                     Ｅ          Map over digits
                             λ  Current digit
                         Π…·¹   Factorial
                    Π           Take the product
     Σ                          Take the sum
    Ｉ                           Cast to string
                                Implicitly print

Imagine if I had the builtins to write Print(Cast(Sum(Map(CartesianPower(InclusiveRange(0, InputNumber()), InputNumber()), IntDivide(Factorial(Sum(i)), Product(Factorial(i)))))));. Still would be 17 bytes though...

28 bytes for a (relatively slow) port of @WheatWizard's Haskell answer:

Ｎθ⊞υＥＮθＦυＦ⌕ＡＸ⁰ι⁰⊞υＥι⁻λ⁼μκＩＬυ

Try it online! Link is to verbose version of code. Takes inputs in the order k, n. Explanation:

Ｎθ

Input k.

⊞υＥＮθＦυ

Start a breadth-first traversal with a list containing n copies of k.

Ｆ⌕ＡＸ⁰ι⁰

Enumerate the non-zero indices of the list.

⊞υＥι⁻λ⁼μκ

For each such index, push the result of decrementing that index to the traversal list.

ＩＬυ

Output the total of lists traversed.

JavaScript (V8), 140 138 bytes

a=(n,k,l=[],f=n=>n?f(n-1)*n:1)=>n?[...Array(k+1).keys()].reduce((s,v)=>s+a(n-1,k,[...l,v]),0):f(eval(l.join`+`))/l.reduce((s,v)=>s*f(v),1)

Try it online!

A first attempt that can probably be simplified. This code runs slowly due to eval() being called after every recursion, but it saves some bytes.

Haskell, 68 bytes

import Data.List
n#k=length$nub$inits=<<permutations([1..k]>>[1..n])

Try it online!

Naïve solution that generates all possible strings and counts them.

[1..k]>>[1..n] is monad magic that concatenates k copies of [1..n] which is the entire pool of symbols a word is allowed to have. We pass that to permutations which will rearrange those symbols every possible way, then we use inits=<< to generate every prefix of those words, then we use nub to remove duplicates, and finally count the length.

Haskell + hgl, 44 bytes

f x=1+fo[f$a++b-1:c|(a,b:c)<-spe x,b>0]

f<<rl

Explanation

For input \$n\$ and \$k\$ we make a list of \$n\$ copies of \$k\$ then we count the ways to decrement the list without passing below zero.