Consistent estimator of conditional expectation, when conditioning on binary variable Suppose I have a sequence of i.i.d. random variables $\{Y_i,X_i,Z_i\}_{i=1}^n$ and $Z_i$ is binary.
Is the following a consistent estimator of $E(Y_i*X_i|Z_i=1)$ as $n\rightarrow \infty$? Under which theorem or law of large numbers?
$$
\frac{\sum_{i=1}^n Y_i*X_i*Z_i}{\sum_{i=1}^n Z_i}
$$

Now, let $W$ be a random variable and suppose that I have a sequence of random variables $\{Y_i\}_{i=1}^n$ that are i.i.d.  conditional on $W$. Consider some functions $\{f_i(W)\}_{i=1}^n$ and $\{g_i(W)\}_{i=1}^n$, where $g_i(\cdot)$ is binary for $i=1,...,n$. Is
$$
\frac{\sum_{i=1}^n Y_i*f_i(W)*g_i(W)}{\sum_{i=1}^n g_i(W)}
$$
a consistent estimator of $E(Y_i*f_i(W)| g_i(W))$? Under which theorem or law of large numbers?
 A: It is useful to generalise your problem to deal with the conditional expectation of an arbitrary function of a random vector.  To this end, suppose you have an IID sequence $\{ \mathbf{X}_i, Z_i | i \in \mathbb{N} \}$ where each $Z_i$ is a binary variable (taking values zero or one) and $\mathbf{X}_i = (X_{i,1},...,X_{i,m})$ is a random vector with some distribution.  Given some continuous measureable function $g:\mathbb{R}^m \rightarrow \mathbb{R}$ we will denote the following expectations:
$$\Gamma \equiv \mathbb{E}( Z_i \cdot g(\mathbf{X}_i))
\quad \quad \quad \quad \quad
\gamma \equiv \mathbb{P}( Z_i=1 ),$$
and note that their ratio gives the conditional expectation:
$$\Pi \equiv \frac{\Gamma}{\gamma} = \mathbb{E}( g(\mathbf{X}_i) | Z_i=1).$$
We also define the following estimators with $n$ data points:
$$\hat{\Gamma}_n \equiv \frac{1}{n} \sum_{i=1}^n Z_i \cdot g(\mathbf{X}_i)
\quad \quad \quad \quad \quad
\hat{\gamma}_n \equiv \frac{1}{n} \sum_{i=1}^n Z_i,$$
and note that their ratio gives the estimator:$^\dagger$
$$\hat{\Pi}_n \equiv \frac{\hat{\Gamma}_n}{\hat{\gamma}_n} = \frac{\sum_{i=1}^n Z_i \cdot g(\mathbf{X}_i)}{\sum_{i=1}^n Z_i}.$$
We will assume that $\gamma>0$ so that the conditioning event occurs with some non-zero probability.  (In the contrary case where the conditioning event occurs with probability zero, the conditional probability of interest can be set to any value, so any estimator that converges to a value is arguably a consistent estimator.)  From the law of large numbers (LLN) for IID random variables/vectors, we have $\hat{\Gamma}_n \rightarrow \Gamma$ and $\hat{\gamma}_n \rightarrow \gamma$.  Application of the continuous mapping theorem (CMT) then gives us the desired result:
$$\hat{\Pi}_n = \frac{\hat{\Gamma}_n}{\hat{\gamma}_n} \rightarrow \frac{\Gamma}{\gamma} = \Pi = \mathbb{E}(g(\mathbf{X}_i)|Z_i=1).$$
Since the LLN and the cCMT both hold for strong convergence, the estimator $\hat{\Pi}_n$ is strongly consistent for $\Pi$ (weak consistency also follows).

Application to your problem: Your question is a special case of this general result, where you are using a bivariate random vector $\mathbf{X}_i$ (which you have denoted as $(X_i,Y_i)$) and you use the function $g(x,y) = x \cdot y$ to take the product of its elements.  Application of the above result to this special case gives:
$$\hat{\Pi}_n = \frac{\sum_{i=1}^n X_i \cdot Y_i \cdot Z_i}{\sum_{i=1}^n Z_i}
\rightarrow \mathbb{E}(X_i \cdot Y_i|Z_i=1).$$
This confirms the asserted consistency of your estimator.  Extension to the second case where you take a function of your variables still gives the same result, so long as the functions you are using are continuous.  The relevant theorems are the law of large numbers and the continuous mapping theorem.  Depending on whether you wish to assert weak or strong consistency, the relevant theorems are the particular versions of the LLN and CMT for those cases.

$^\dagger$ In the case where $\hat{\gamma}_n = 0$ we will set $\hat{\Pi}_n \equiv 0$ by convention.  It can be set to other values without affecting the result of our analysis.
