  We present two new connections between the inhomogeneous stochastic higher
spin six vertex model in a quadrant and integrable stochastic systems from the
Macdonald processes hierarchy. First, we show how Macdonald $q$-difference
operators with $t=0$ (an algebraic tool crucial for studying the corresponding
Macdonald processes) can be utilized to get $q$-moments of the height function
$\mathfrak{h}$ in the higher spin six vertex model first computed in
arXiv:1601.05770 using Bethe ansatz. This result in particular implies that for
the vertex model with the step Bernoulli boundary condition, the value of
$\mathfrak{h}$ at an arbitrary point
$(N+1,T)\in\mathbb{Z}_{\ge2}\times\mathbb{Z}_{\ge1}$ has the same distribution
as the last component $\lambda_N$ of a random partition under a specific $t=0$
Macdonald measure. On the other hand, it is known that
$\mathbf{x}_N:=\lambda_N-N$ can be identified with the location of the $N$th
particle in a certain discrete time $q$-TASEP started from the step initial
configuration. The second construction we present is a coupling of this
$q$-TASEP and the higher spin six vertex model (with the step Bernoulli
boundary condition) along time-like paths providing an independent
probabilistic explanation of the equality of $\mathfrak{h}(N+1,T)$ and
$\mathbf{x}_N+N$ in distribution. Combined with the identification of averages
of observables between the stochastic higher spin six vertex model and Schur
measures (which are $t=q$ Macdonald measures) obtained recently in
arXiv:1608.01553, this produces GUE Tracy--Widom asymptotics for a discrete
time $q$-TASEP with the step initial configuration and special jump parameters.
