In exactly what sense do MCMC draws approximate the target? Background
We want to sample from some intractable density $\pi(\theta)$. Using an MCMC algorithm, we generate a sample of draws $\{\theta_i\}_{i=1}^N$ from a Markov chain that has $\pi(\theta)$ as its invariant distribution. We typically use the draws $\{\theta_i\}_{i=1}^N$ to "approximate" $\pi(\theta)$ in the following ways:


*

*Use sample averages $(1/N)\sum_{i=1}^Nf(\theta_i)$ to approximate integrals like
$$\mathbb{E}[f(\theta)] = \int f(\theta)\pi(\theta)d\theta;$$

*Generate histograms and kernel density estimates to visualize $\pi(\theta)$ or its marginals;

*Compute the ECDF $\hat{F}(\theta)=(1/N)\sum_{i=1}^N\mathbf{1}_{\theta_i\leq\theta}$ to estimate quantiles of $\pi(\theta)$.


If $\{\theta_i\}_{i=1}^N$ were direct, iid draws from $\pi(\theta)$, then there is theory that justifies all of these approximations. 
Question
What exactly are the theoretical justifications for the different ways that we use MCMC draws to approximate the target distribution?
As an example, the theorems in Sections 4.5 and 4.7 of Geweke (2005) establish that the sample averages $(1/N)\sum_{i=1}^Nf(\theta_i)$ satisfy a central limit theorem with respect to the true value $\mathbb{E}[f(\theta)]$. Great. That takes care of (1).
What about (2) or (3)? What justifies using a kernel density estimator with MCMC draws, or computing the ECDF of MCMC draws?
 A: The main theoretical justification of (1), (2), and (3) is ergodicity of a Markov chain. Let $\mathcal{X}$ be the state space of the Markov chain and let $\mathcal{B}(\mathcal{X})$ be the Borel $\sigma$-algebra for $\mathcal{X}$. A Markov chain is characterized by its Markov transition kernel $P: \mathcal{X} \times \mathcal{B}(\mathcal{X}) \to [0,1]$. That is for element $x \in \mathcal{X}$ and set $A \in \mathcal{B}(\mathcal{X})$, 
$$P(x, A) = \Pr(X_{1} \in A \mid X_0 = x)\,.$$ 
After taking $t$ steps, the $t$-step transition for the Markov chain is 
$$P^{t}(x, A) = \Pr(X_{t} \in A \mid X_0 = x)$$
Under certain regularity conditions (aperiodicity, irreducibility and Harris recurrence), the $t$-step transition kernel converges to the stationary distribution $\pi$ in total-variation norm. That is,
$$ \|P^t(x, \cdot) - \pi(\cdot)\|_{TV} \to 0 \text{ as } t \to \infty\,.$$
Convergence in total variation is strong than convergence in distribution. Thus, for large $t$, drawing from $P^{t}(x, \cdot)$ is approximately drawing from $\pi$. As a consequence of this you get (1) by the Birkhoff Ergodic Theorem for Markov chains, (2) by the definition of convergence in distribution, and (3) because it is essentially an expectation so it is a special case of (1).
