  Let $(h_t)_{t\in \mathbb{R}}$ be the horocycle flow acting on
$(M,\mu)=(\Gamma \backslash \text{SL}(2,\mathbb{R}),\mu)$, where $\Gamma$ is a
co-compact lattice in $\text{SL}(2,\mathbb{R})$ and $\mu$ is the homogeneous
probability measure locally given by the Haar measure on
$\text{SL}(2,\mathbb{R})$. Let $\tau\in W^6(M)$ be a strictly positive function
and let $\mu^{\tau}$ be the measure equivalent to $\mu$ with density $\tau$. We
consider the time changed flow $(h_t^\tau)_{t\in \mathbb{R}}$ and we show that
there exists $\gamma=\gamma(M,\tau)>0$ and a constant $C>0$ such that for any $
f_0, f_1, f_2\in W^6(M)$ and for all $0=t_0<t_1<t_2$, we have $$\ \left|\int_M
\prod_{i=0}^{2} f_i\circ h^\tau_{t_i} d \mu^\tau -\prod_{i=0}^{2}\int_M f_i d
\mu^\tau \right|\leq C \left(\prod_{i=0}^{2} \|f_i\|_6\right) \left(\min_{0\leq
i<j\leq 2} |t_i-t_j|\right)^{-\gamma}.$$ With the same techniques, we establish
polynomial mixing of all orders under the additional assumption of $\tau$ being
fully supported on the discrete series.
