Let $M,N$ be continuous local martingales, and let $Y$ be a bounded continuous adapted process. Let $t\in \mathbb{R^+},$ and let $(\pi_{t}^n)_{n\in\mathbb{N}}$ be a sequence of partition of $[0,t]$ with $\lim\limits_{n\rightarrow \infty}\delta\pi_t^n=0$ where $\delta$ is the mesh of $\pi$. Let $$ Z_n=\sum\limits_{i=0}^{n-1}Y_{t_i}(M_{t_{i+1}}-M_{t_i})(N_{t_{i+1}}-N_{t_i}). $$ Prove that $\left\{Z_n\right\}_{n\in\mathbb{N}}$ converges to $\int_{0}^t Y_s\,d[M,N]_s$ in probability.
My idea is using the fact $ \sum\limits_{i=0}^{n-1}(M_{t_{i+1}}-M_{t_i})(N_{t_{i+1}}-N_{t_i}) $ converges to $[M,N]_t$ in probability.
Given $\alpha>0$, we consider $$\begin{multline} P\Biggl(\left\lvert Z_n-\int_{0}^t Y_s\,d[M,N]_s\right\rvert >\alpha\Biggr) \\ \begin{aligned}&= P\Biggl(\left\lvert\sum\limits_{i=0}^{n-1}Y_{t_i}(M_{t_{i+1}}-M_{t_i})(N_{t_{i+1}}-N_{t_i})-\sum\limits_{i=0}^{n-1}\int_{t_i}^{t_{i+1}}Y_s\,d[M,N]_s\right\rvert>\alpha\Biggr) \\ &=P\Biggl(\left\lvert\sum\limits_{i=0}^{n-1}\left\{Y_{t_i}(M_{t_{i+1}}-M_{t_i})(N_{t_{i+1}}-N_{t_i})-\int_{t_i}^{t_{i+1}}Y_s\,d[M,N]_s\right\}\right\rvert>\alpha\Biggr). \end{aligned}\end{multline}$$ If we can move $Y_{t_i},Y_s$ in front of the sum techniquely, then since $Y$ is bounded, and let $K$ be an upper bound of $Y$. Thus, we have $$\begin{multline} P\Biggl(\left\lvert Z_n-\int_{0}^t Y_s\,d[M,N]_s\right\rvert>\alpha\Biggr) \\ \begin{aligned}&\leq P\Biggl(\left\lvert K\right\rvert\left\lvert\sum\limits_{i=0}^{n-1}\left\{(M_{t_{i+1}}-M_{t_i})(N_{t_{i+1}}-N_{t_i})-\int_{t_i}^{t_{i+1}}d[M,N]_s\right\}\right\rvert>\alpha\Biggr) \\ &=P\Biggl(\left\lvert\sum\limits_{i=0}^{n-1}\left\{(M_{t_{i+1}}-M_{t_i})(N_{t_{i+1}}-N_{t_i})-\int_{t_i}^{t_{i+1}}d[M,N]_s\right\}\right\rvert>\frac{\alpha}{K}\Biggr) \\ &=P\Biggl(\left\lvert\sum\limits_{i=0}^{n-1}(M_{t_{i+1}}-M_{t_i})(N_{t_{i+1}}-N_{t_i})-[M,N]_t\right\rvert>\frac{\alpha}{K}\Biggr). \end{aligned}\end{multline}$$ Then, the last probability converges to $0$ by the fact $ \sum\limits_{i=0}^{n-1}(M_{t_{i+1}}-M_{t_i})(N_{t_{i+1}}-N_{t_i}) $ converges to $[M,N]_t$ in probability.
Thus, my problem is how to deal with $Y_{t_j}$ and $Y_s$? Is there any good way to move it out of the sum? or, any other suggestion to solve this problem. Thanks.
For simplicity of notation I'll consider the case $M=N$; the proof can be easily generalized to the general setting (alternatively, you can use the polarization identity).
Denote by $(\tau_k)_{k \in \mathbb{N}}$ a localizing sequence of $(M_t)_{t \geq 0}$. Since $(M_t)_{t \geq 0}$ has continuous sample paths, we may assume without loss of generality that $|M_{t \wedge \tau_k}| \leq k$ for all $t \geq 0$, $k \in \mathbb{N}$ (otherwise replace $\tau_k$ by $\tau_k \wedge \inf\{t; |M_t| \geq k\}$)). Writing
$$\int_{t_i}^{t_{i+1}} Y_s \, d[M]_s = \sum_{i=0}^{n-1} Y_{t_i} ([M]_{t_{i+1}}-[M]_{t_i}) + \int_{t_i}^{t_{i+1}} (Y_s-Y_{t_i}) \, d[M]_s$$
and following the reasoning in the first part of your answer, we find
$$\begin{align*} &\quad \mathbb{P} \left( \left| Z_n- \int_0^t Y_s \, d[M]_s \right|> \alpha \right) \\ &\leq \mathbb{P} \left( \left| \sum_{i=0}^{n-1} Y_{t_i} \big( (M_{t_{i+1}}-M_{t_i})^2 - ([M]_{t_{i+1}}-[M]_{t_i}) \big) \right| > \frac{\alpha}{2} \right) \\ &\quad + \mathbb{P} \left( \left| \sum_{i=0}^{n-1} \int_{t_i}^{t_{i+1}} (Y_s-Y_{t_i}) \, d[M]_s \right|> \frac{\alpha}{2} \right) \\ &=: I_1 + I_2. \end{align*}$$
Clearly,
$$I_2 \leq \mathbb{P} \left( [M]_t \sup_{|s-r| \leq \delta \pi^n} |Y_r-Y_s| > \frac{\alpha}{2} \right);$$
since the mesh size $\delta \pi^n$ tends to $0$ and $t \mapsto Y_t(\omega)$ is uniformly continuous on compact sets, this shows that $I_2$ converges to $0$ as $n \to \infty$. It remains to estimate the first term; this is a little bit more tricky. We have
$$\begin{align*} I_1 &\leq \mathbb{P}(\tau_k \leq t) + \mathbb{P} \left( \left| \sum_{i=0}^{n-1} Y_{t_i} \big( (M_{t_{i+1} \wedge \tau_k}-M_{t_i \wedge \tau_k})^2 - ([M]_{t_{i+1} \wedge \tau_k}-[M]_{t_i \wedge \tau_k}) \big) \right| > \frac{\alpha}{2} \right)\\ & =: I_{11} + I_{12}. \tag{1} \end{align*}$$
(Recall that $(\tau_k)_k$ is the localizing sequence.) By Markov's inequality,
$$I_{12} \leq \frac{4}{\alpha^2} \mathbb{E} \left| \sum_{i=0}^{n-1} Y_{t_i} \big( (M_{t_{i+1} \wedge \tau_k}-M_{t_i \wedge \tau_k})^2 - ([M]_{t_{i+1} \wedge \tau_k}-[M]_{t_i \wedge \tau_k}) \big) \right|^2.$$
Using that $(M_{t \wedge \tau_k}^2-[M]_{t \wedge \tau_k})_{t \geq 0}$ is a martingale and the tower property of conditional expectation (see the Lemma below), it follows that
$$\begin{align*} I_{12} &\stackrel{\text{Lemma}}{\leq} \frac{4}{\alpha^2} \mathbb{E}\sum_{i=0}^{n-1} \left| Y_{t_i} \big( (M_{t_{i+1} \wedge \tau_k}-M_{t_i \wedge \tau_k})^2 - ([M]_{t_{i+1} \wedge \tau_k}-[M]_{t_i \wedge \tau_k}) \big) \right|^2 \\ &\leq \frac{4C^2}{\alpha^2} \mathbb{E}\sum_{i=0}^{n-1} \left| \big( (M_{t_{i+1} \wedge \tau_k}-M_{t_i \wedge \tau_k})^2 - ([M]_{t_{i+1} \wedge \tau_k}-[M]_{t_i \wedge \tau_k}) \big) \right|^2 \\ &\stackrel{\text{Lemma}}{=} \frac{4 C^2}{\alpha^2} \mathbb{E} \left| \sum_{i=0}^{n-1} \big( (M_{t_{i+1} \wedge \tau_k}-M_{t_i \wedge \tau_k})^2 - ([M]_{t_{i+1} \wedge \tau_k}-[M]_{t_i \wedge \tau_k}) \big) \right|^2 \\ &= \frac{4C^2}{\alpha^2} \mathbb{E} \left| \sum_{i=0}^{n-1} (M_{t_{i+1} \wedge \tau_k}-M_{t_i \wedge \tau_k})^2 - [M]_{t \wedge \tau_k} \right|^2. \end{align*}$$
Since $(M_{t \wedge \tau_k})_{t \geq 0}$ is a bounded martingale, the expression on the right-hand side converges to $0$ as $n \to \infty$. Letting first $n \to \infty$ and then $k \to \infty$ in $(1)$, we find that $I_2$ converges to $0$ as $n \to \infty$. This finishes the proof.
To prove this, show first that
$$\mathbb{E} \big( (N_t-N_s)^2 - ([N]_t-[N]_s) \mid \mathcal{F}_s \big)=0$$
for all $s \leq t$. Use this to show that in the squared sum
$$\left( \sum_{i=0}^n \dots \right)^2 = \sum_{i} \sum_j \dots$$
all terms with $i \neq j$ vanish.