Mathematical motivation of OPE? In Peskin & Schroeder (and also Cheng which I have skimmed through) they motivate the Operator Product Expansion with a lot of words. 
Is there any way to motivate it mathematically, e.g. Taylor expansion or similar?
 A: Firstly, I'd like to recommend Conformal Field Theory by Di-Francesco, it is a comprehensive text which is thorough and contains many applications of conformal field theory. The text is indispensable.

In conformal field theory, it is often characteristic of correlation functions to diverge as points of two or more fields coincide. The operator product expansion is essentially a Laurent series, and it represents a series of operators (always understood to be in correlation functions) as the sum of well-defined operators, times a factor which diverges as points coincide.  Recall under a conformal map $z\to w(z), \bar z \to \bar w(\bar z)$, a quasi-primary field transforms as,
$$\phi'(w,\bar w) = \left( \frac{\partial w}{\partial z} \right)^{-h}\left( \frac{\partial \bar w}{\partial \bar z} \right)^{-\bar h} \phi(z,\bar z)$$
where $(h,\bar h)$ are the conformal dimensions. The scaling dimension of the field is $\Delta = h+\bar h$ and planar spin $s = h - \bar h$. An example of an operator product expansion would be,
$$T(z) \phi(w,\bar w) \sim \frac{h}{(z-w)^2} \phi(w,\bar w) + \frac{1}{z-w} \partial_w \phi(w,\bar w) + \dots$$
with stress-energy tensor $T(z)$ where it is common practice to omit terms not singular as $z\to w$. This OPE in fact defines a primary operator. The OPE can also tell us other things. Consider a free boson,
$$\mathcal{S} = \frac{1}{2}g\int d^2x \, \partial_\mu \varphi \partial^\mu \varphi$$
wherein the propagator is given by,
$$\langle \varphi(x) \varphi(y)\rangle = -\frac{1}{4\pi g} \ln(x-y)^2$$
which in complex coordinates is,
$$\langle \varphi(z,\bar z)\varphi(w,\bar w)\rangle =  -\frac{1}{4\pi g} \left[ \ln(z-w) + \ln(\bar z - \bar w)\right]$$
By differentiating we can separate the holomorphic and anti-holomorphic parts:
$$\langle \partial_z \varphi(z,\bar z) \partial_w \varphi(w,\bar w) \rangle = -\frac{1}{4\pi g} \frac{1}{(z-w)^2}$$
and similarly with $(z \leftrightarrow \bar z)$, etc. The OPE of the field with itself is then,
$$\partial \varphi(z) \partial \varphi(w) \sim -\frac{1}{4\pi g} \frac{1}{(z-w)^2}$$
Notice as $z\to w$, i.e. if we exchange fields, the OPE does not change sign; this reflects the bosonic character of $\varphi$. A similar calculation for a fermion field reveals the OPE with itself is,
$$\psi(z) \psi(w) \sim \frac{1}{2\pi g}\frac{1}{z-w}$$
Exchanging $z$ and $w$ picks up a minus sign which reflects the anti-commuting or fermionic character of the field $\psi$.

As aforementioned, the OPE can be thought of as a Laurent series, which means we can compute residues from it. I won't go through the whole derivation, but the homomorphic component $J_z$ of the conserved current $J_a$ under conformal transformations has an OPE,
$$J_z(z)\mathcal{O}(w,\bar w)=  \dots + \frac{\mathrm{Res}[J_z(z)\mathcal{O}(w,\bar w)]}{z-w} + \dots$$
for a generic operator $\mathcal{O}$. It's analogous to when one tries to compute residues by computing the Laurent expansion of a function, rather than through the tedious limit formula. For example, for $\sin z / z$, we have,
$$\frac{\sin z}{z} = \frac{1}{z} \left( z-\frac{1}{3!} z^3 + \dots\right) = 1 -\frac{1}{3!} z^2 + \dots$$
The $z=0$ pole is a simple pole, and from the Laurent expansion, we can conclude the residue is zero. A second example: 
$$\frac{e^z}{\sin z} = \frac{1}{z} + 1 + \frac{2z}{3} + \frac{z^2}{3} + \dots$$
The pole $z=0$ is likewise simple, and this time the residue is one.
