Consider a Markov chain on a state space V with size N, and let $\pi(v_j) = \sum_{v_i \in V} \pi (v_i)P(v_i,v_j)$ be the stationary distribution, where $P(v_i,v_j)$ is the transition probability. Moreover, let $X_0, X_1, \cdots$ be a realization of the chain with $X_i \sim \pi$ for $i \geq 0$.
Then, for any function $f : V \rightarrow R$, the following estimator is the asymptotically unbiased estimator of $(1/N) \sum_{i=1}^N f(v_i)$.
$\hat{\mu} = \dfrac{1}{\sum_{i=0}^{n-1} 1/\pi(X_i)} \sum_{i=0}^{n-1} \dfrac{f(X_i)}{\pi(X_i)}$
But I can't understand why the estimator above can be the asymptotically unbiasdness. I tried to find some documents related to this question, and it seems like we need to use delta method to approximate the expectation of $\hat{\mu}$. However, because I am not familiar with delta method, it makes me more confused. Could you please help me with proving the asymptotically unbiasedness?
Write $\hat\mu_n=\frac {F_n} {G_n}$ where $$F_n=\frac 1 n\sum_{i=0}^{n-1} \frac{f(X_i)}{\pi(X_i)}$$ $$G_n=\frac 1 n \sum_{i=0}^{n-1} \frac 1 {\pi(X_i)}$$ By weak law of large numbers $$F_n\to E\Big(\frac {f(X_i)}{\pi(X_i)}\Big)=\sum_{i=1}^N f(v_i)$$ $$G_n\to E\Big(\frac 1 {\pi (X_i)}\Big)=N$$ in probability. Note that $G_n\geq 1$ to get convergence of $$\frac {F_n}{G_n}\to\frac{\sum_{i=1}^N f(v_i)}N$$ in probability and note that $\frac {F_n}{G_n}$ is bounded to get convergence in mean ($L^1$) as required.
Sidenote: The only reason to use this estimator is if you don't have space to store unique values of $(X_0,\dots,X_n)$. Otherwise, averaging $f$ over the set of these unique values will give a more precise estimator which will be equal to the quantity estimated as soon as all states are visited (which can assured when the sum of $\pi$ over this set of unique values is equal to 1). Estimator in this problem would take many passes over the whole set for empirical probabilities of visiting a state to approach stationary.