(Ito lemma proof): convergence of $\sum_{i=0}^{n-1}f(W(t_{i}))(W(t_{i+1})-W(t_{i}))^{2}.$ The purpose of this question is to complete my personal exposition on the rigorous proof of Ito's lemma.  I have consulted more than half a dozen mathematical finance texts and not a single one, for all its importance, rigorously proves this.  The closest I've seen is Shreve's treatment in Stochastic Calculus for Finance Vol II

Let $\{W(t)\}_{t\geq0}$ be a standard Brownian motion.  Let $f$ be a smooth function, $T>0$, and $\Pi=\{t_{0}=0,t_{1},\ldots,t_{n}=T\}$ be a partition of $[0,T]$.
Consider the random variable
$$X(\omega)=\lim_{||\Pi||\to0}\;\sum_{i=0}^{n-1}f(W(t_{i}))(W(t_{i+1})-W(t_{i}))^{2}.$$
If $f\equiv1$, then almost surely $X\equiv T$.  This is the standard result that Brownian motion accumulates quadratic variation at a unit rate.
If we make the substitution $(W(t_{i+1})-W(t_{i}))^{2}\approx(t_{i+1}-t_{i})$, then after $||\Pi||\to0$ and summing, the errors committed with this approximation go to $0$ and we recover almost surely
$$X\equiv \int_{0}^{T}f(W(s))\;ds.$$
Is there a way to make this more precise?  I would like to adapt the usual proof that $$\lim_{||\Pi||\to0}\;\sum_{i=0}^{n-1}(W(t_{i+1})-W(t_{i}))^{2}=T,$$
which proceeds by showing that the expectation of the sampled quadratic variation (i.e. the sum on any partition) is $T$ and that the variance tends to $0$ as the partition becomes finer.  Indeed, $$\mathbb{E}[(W(t_{i+1})-W(t_{i}))^{2}]=\text{Var}[W(t_{i+1})-W(t_{i})]=t_{i+1}-t_{i}$$ and $$\text{Var}[(W(t_{i+1})-W(t_{i}))^{2}]=2(t_{i+1}-t_{i})^{2}$$ and the claim follows after summing, taking limits, and performing a simple estimate on the variance.  Based on my efforts, it does not seem like such simple formulas hold for the expectation and variance when the quantity is multiplied by $f(W(t_{i})).$
A heuristic for the approximation above can be obtained by noting that (if $\Pi$ is a uniformly spaced partition)
$$(W(t_{i+1})-W(t_{i}))^{2}=\frac{T}{n}\left(\frac{W(t_{i+1})-W(t_{i})}{\sqrt{t_{i+1}-t_{i}}}\right)^{2}.$$
Thus if
$$Y_{i}:=\left(\frac{W(t_{i+1})-W(t_{i})}{\sqrt{t_{i+1}-t_{i}}}\right)^{2},$$
then (being a collection of squared normal iid random variables with mean $1$) the law of large numbers implies
$$\sum_{i=0}^{n-1}(W(t_{i+1})-W(t_{i}))^{2}=T\sum_{i=0}^{n-1}\frac{Y_{i+1}}{n}\to T\;\text{as}\;n\to\infty.$$
Thus, while $(W(t_{i+1})-W(t_{i}))^{2}\neq t_{i+1}-t_{i}$, the law of large numbers says that the error committed with this approximation cancels out in the limit.  In other words, we may write
$$(W(t_{i+1})-W(t_{i}))^{2}=(t_{i+1}-t_{i})+\phi(t_{i})$$
where $\sum_{i=0}^{n-1}\phi(t_{i}))\to0\;\text{as}\;n\to\infty.$  What I'm trying to rigorously establish is that the above equations/limits all remain valid if we multiply them by $f(W(t_{i}))$.
In other words, is it correct to write $$f(W(t_{i}))(W(t_{i+1})-W(t_{i}))^{2}=f(W(t_{i}))(t_{i+1}-t_{i})+\phi(t_{i}))$$ where $\sum_{i=0}^{n}f(W(t_{i}))\phi(t_{i})\to0\;\text{as}\;n\to\infty?$
The obvious problem with the above heuristic is that the random variables
$$Y_{i}':=f(W(t_{i}))Y_{i}:=f(W(t_{i}))\left(\frac{W(t_{i+1})-W(t_{i})}{\sqrt{t_{i+1}-t_{i}}}\right)^{2}$$
are no longer independent, so the law of large numbers cannot be used in the manner above.
 A: Yes, there is a way to make this precise; in fact, it holds even in a more general framework. The following proof is adapted from René L. Schilling/Lothar Partzsch: Brownian Motion - An Introduction to Stochastic Processes, Lemma 17.4 (2nd edition).

Let $f \in C_b$ be a bounded continuous function and $(W_t)_{t \geq 0}$ a (one-dimensional) Brownian motion. Then $$\sum_{j=1}^n f(W_{t_{j-1}}) \cdot (W_{t_{j}}-W_{t_{j-1}})^2 \stackrel{L^2(\mathbb{P})}{\to} \int_0^T f(W_s) \, ds \quad \text{as $|\Pi| \to 0$}$$ where $|\Pi| := \max_j |t_j-t_{j-1}|$ is the mesh size of the partition $\Pi = \{0=t_0<\ldots<t_n=T\}$.

Proof: Since $t \mapsto f(W_t)$ is (almost surely) continuous, we know that $$\sum_{j=1}^n f(W_{t_{j-1}}) \cdot (t_j-t_{j-1}) \to \int_0^T f(W_s) \, ds$$ as $|\Pi| \to 0$. Therefore, it suffices to show
$$I_n := \sum_{j=1}^n f(W_{t_{j-1}}) \bigg[ (W_{t_j}-W_{t_{j-1}})^2 - (t_j-t_{j-1}) \bigg] \to 0 \tag{1}$$
as $|\Pi| \to 0$. We write $\Delta t_j := t_j-t_{j-1}$ and $\Delta W_j := W_{t_j}-W_{t_{j-1}}$. Recall  that $(W_t^2-t)_{t \geq 0}$ is a martingale with respect to the canonical filtration $(\mathcal{F}_t)_{t \geq 0}$, and therefore it is not difficult to see that
$$\begin{align*} &\quad \mathbb{E} \bigg( f(W_{t_{j-1}}) f(W_{t_{i-1}}) \cdot [\Delta W_j^2 - \Delta_j] \cdot [\Delta W_i^2-\Delta_i]\bigg)\\ &= \mathbb{E} \bigg( \mathbb{E} \bigg( f(W_{t_{j-1}}) f(W_{t_{i-1}}) \cdot [\Delta W_j^2 - \Delta_j] \cdot [\Delta W_i^2-\Delta_i] \mid \mathcal{F}_{t_{j-1}} \bigg) \bigg) \\ &= \mathbb{E} \bigg( f(W_{t_{j-1}}) f(W_{t_{i-1}}) \cdot [\Delta W_i^2-\Delta_i] \cdot \underbrace{\mathbb{E} \bigg( \Delta W_j^2 - \Delta_j \mid \mathcal{F}_{t_{j-1}} \bigg)}_{\mathbb{E}(\Delta W_j^2-\Delta j)=0} \bigg) = 0. \tag{2} \end{align*}$$
From the definition of $I_n$ and $(2)$ we see that
$$\mathbb{E}(I_n^2) = \mathbb{E}\left(\sum_{j=1}^n f(W_{t_{j-1}})^2 (\Delta W_j^2-\Delta_j)^2 \right).$$
(Note that the mixed terms vanish!) Using that $f$ is bounded and $W_t-W_s \sim W_{t-s} \sim \sqrt{t-s} W_1$ we find
$$\begin{align*} \mathbb{E}(I_n^2) &\leq \|f\|_{\infty}^2 \sum_{j=1}^n \mathbb{E}\bigg[(\Delta W_j^2-\Delta_j)^2\bigg] \\ &= \|f\|_{\infty}^2 \sum_{j=1}^n \Delta_j^2 \cdot \mathbb{E}[(W_1^2-1)^2] \\ &\leq C |\Pi| \sum_{j=1}^n \Delta_j = C |\Pi| T \end{align*}$$
for $C := \|f\|_{\infty}^2 \mathbb{E}[(W_1^2-1)^2]$. Letting $|\Pi| \to 0$, the claim follows.
