$\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{star} \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{:=} [1]_q [2]_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.

Question: Is there some alternative (more recent?) reference where \eqref{star} is stated and easily referenced? Alternatively, someone who can see exactly where in the paper \eqref{star} is proven?

Garsia, A. M.; Gessel, I., Permutation statistics and partitions, Adv. Math. 31, 288-305 (1979). ZBL0431.05007.