# Number of permutations in $S_{a+b}$ with $\operatorname{maj}(\pi)=a$ and $\operatorname{maj}(\pi^{-1})=b$
$$\DeclareMathOperator\maj{maj}\DeclareMathOperator\inv{inv}$$Let $$t_{a,b}$$ be the numbers $$t_{a,b} \mathrel{:=} \lvert\{ \pi \in S_{a+b} : \maj(\pi)=a \text{ and } \maj(\pi^{-1})=b \}\rvert.$$ Here, $$S_{a+b}$$ denotes the set of permutations of $$1,2,\dotsc,a+b$$. By a result of Foata, one can also look at the pair of statistics $$(\maj, \inv)$$, and a few other combinations — these pairs of statistics will produce the same numbers.
Now, according to the OEIS entry A090806, it is proved by Garsia–Gessel that $$\begin{equation} \label{*} \tag{*} \sum_{a,b} t_{a,b} q^a t^b = \prod_{i,j \geq 1} \frac{1}{1-q^i t^j}. \end{equation}$$ I cannot see exactly where in their paper one can deduce this.
My attempt: I have tried to prove this myself (mainly by resorting to RSK, the Cauchy identity, and some symmetric function identities). This leads to the following (which appears in Stanley's EC2): $$\begin{equation} \sum_{n \geq 0} \frac{z^n}{(1-q)^n[n]_q!(1-t)^n [n]_t!} \sum_{\pi \in S_n} t^{\maj(\pi)} q^{\maj(\pi^{-1})} = \prod_{i,j \geq 0} \frac{1}{1-z q^i t^j}, \end{equation}$$ where $$[n]_q! \mathrel{:=} _q _q \dotsm [n]_q$$, and $$[n]_q = 1+q+q^2+\dotsb + q^{n-1}$$. However, I do not see some short way to deduce the above generating function from this.