Let $T_n$ denote the $n^\text{th}$ Chebyshev polynomial, defined by the recursion
\begin{align} T_0(x) &= 1,\\ T_1(x) &= x,\\ T_n(x) &= 2x \cdot T_{n-1} (x) - T_{n-2} (x). \end{align}
Now, supposed that I have $m \geqslant n$ samples $X_1, \ldots, X_m$ from a probability measure $p$, and I want to construct an unbiased estimator of the quantity
\begin{align} \tau_n = T_n ( \mathbf{E}_p [X]). \end{align}
I know that this is possible: one can easily construct rudimentary unbiased estimators for $\{ \mu_d = \mathbf{E}_p [X]^d \}_{0 \leqslant d \leqslant n}$ by writing $\mathbf{E}_p [X]^d = \mathbf{E}_{p^{\otimes d}}[X_1 X_2 \cdots X_d]$ , and then use a Rao-Blackwell-type argument to average over permutations of the various $X_i$ to obtain estimators of minimal variance. By linearity, it follows that unbiased estimation of $\tau_n$ is also possible.
Now, I could apply this linearity algorithmically as well, i.e. write out $T_n$ in terms of monomials, estimate the monomials as applied to $\mathbf{E}[X]$, and then aggregate them into a final estimator for $\tau_n$. However, this could get a bit hefty, particularly if I am interested in estimating $\tau_n$ for several values of $n$.
As such, I want to ask whether there is a known way of forming unbiased estimates of $\tau_n$ recursively. I know that this can be done for the sequence $\mu_n = ( \mathbf{E}_p [X])^n$, using e.g. the Newton-Girard identities, which give a highly efficient solution. My hope is that the recursive definition of the Chebyshev polynomials will enable a spiritually similar solution for $\tau_n$, but I am open to the possibility that this is too optimistic.
