Here is a derivation of $(\ast)$ from the displayed equation
$$\sum \frac{z^n}{(1-q)^n [n]_q! (1-t)^n [n]_t!} \sum_{\pi \in S_n} q^{\mathrm{maj}(\pi)} t^{\mathrm{maj}(\pi^{-1})} = \prod_{i,j \geq 0} \frac{1}{1-zt^i q^j}.$$
Taking the coefficient of $z^n$ on both sides, we have
$$\frac{\sum_{\pi \in S_n} q^{\mathrm{maj}(\pi)} t^{\mathrm{maj}(\pi^{-1})}}{\prod_{i=1}^n (1-q^i) \prod_{j=1}^n (1-t^j)} = h_n(\{t^i q^j : i,j \geq 0\}).$$
Here the RHS is the complete homogenous symmetric function evaluated on the set of monomials $\{t^i q^j : i,j \geq 0\}$.
Now, $h_n(1, u_1, u_2, u_3, \ldots) = \sum_{k=0}^n h_k(u_1, u_2,\ldots).$ So we can rewrite this RHS to get
$$\frac{\sum_{\pi \in S_n} q^{\mathrm{maj}(\pi)} t^{\mathrm{maj}(\pi^{-1})}}{\prod_{i=1}^n (1-q^i) \prod_{j=1}^n (1-t^j)} = \sum_{k=0}^n h_k(\{ q^i t^j : i,j \geq 0,\ (i,j) \neq (0,0) \}). \qquad (\clubsuit)$$
According to the OEIS entry, the quantity $\# \{ \pi \in S_n : \mathrm{maj}(\pi) = a,\ \mathrm{maj}(\pi^{-1}) = b \}$ stabilizes at $t_{ab}$ as $n \to \infty$; the OP's mention of $n=a+b$ is just a particular value in the stable range. So we can take the limit of both sides of $(\clubsuit)$ as $n \to \infty$ to get
$$\frac{\sum_{a,b \geq 0} t_{ab} q^a t^b}{\prod_{i=1}^{\infty} (1-q^i) \prod_{j=1}^{\infty} (1-t^j)} =
\sum_{k=0}^{\infty} h_k(\{ q^i t^j : i,j \geq 0,\ (i,j) \neq (0,0) \}) = \prod_{\substack{i,j \geq 0 \\ (i,j) \neq (0,0)}} \frac{1}{1-q^i t^j} . \qquad (\diamondsuit)$$
Now cancel common factors from both sides of $(\diamondsuit)$ to get the claim.
Here is another approach, which suggests that something more interesting may be going on. Recall that the RS correspondence is a bijection between the symmetric group $S_n$ and pairs of SYT $(T,U)$ of the same shape $\lambda$, where $|\lambda| = n$. We define $\mathrm{maj}(T)$ for an SYT $T$ to be the sum of those $i$ such that $i$ occurs in a strictly higher row of $T$ than $i+1$ does. Then, if RS maps $\pi$ to $(T,U)$, we have
$\mathrm{maj}(\pi)= \mathrm{maj}(T)$ and $\mathrm{maj}(\pi^{-1})= \mathrm{maj}(U)$.
Define $f^{\lambda}(q)$ to be the sum, over $SYT$ of shape $\lambda$, of $q^{\mathrm{maj}(T)}$. So
$$\sum_{\pi \in S_n} q^{\mathrm{maj}(\pi)} t^{\mathrm{maj}(\pi^{-1})} = \sum_{|\lambda| = n} f^{\lambda}(q) f^{\lambda}(t).$$
Given a permutation $\mu$, and $n > |\mu| + \mu_1$, let $\mu[n]$ be the partition $(n-|\mu|, \mu_1, \mu_2, \ldots, \mu_k)$. So every partition of $n$ is of the form $\mu[n]$ for a unique $\mu$ and that, for any $\mu$, the partition $\mu[n]$ is well defined for $n$ large enough.
It is easy to see that $\lim_{n \to \infty} f^{\mu[n]}(q)$ exists. Set $f^{\mu[\infty]}(q)$ to be $\lim_{n \to \infty} f^{\mu[n]}(q)$.
So we have
$$\lim_{n \to \infty} \sum_{\pi \in S_n} q^{\mathrm{maj}(\pi)} t^{\mathrm{maj}(\pi^{-1})} = \sum_{\mu} f^{\mu[\infty]}(q) f^{\mu[\infty]}(t).$$
On the other hand, the Cauchy identity gives
$$\prod_{i,j \geq 1} \frac{1}{1-q^i t^j} = \sum_{\mu} s_{\mu}(q,q^2, \cdots) s_{\mu}(t,t^2, \cdots).$$
So here is the weird thing: It turns out that $f^{\mu[\infty]}(q)$ and $s_{\mu}(q)$ are the same thing! This seems like it should have a combinatorial proof. We can clearly think of $f^{\mu[\infty]}(q)$ as a generating function for tableau of shape $\mu$ with distinct entries, counted by a variant of major index. And $s_{\mu}(q)$ is the generating function for semistandard tableau of shape $\mu$, counted by weight. It feels like there should be an easy bijection here.
Well, I couldn't find one. But it isn't hard to prove the equality using hook length formulas. Let $m = |\mu|$ and let $h_1$, $h_2$, ..., $h_m$
be the hook lengths of $\mu$. Set $M = \sum (i-1) \mu_i$. We have
$$s_{\mu}(q) = \frac{q^M}{\prod_i (1-q^{h_i})} .$$
On the other hand, the hook lengths of $\mu[n]$ are $h_1$, $h_2$, ..., $h_m$ together with an $n-m$ additional hook lengths $S$. The exact set $S$ doesn't matter; what is important is $\{1,2,\ldots, n-\mu_1-m \} \subseteq S \subseteq \{1,2,\ldots, n \}$. So, by a formula of Stanley (EC2, Cor. 7.21.5), we have
$$f^{\mu[n]}(q) = \frac{q^M (1-q)(1-q^2) \cdots (1-q^N)}{\prod_{s \in S} (1-q^s) \prod_i (1-q^{h_i})}.$$
In the limit as $n \to \infty$, both $(1-q)(1-q^2) \cdots (1-q^N)$ and $\prod_{s \in S} (1-q^s)$ approach $\prod_{k=1}^{\infty} (1-q^k)$, so they cancel and we are left with
$$ f^{\mu[\infty]}(q) = \frac{q^M}{\prod_i (1-q^{h_i})} = s_{\mu}(q,q^2, \cdots).$$
