First, let's note that there's no reason why the number of ``bags'' should be the same as $|A|$. The problem is nicer if we let the number of bags $m$ be independent of $n$.

As Brendan noted, the problem can be viewed as counting orbits of a certain set under an action of $S_m\times S_n$. However, the action of $S_m$ just reduces counting sequences to counting multisets, so if we can count orbits of multisets in some other way, we don't need this action and we can just count orbits of multisets under an action of $S_n$.

More generally, we could
start with a finite group $G$ acting on a set $T$. Then $G$ acts in a natural way on multisets of elements of $T$. We would like to count orbits under this action of multisets of $m$ elements of $T$.
We can do this using (the weighted form of) Burnside's lemma.

it is sufficient to count, for each $g\in G$, the number of $m$-multisets fixed by $g$. But a multiset is fixed by $g$ if and only if it is a multiset union of orbits of $g$ on $T$. In our example, $G$ is the symmetric group $S_n$ and $T$ is also $S_n$, and $G$ acts on $T$ by multiplication. The orbits of $g$ are easy to determine, as described below.

Let's look first at the case $n=2$. There are two permutations of $\{1,2\}$, which I'll denote as $\iota=\binom{1\ 2}{1\ 2}$ and $\tau=\binom{1\ 2}{2\ 1}$. The generating function for multisets of these permutations is $1/(1-x)^2$ (where the weight of a multiset of size $j$ is $x^j$) and this is the same as the generating function for multisets fixed by $\iota$. There is only one orbit of $\tau$, which is all of $S_2$, so multisets fixed by $\tau$ have the same number of copies of $\iota$ and $\tau$, and the generating function for these is $1/(1-x^2)$. So by Burnside's lemma, the number of orbits of multisets is
$$\frac12\left(\frac{1}{(1-x)^2} +\frac{1}{1-x^2}\right)
=\frac{1}{(1-x)(1-x^2)}.
$$

In the general case we need to find the generating function for multisets of elements of $S_n$ fixed by a permutation $\pi\in S_n$. A multiset fixed by $\pi$ must be a multiset union of orbits $\pi$ on $S_n$. Let $o(\pi)$ be the order of $\pi$, that is, the least $k$ such that $\pi^k$ is the identity permutation. Then every orbit of $\pi$ on $S_n$ has $o(\pi)$ elements and there are $n!/o(\pi)$ orbits. Thus the generating function for multisets of these orbits, which is the generating function for multisets of elements of $S_n$ fixed by $\pi$, is
$$
\frac{1}{(1-x^{o(\pi)})^{n!/o(\pi)}}
$$
and thus the generating function for orbits of multisets under the action of $S_n$
is
$$\frac{1}{n!} \sum_{\pi\in S_n} \frac{1}{(1-x^{o(\pi)})^{n!/o(\pi)}}.$$
(If $\pi$ has cycles of length $\lambda_1,\dots, \lambda_k$ then $o(\pi)$ is the least common multiple of $\lambda_1,\dots, \lambda_k$.)

Since $o(\pi)$ depends only on the cycle type of $\pi$ we can get a sum with fewer terms by summing over cycle types, which are partitions of $n$. For a partition $\lambda=(\lambda_1, \lambda_2,\dots, \lambda_k)$ of $n$ in which $i$ occurs $m_i$ times as a part, let $z_\lambda=\prod_i i^{m_i}m_i!$, so there are $n!/z_\lambda$ permutations in $S_n$ of cycle type $\lambda$. Let $o(\lambda) = {\rm lcm}(\lambda_1,\dots, \lambda_k)$. Then the generating function may be written
$$\sum_{\lambda\vdash n} \frac{1}{z_\lambda(1-x^{o(\lambda)})^{n!/o(\lambda)}}$$
where the sum is over all partitions of $n$. The answer to the OP's original question is the coefficient of $x^n$ in this generating function; here any partitions $\lambda$ with $o(\lambda)>n$ can be ignored.

For $n=3$ the generating function is
\begin{multline*}
\frac16\left(\frac{1}{(1-x)^6} +\frac{3}{(1-x^2)^3}+\frac{2}{(1-x^3)^2}
\right)
= \frac{1+{x}^{2}+3{x}^{3}+5{x}^{4}+{x}^{5}+{x}^{6}}{(1-x)(1-x^2)^3(1-x^3)^2}\\
= 1+x+5{x}^{2}+10{x}^{3}+24{x}^{4}+42{x}^{5}+83{x}^{6}+132{x}^{7}+\cdots
\end{multline*}
(OEIS sequence A037240).
The coefficient of $x^3$ is 10, which agrees with bof.

For $n=4$ the generating function is
\begin{multline*}
\frac1{24} \left( 1-x \right) ^{-24}+\frac38 \left( 1-{x}^{2} \right) ^{-12}
+\frac13 \left( 1-{x}^{3} \right) ^{-8}+\frac14 \left( 1-{x}^{4} \right) ^
{-6}\\
=1+x+17{x}^{2}+111{x}^{3}+762{x}^{4}+4095{x}^{5}+19941{x}^{6}
+84825{x}^{7}+\cdots
\end{multline*}
and for $n=5$ it is
\begin{multline*}
\frac {1}{120} \left( 1-x \right) ^{-120}+{\frac {5}{24}}
\left( 1-{x}^{2} \right) ^{-60}+\frac16 \left( 1-{x}^{3} \right) ^{-40}
+\frac16 \left( 1-{x}^{6} \right) ^{-20}\\
+\frac14 \left( 1-{x}^{4} \right)
^{-30}+\frac15 \left( 1-{x}^{5} \right) ^{-24}\\
=1+x+73{x}^{2}+2467{x}^{3}+76044{x}^{4}\\+1876255{x}^{5}+39096565
{x}^{6}+703593825{x}^{7}+\cdots
\end{multline*}

The first few values of the numbers $a_n$ corresponding to the OP's original problem are $a_1=1$, $a_2=2$, $a_3=10$, $a_4 = 762$, $a_5=1876255$, $a_6=274382326290$, $a_7=3265588553925722827, \dots$
It's easy to go much further ($n=30$ takes just over a second in Maple) but the numbers get very big quickly.

A different formula for counting orbits of multisets (useful in some other problems, but more complicated for this problem) is described here and used here.
Some related papers are
Marko V. Jarić and Joseph L. Birman,
New algorithms for the Molien function.
J. Mathematical Phys. 18 (1977), no. 7, 1456-1458,
and
Marko V. Jarić and Joseph L. Birman, Calculation of the Molien generating function for invariants of space groups.
J. Mathematical Phys. 18 (1977), no. 7, 1459-1465. The second of these papers contains the generating function for $n=3$ given above, though I'm not sure that the interpretation is the same.

11more comments