$\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{:=} [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{*} is
stated and easily referenced? Alternatively, someone who can see exactly where in the paper \eqref{*} is proven?

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