The Quick Answer
The big yellow book (namely Di Francesco et al) that the OP quotes, largely obscures the distinction between what I call (b) and (c) below. If the OP is just interested in deriving the result in the fastest way he can Taylor-expand in $\delta$ the quantities $w(z+\delta/2)$, etc., and take the limit $\delta\rightarrow 0$. E.g.,
$$
w(z+\delta/2)\simeq w(z)+\frac{\delta}{2}\partial_zw(z)+\frac{1}{2!}\Big(\frac{\delta}{2}\Big)^2\partial_z^2w(z)+\frac{1}{3!}\Big(\frac{\delta}{2}\Big)^3\partial_z^3w(z)+\dots
$$
$$
\partial_zw(z+\delta/2)\simeq \partial_zw(z)+\frac{\delta}{2}\partial_z^2w(z)+\frac{1}{2!}\Big(\frac{\delta}{2}\Big)^2\partial_z^3w(z)+\dots
$$
The higher order terms do not contribute.
Using messy but straight-forward algebraic manipulations one indeed finds that equation (3) implies (4) in the OP's question. For example, since the denominator in (3) seems to have been causing the OP some trouble I'll also note that (from the above Taylor expansion it follows that):
$$
\big(w(z+\delta/2)-w(z-\delta/2)\big)^2=\big(\partial_zw(z)\big)^2\delta^2+\frac{1}{12}\big( \partial_z^3w\,\partial_zw(z)\big)\delta^4+\mathcal{O}(\delta^6),
$$
so the inverse is then,
$$
\frac{1}{[w(z+\delta/2)-w(z-\delta/2)]^2}=\frac{1}{\delta^2}\frac{1}{(\partial_zw(z))^2}-\frac{1}{12}\frac{\partial_z^3w(z)}{(\partial_zw(z))^3}+\mathcal{O}(\delta^2).
$$
Making use of the above Taylor expansions for the numerator, subtracting $1/\delta^2$ from the result and multiplying by a factor $1/2$ yields precisely (4).
However, I don't believe this is the right way to think about it: this derivation may be fast but it also hides many subtleties under the rug, so that in fact one has learnt very little (if not a negative amount) by following the big yellow book's derivation.
So I want to rather discuss a much more pedagogical (but also longer) derivation: in what follows we show how to derive an explicit expression for a normal ordered operator under any holomorphic change of coordinates in detail. (Hopefully, future readers interested in related questions will also benefit.) We take the energy-momentum tensor as our basic example.
The Long Answer
The OP asks to show that the energy-momentum tensor, $T(z)$, of a free scalar, $\phi(z)$, in 2 dimensions transforms with the funny Schwarzian derivative term under a conformal change of coordinates, $z\rightarrow w(z)$,
$$
\boxed{(\partial_{z_2}w_2)^2T^{'(w)}(w_2) = T^{(z)}(z_2)-\frac{1}{12}\Big[\frac{\partial_{z_2}^3w_2}{\partial_{z_2}w_2}-\frac{3}{2}\Big(\frac{\partial_{z_2}^2w_2}{\partial_{z_2}w_2}\Big)^2\Big]\,\,}
$$
taking as a starting point the defining equation for the normal ordered energy-momentum tensor for a free scalar which in my conventions reads:
$$
T^{(z)}(z_2) = \lim_{z_1\rightarrow z_2}-\frac{1}{2}\!:\!\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)\!:_z
$$
$$
T^{'(w)}(w_2) = \lim_{w_1\rightarrow w_2}-\frac{1}{2}\!:\!\partial_{w_1}\phi'(w_1)\partial_{w_2}\phi'(w_2)\!:_w
$$
I have deliberately cluttered the notation slightly (included some superscripts on $T$ and corresponding subscripts on $:\!(\dots)\!:$ and wrote $z_2,w_2$ rather than $z,w$ etc.) for reasons that will become clear momentarily. (In fact, this notation exposes the precise data on which these operators depend and will therefore allow us to track how these objects change as we change this data, one piece at a time. E.g., dropping the superscript from the energy-momentum tensor makes it impossible to distinguish between the quantities:
$$
T^{(z)}(z_2), \qquad\longleftrightarrow\qquad T^{(w)}(z_2),
$$
but this distinction will in turn play a crucial role below since it corresponds to changing normal ordering keeping fixed the coordinates - this is where the Schwarzian derivative makes its appearance. We can conversely also change coordinates keeping fixed the normal ordering,
$$
T^{(w)}(z_2), \qquad\longleftrightarrow\qquad T^{'(w)}(w_2),
$$
and this corresponds to the classical or ordinary change of coordinates (where one transforms coordinates assuming the object transforms as a conformal tensor) that is also used in the path integral. Incidentally, from an honest path integral viewpoint these points are manifest, and this is why people say the 'path integral is useful primarily because it provides the useful understanding/intuition', but I won't elaborate on this connection further. But let's go through the reasoning slowly and carefully.)
We will break down the computation into three independent steps:
(a) Normal ordering
(b) Change of normal ordering keeping coordinates fixed
(c) Change of coordinates keeping normal ordering fixed
We will proceed by exposing these three steps, (a), (b) and (c), (one at a time and in this order). Then, to change coordinates in any given normal ordered operator is to derive the map associated to the following composition:
$$
{\bf (c)}\circ{\bf (b)}\circ{\bf (a)}:\mathcal{O}(\phi)\,\longrightarrow \,\,??
$$
and when, e.g., $\mathcal{O}(\phi)$ is identified with the (non-normal-ordered) energy-momentum tensor then the "codomain" of this map will correspond to the coordinate-transformed normal-ordered energy-momentum tensor (given in terms of the Schwarzian derivative term above).
Let me add that the OP's question is a good question, in that I'm not even aware of a transparent and explicit derivation along these lines in the literature$^\dagger$ (but that doesn't mean it doesn't exist, somewhere..). The only paper that I know of that really exposes these issues is a paper by Polchinski (from 1987) on vertex operators, but there are intermediate steps between that paper and what follows that I'm not including here. Finally, I will focus on $c=1$ bulk scalars, $\phi(z)$, the generalisation to tensors (Grassmann-even or odd ghosts, matter fermions, etc.) and boundary operators being similar. For tensors the change of normal ordering with fixed coordinates then acquires an additional factor in the propagator but is otherwise entirely parallel.
$^\dagger$ UPDATE (Feb 2020): a derivation along the lines discussed below, including generalisation to the BRST framework, now exists (see Sec. 4).
We must first understand what it means to normal order an operator. We will use the path integral definition (although this is implicit).
(a) Normal Ordering
A normal ordering prescription is a prescription for subtracting infinities arising from self contractions within a (possibly composite) operator. In a free theory, such as the case of interest here, Wick's theorem gives all the self contractions and we therefore have succinctly that:$^*$
$$
\boxed{:\mathcal{O}(\phi)\!:_z \,\,= \mathcal{O}(\delta_J)\,\exp\Big(-\frac{1}{2}\int_{z'}\int_z\,J(z')J(z)\,G(z',z)+\int_z J(z)\phi(z)\Big)\Bigg|_{J=0}\,\,}
$$
where $G(z',z)=\langle\phi(z')\phi(z)\rangle$ is the free propagator used in the $z$ normal ordering, e.g., it will be sufficient to consider the standard expression for scalars:
$$
G(z',z) = -\ln |z'-z|^2,
$$
where (to justify the name '$z$ normal ordering') by $z,z'$ we implicitly mean here $z(p),z(p')$, where $p,p'$ are points on the surface, so $z$ is really a holomorphic chart coordinate.$^{**}$ Note also that I'm using a more traditional normalisation than the OP (obtained by taking $g=1/(4\pi)$). The integration measures, e.g. $d^2z$, are implicit in the above boxed expression (and we could more fully write $J(z,\bar{z})$ instead of $J(z)$, etc.).$^{***}$
The quantity $\mathcal{O}(\phi)$ is any (typically infinite if the elementary constituents are evaluated at coincident points) operator of interest, such as:
$$
\mathcal{O}(\phi) = \lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)\Big],
$$
where we will take $z_1,z_2$ to be coordinate points specified in the $z$ coordinate system, e.g., $z_1\equiv z(p_1)$, where $p_1$ is a marked point on the surface. Let us check that the boxed equation makes sense,
\begin{equation}
\begin{aligned}
T^{(z)}(z_2)&\equiv\,:\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)\Big]\!:_z \\
&=\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\frac{\delta}{\delta J(z_1)}\partial_{z_2}\frac{\delta}{\delta J(z_2)}\Big]\,\exp\Big(-\frac{1}{2}\int_{z'}\int_z\,J(z')J(z)\,G(z',z)+\int_z J(z)\phi(z)\Big)\Bigg|_{J=0}\\
&=\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\frac{\delta}{\delta J(z_1)}\partial_{z_2}\Big]\,\Big(-\int_{z'}\,J(z')\,G(z',z_2)+\phi(z_2)\Big)\\
&\qquad\qquad\times\exp\Big(-\frac{1}{2}\int_{z'}\int_z\,J(z')J(z)\,G(z',z)+\int_z J(z)\phi(z)\Big)\Bigg|_{J=0}\\
&=\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\partial_{z_2}\Big]\,\Big(\phi(z_1)\phi(z_2)-G(z_1,z_2)\Big)\\
&=\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\,\Big(\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)-\partial_{z_1}\partial_{z_2}G(z_1,z_2)\Big)\Big]\\
&=\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\,\Big(\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)+\frac{1}{z_{12}^2}\Big)\Big]\\
\end{aligned}
\end{equation}
where in the second equality we used the boxed equation above, in the third we carried out one of the two functional derivatives using the defining property,
$$
\int_z \frac{\delta J(z)}{\delta J(z_2)}f(z)=\int_z \delta^2(z-z_2)f(z)=f(z_2),
$$
in the fourth equality we carried out the remaining functional derivative and set $J=0$, and in the sixth we made use of the definition of $G(z',z)$ above (with $z_{12}\equiv z_1-z_2$).
So this defines what we mean by 'the energy momentum for a scalar in the $z$ normal ordering'.
Incidentally, inside the normal ordering we can freely take the limit as it is non-singular,
$$
:\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)\Big]\!:_z\,=\,:\!\Big[-\frac{1}{2}\partial_{z_2}\phi(z_2)\partial_{z_2}\phi(z_2)\Big]\!:_z\,.
$$
$^*$ If you are curious and happen not to know that this is equivalent to Joe Polchinski's conformal normal ordering definition, namely (2.2.7) in his volume 1 (or his vertex operator paper where he introduced it), the hint is on p.152 in Coleman's book 'Aspects of Symmetry'. (As a historical note, Joe once mentioned that he learnt all about normal ordering in 2-d quantum field theories from Coleman's lectures.)
$^{**}$ To avoid confusion let me be pedantic and mention that the integrals over $z,z'$ integrate over images of all points $p,p'$ in the manifold using $z$ chart coordinates, rather than integrating over all chart coordinates for fixed $p,p'$! (Had I not made the notation so explicit it would most likely not have been exposed how subtle but sharp all of these steps actually are; and there's more I'm not even mentioning for the sake of "brevity", otherwise this post would turn into a book..)
$^{***}$ The boxed equation above that defines normal ordering is actually a "baby version" of equation (3.1) in this paper; the latter provides the natural generalisation of the notion of normal ordering to interacting theories where it is termed complete normal ordering. For free theories (the case of interest here) the two notions are indistinguishable.
(b) Change of Normal Ordering (keeping coordinates fixed)
Very generally, we obtain different normal ordering prescriptions by replacing $G(z',z)$ in the above boxed equation by $G(z',z)+\Delta(z',z)$. We want to do something more specific here, namely we want to go through the exact same computation as we did above but in (what we will call) the '$w$ normal ordering'. We define the latter to be related to the $z$ normal ordering by a conformal transformation, $z\rightarrow w(z)$, by which we mean precisely the following: we are to simply$^{****}$ replace $G(z',z)$ by $G(w(z'),w(z))$ on the right-hand side in the above boxed equation keeping everything else fixed,
$$
\boxed{:\mathcal{O}(\phi)\!:_w \,\,= \mathcal{O}(\delta_J)\,\exp\Big(-\frac{1}{2}\int_{z'}\int_z\,J(z')J(z)\,G(w(z'),w(z))+\int_z J(z)\phi(z)\Big)\Bigg|_{J=0}\,\,}
$$
The subscript $w$ on the left-hand side is the reminder that this is $w$ normal ordering and the corresponding $w$ dependence on the right-hand side is entirely explicit (and contained solely in $G(w(z'),w(z))$). This is the definition of '$w$ normal ordering'. Notice that it is defined with respect to the reference/auxiliary '$z$ normal ordering'. (Clearly, we can similarly define a, say, '$u$ normal ordering' in precisely the same way, namely we simply replace $w$ by $u$, and that also will then be defined with respect to the reference '$z$ normal ordering', or we can consider $w(u(z))$ normal ordering, etc.., depending on context.)
Let us apply $w$ normal ordering to the case of interest,
\begin{equation}
\begin{aligned}
T^{(w)}(z_2)&\equiv\,:\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)\Big]\!:_w \\
&=\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\frac{\delta}{\delta J(z_1)}\partial_{z_2}\frac{\delta}{\delta J(z_2)}\Big]\,\exp\Big(-\frac{1}{2}\int_{z'}\int_z\,J(z')J(z)\,G(w(z'),w(z))+\int_z J(z)\phi(z)\Big)\Bigg|_{J=0}\\
&=\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\frac{\delta}{\delta J(z_1)}\partial_{z_2}\Big]\,\Big(-\int_{z'}\,J(z')\,G(w(z'),w(z_2))+\phi(z_2)\Big)\\
&\qquad\qquad\times\exp\Big(-\frac{1}{2}\int_{z'}\int_z\,J(z')J(z)\,G(w(z'),w(z))+\int_z J(z)\phi(z)\Big)\Bigg|_{J=0}\\
&=\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\partial_{z_2}\Big]\,\Big(\phi(z_1)\phi(z_2)-G(w(z_1),w(z_2))\Big)\\
&=\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\,\Big(\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)-\partial_{z_1}\partial_{z_2}G(w(z_1),w(z_2))\Big)\Big],
\end{aligned}
\end{equation}
the steps being identical to the above. We next consider the last term in detail. We are interested in the limit $z_1\rightarrow z_2$. Since $w(z_1)$ is by definition a holomorphic function of $z_1$ this means we can Taylor expand it around $z_2$ in $G(w(z_1),w(z_2))$,
\begin{equation}
\begin{aligned}
G(w(z_1),w(z_2))&=-\ln\big|w(z_1)-w(z_2)\big|^2\\
&=-\ln\Big|\sum_{n=0}^{\infty}\frac{1}{n!}z_{12}^n\partial_{z_2}^nw(z_2)-w(z_2)\Big|^2\\
&=-\ln\Big|\sum_{n=1}^{\infty}\frac{1}{n!}z_{12}^n\partial_{z_2}^nw(z_2)\Big|^2\\
&=-\ln\Big|z_{12}\sum_{n=1}^{\infty}\frac{1}{n!}z_{12}^{n-1}\partial_{z_2}^{n}w(z_2)\Big|^2\\
&=G(z_1,z_2)-\ln\Big|\sum_{n=1}^{\infty}\frac{1}{n!}z_{12}^{n-1}\partial_{z_2}^{n}w(z_2)\Big|^2\\
\end{aligned}
\end{equation}
Now I will leave the following as a fun
EXERCISE: Let us write $w_1\equiv w(z_1)$ and $w_2\equiv w(z_2)$. Show that for $|z_{12}|=|z_1-z_2|$ small:
$$
\partial_{z_1}\partial_{z_2}\ln\Big|\sum_{n=1}^{\infty}\frac{1}{n!}z_{12}^{n-1}\partial_{z_2}^{n}w(z_2)\Big|^2=\frac{2}{12}\bigg[\frac{\partial_{z_2}^3w_2}{\partial_{z_2}w_2}-\frac{3}{2}\Big(\frac{\partial_{z_2}^2w_2}{\partial_{z_2}w_2}\Big)^2\bigg]+\frac{1}{12}\bigg[3\Big(\frac{\partial_{z_2}^2w_2}{\partial_{z_2}w_2}\Big)^3+\frac{\partial_{z_2}^4w_2}{\partial_{z_2}w_2}-4\frac{\partial_{z_2}^3w_2\,\partial_{z_2}^2w_2}{(\partial_{z_2}w_2)^2}\bigg]\,z_{12}+\mathcal{O}(z_{12}^2).
$$
This follows directly by using the chain rule, taking into account that only the $z_{12}^{n-1}$ terms depend on $z_1$ and that both $z_{12}^{n-1}$ and $\partial_{z_2}^{n}w(z_2)$ depend on $z_2$. Since only the $z_{12}\rightarrow 0$ limit is of interest we can drop all terms on the right-hand side that vanish in this limit.
Substituting the result of this exercise into the above we learn that:
$$
\boxed{\lim_{z_1\rightarrow z_2}\partial_{z_1}\partial_{z_2}G(w(z_1),w(z_2))=\lim_{z_1\rightarrow z_2}\partial_{z_1}\partial_{z_2}G(z_1,z_2)-\frac{2}{12}\bigg[\frac{\partial_{z_2}^3w_2}{\partial_{z_2}w_2}-\frac{3}{2}\Big(\frac{\partial_{z_2}^2w_2}{\partial_{z_2}w_2}\Big)^2\bigg]}
$$
Let us in turn substitute this into the above expression for $T^{(w)}(z_2)$,
\begin{equation}
\begin{aligned}
T^{(w)}(z_2)&\equiv\,:\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)\Big]\!:_w \\
&=\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\,\Big(\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)-\partial_{z_1}\partial_{z_2}G(w(z_1),w(z_2))\Big)\Big]\\
&=\lim_{z_1\rightarrow z_2}\bigg\{-\frac{1}{2}\,\Big(\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)-\partial_{z_1}\partial_{z_2}G(z_1,z_2)+\frac{2}{12}\Big[\frac{\partial_{z_2}^3w_2}{\partial_{z_2}w_2}-\frac{3}{2}\Big(\frac{\partial_{z_2}^2w_2}{\partial_{z_2}w_2}\Big)^2\Big]\Big)\bigg\}\\
&=\lim_{z_1\rightarrow z_2}\bigg\{-\frac{1}{2}\,\Big(\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)-\partial_{z_1}\partial_{z_2}G(z_1,z_2)\Big)\bigg\}-\frac{1}{12}\Big[\frac{\partial_{z_2}^3w_2}{\partial_{z_2}w_2}-\frac{3}{2}\Big(\frac{\partial_{z_2}^2w_2}{\partial_{z_2}w_2}\Big)^2\Big]\\
&=:\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)\Big]\!:_z-\frac{1}{12}\Big[\frac{\partial_{z_2}^3w_2}{\partial_{z_2}w_2}-\frac{3}{2}\Big(\frac{\partial_{z_2}^2w_2}{\partial_{z_2}w_2}\Big)^2\Big]\\
&=T^{(z)}(z_2)-\frac{1}{12}\Big[\frac{\partial_{z_2}^3w_2}{\partial_{z_2}w_2}-\frac{3}{2}\Big(\frac{\partial_{z_2}^2w_2}{\partial_{z_2}w_2}\Big)^2\Big]\\
\end{aligned}
\end{equation}
where we noted in the last two lines that:
\begin{equation}
\begin{aligned}
T^{(z)}(z_2)&\equiv \,:\lim_{z_1\rightarrow z_2}\Big[-\frac{1}{2}\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)\Big]\!:_z\\
&=
\lim_{z_1\rightarrow z_2}\bigg\{-\frac{1}{2}\,\Big(\partial_{z_1}\phi(z_1)\partial_{z_2}\phi(z_2)-\partial_{z_1}\partial_{z_2}G(z_1,z_2)\Big)\bigg\}
\end{aligned}
\end{equation}
as shown above.
So we learn that a finite holomorphic change in normal ordering, $z\rightarrow w(z)$, with fixed coordinates, $z_2$, of the energy-momentum tensor is given by:
$$
\boxed{T^{(w)}(z_2)=T^{(z)}(z_2)-\frac{1}{12}\Big[\frac{\partial_{z_2}^3w_2}{\partial_{z_2}w_2}-\frac{3}{2}\Big(\frac{\partial_{z_2}^2w_2}{\partial_{z_2}w_2}\Big)^2\Big]}\qquad\qquad (*)
$$
Notice that we have not actually changed the coordinates to derive the Schwarzian derivative! Evidently, the entire content of the Schwarzian derivative lives entirely in the change of normal ordering of the energy-momentum tensor keeping the coordinates fixed.
$^{****}$ As mentioned above, this procedure is as simple as stated in the case of scalars; it is slightly more complicated for ghosts and matter fermions or tensors more generally.
(c) Change of Coordinates (keeping normal ordering fixed)
All that remains is to change coordinates, $z_2\rightarrow w_2\equiv w(z_2)$. Since $\phi(z_2)$ transforms as a scalar and its derivative as a weight-1 operator we have that,
$$
\partial_{z_2}\phi(z_2)dz_2 = \partial_{w_2}\phi'(w_2)dw_2.
$$
Furthermore, since we have treated the change in normal ordering separately from the change of coordinates we can now perform the change of coordinates just as we would do naively, and I want to emphasise the following statement (which follows from the defining equation of normal ordering above) as strongly as possible:
$T^{(w)}(z_2)$ does transform as a (weight-2) holomorphic tensor under a holomorphic change of coordinates provided we keep the normal ordering fixed:
$$
T^{(w)}(z_2)dz_2^2 = T^{'(w)}(w_2)dw_2^2\qquad\Rightarrow\qquad \boxed{T^{(w)}(z_2) = (\partial_{z_2}w_2)^2T^{'(w)}(w_2)}
$$
independently of the fact that the central charge of a free scalar $c=1$.
So you see why I insisted on using the cluttered notation above. Omitting the normal ordering symbol, $(w)$, from $T^{(w)}(z_2)$ clearly obscures the meaning of this local operator, while also leading to the perception that the Schwarzian derivative is somehow generated by a change of coordinates - as we have just seen, it is the change in normal ordering that is doing all of the magic. Unfortunately, almost all of the CFT literature (as do I most of the times) drops the normal ordering from the notation causing all sorts of unnecessary confusion.
Summarising
The final step of the derivation is trivial, we simply gather what we have learnt. In particular, we substitute the relation we have just derived into (*), leading precisely to the final answer for the transformation of the energy-momentum tensor under a change of coordinates (with a corresponding change of normal ordering),
$$
\boxed{(\partial_{z_2}w_2)^2T^{'(w)}(w_2) =T^{(z)}(z_2)-\frac{1}{12}\Big[\frac{\partial_{z_2}^3w_2}{\partial_{z_2}w_2}-\frac{3}{2}\Big(\frac{\partial_{z_2}^2w_2}{\partial_{z_2}w_2}\Big)^2\Big]}
$$
Note also that using the OPE with the energy momentum tensor generates the infinitesimal version of this last relation: it automatically merges steps (b) and (c) above.
UPDATE (Feb 2020): Further details along the lines of this derivation (for arbitrary operators and in the BRST framework) can be found in Sec. 4 in:
D. Lüst and D. Skliros, ``Handle Operators in String Theory,'' arXiv:1912.01055 [hep-th]
Best Answer
I think the original source of this claim is the famous unpublished paper of Luescher and Mack. Everyone's citing it. It is more rigorous mathematically and general (they don't assume parity) than Di Francesco. It starts on pages 1-2 of the manuscript. The proof below is basically the same proof, just with added details and a little bit different notation.
Assumptions: $T_{\mu\nu}$ are local covariant fields of dimension 2 and $$T^{\dagger}_{\mu\nu} = T_{\mu\nu}, \quad T_{\mu\nu} = T_{\nu\mu},\quad \partial^{\mu}T_{\mu\nu} = 0.$$ Moreover, one assumes that the correlation functions are well-defined and behaved, and satisfy usual axioms of CFT as stated e.g. in this page of nLab.
Proof: First we define the 2-point functions $$ S_{\mu\nu\rho\sigma}(\vec{x}, \vec{y}) := \text{analytic continuation of } \langle \Omega \mid T_{\mu\nu}(\vec{x}) T_{\rho\sigma}(\vec{y}) \mid \Omega \rangle. $$ All $S$'s are translationally invariant, so we may as well set $\vec{y}=0$ so that we are left with $$ S_{\mu\nu\rho\sigma}(\vec{x}) := \langle \Omega \mid T_{\mu\nu}(\vec{x}) T_{\rho\sigma}(0) \mid \Omega \rangle. $$ This function is real analytic for $\vec{x}\in\mathbb{R}^2,\, \vec{x}\neq 0$, and because of invariance under dilations, and the assumption that $T_{\mu\nu}$'s are of conformal weight 2, we have
\begin{equation} S_{\mu\nu\rho\sigma}(\vec{x})=\lambda^4 S_{\mu\nu\rho\sigma}(\lambda \vec{x})\quad\forall\lambda\in\mathbb{R}.\quad\quad (1) \end{equation} Moreover, from $T_{\mu\nu}=T_{\nu\mu}$ we get $$ S_{\mu\nu\rho\sigma}(\vec{x})=S_{\nu\mu\rho\sigma}(\vec{x})=S_{\mu\nu\sigma\rho}(\vec{x}).\quad\quad (2) $$ Furthermore, by locality, invariance under translations and Equation 1 with $\lambda=-1$ we get $$ S_{\mu\nu\rho\sigma}(\vec{x},0)=S_{\rho\sigma\mu\nu}(0,\vec{x})=S_{\rho\sigma\mu\nu}(-\vec{x},0)=S_{\rho\sigma\mu\nu}(\vec{x},0), $$ or in shorthand notation $$ S_{\mu\nu\rho\sigma}(\vec{x}) = S_{\rho\sigma\mu\nu}(\vec{x}) .\quad\quad (3) $$ From Equations 2 and 3 we see that $S_{\mu\nu\rho\sigma}(\vec{x})$ has 6 independent components and also taking into account Equation 1 we can write $$ S_{\mu\nu\rho\sigma}(\vec{x}) = (\vec{x}^{\,2})^{-4}\sum_{i=1}^6 \alpha_i f^{\,i} _{\mu\nu\sigma\rho}(\vec{x}),\quad \alpha_i\in\mathbb{C},\quad\quad (4) $$ where \begin{align} &f^{1}_{\mu\nu\rho\sigma}(\vec{x})=(\vec{x}^{\,2})^2 g_{\mu\nu} g_{\rho\sigma},\\\\ &f^{2}_{\mu\nu\rho\sigma}(\vec{x})=(\vec{x}^{\,2})^2 (g_{\mu\rho} g_{\nu\sigma}+g_{\mu\sigma}g_{\nu\rho}),\\\\ &f^{3}_{\mu\nu\rho\sigma}(\vec{x})=\vec{x}^{\,2} (g_{\mu\nu} x_\rho x_\sigma + g_{\rho\sigma} x_\mu x_\nu ),\\\\ &f^{4}_{\mu\nu\rho\sigma}(\vec{x})=x_\mu x_\nu x_\rho x_\sigma,\\\\ &f^{5}_{\mu\nu\rho\sigma}(\vec{x})= \vec{x}^{\,2} \left (g_{\mu\nu} (x_\rho \varepsilon_{\sigma\delta}x^\delta+x_\sigma\varepsilon_{\rho\delta}x^\delta )+ g_{\rho\sigma} (x_\mu\varepsilon_{\nu\delta} x^\delta + x_\nu \varepsilon_{\mu\delta}x^\delta) \right),\\\\ &f^{6}_{\mu\nu\rho\sigma}(\vec{x})= \left( x_\mu \varepsilon_{\nu\delta}x^\delta +x_\nu \varepsilon_{\mu\delta}x^\delta \right) x_\rho x_\sigma + x_\mu x_\nu \left( x_\rho \varepsilon_{\sigma\delta}x^\delta + x_\sigma \varepsilon_{\rho\delta}x^\delta \right), \end{align} and $$ g_{\mu\nu}=\delta_{\mu\nu},\quad \varepsilon_{\mu\nu}=-\varepsilon_{\nu\mu},\quad \varepsilon_{01}=+1,\quad \vec{x}^{\,2}=(x^0)^2+(x^1)^2. $$ The continuity equation $\partial^\mu T_{\mu\nu}=0$ implies that $\partial^\mu S_{\mu\nu\rho\sigma}(\vec{x}) = 0$ and thus reduces the number of independent $\alpha_i$'s to 2 arbitrary constants $\alpha_+$ and $\alpha_-$ (checked myself using Mathematica): \begin{alignat*}{6} &\alpha_1 = 3 \alpha_+,\quad &&\alpha_2 = -\alpha_+,\quad \alpha_3 = -4\alpha_+,\quad \alpha_4 = 8\alpha_+,&\\ &\alpha_5 = \alpha_-,\quad &&\alpha_6 = -2\alpha_-. & \end{alignat*} Inserting these values into Equation 4 we obtain: $$ S^{\mu}_{\;\mu\rho\sigma} (\vec{x}) = S^{\mu\;\rho}_{\;\mu\;\rho} (\vec{x}) = 0, $$ i.e. $\langle\Omega \mid T^\mu_{\;\mu}(\vec{x})\, T^\rho_{\;\rho} (0)\mid\Omega\rangle = 0$, which by the Reeh-Schlieder Theorem implies that $T^\mu_{\;\mu}(\vec{x})=0$.
Now to be able to apply Reeh-Schlieder, we actually have to go to $\vec{x}=0$ as you have correctly pointed out. And the reason why we can do this is analytic continuation, citing p 110 of R. JOST: The General Theory of Quantized Fields, 1965 (the same step in a related lemma):
"Such a Wigtman distribution thus vanishes in the real regularity points and thus, by analytic continuation identically".
tl;dr analytic continuation