Operator Product Expansion (OPE) in Conformal Field Theory We denote local operators of a conformal field theory (CFT) as $\mathcal{O}_i$ where $i$ runs over the set of all operators. Formally, the operator product expansion (OPE) is given by,
$$\mathcal{O}_i(z,\bar{z})\mathcal{O}_j(\omega,\bar{\omega}) = \sum_{k} C^k_{ij}(z-\omega,\bar{z}-\bar{\omega})\mathcal{O}_k(\omega,\bar{\omega})$$ 
where we have left the $\langle...\rangle$ implicit. My question is primarily, in practice, how does one determine the operator product expansion for an explicit case? Is the OPE similar to a Laurent expansion? In which case, can one determine the residues from the OPE? I would appreciate an explicit non-trivial example.

I'm following String Theory and M-Theory by Becker, Becker and Schwarz, and the OPE is dealt with somewhat rapidly.
 A: In the limit of the positions of two operator insertions approaching one another, the product of those two operators can be approximated as a series of local operators. This can be done to arbitrary accuracy and within conformal field theory, the expression is exact. 
The series described above will in general have singular behaviour. At first, this might sound bad, but it is actually the key point of the whole formalism. In fact, all the physical information we are interested in is encoded in the singularities. As a consequence, regular terms are usually not even written down. 
An explicit example would be the OPE of the conformal stress energy tensor with itself, given by
$$T(z)T(z')=\frac{\partial T(z')}{z-z'}+\frac{2T(z')}{(z-z')^2}+\frac{c}{2(z-z')^4}+\dots,$$
where $c$ is the central charge and dots denote regular terms. Furthermore, this expression is to be understood as an operator expression within a vacuum expectation value. This is, for convenience, often omitted in the literature. 
The terms I have written down above are singular in the limit $z\rightarrow z'$ and correspond to the residues of the OPE. As you have suggested correctly, the latter is formally equivalent to a Laurent series. 
But why is the OPE of interest, other than for providing a convenient way of writing products of operators? It turns out that it encodes the transformation behaviour of an operator under conformal transformations. This can be shown by working out the respective Ward identities. Furthermore, the OPE contains information about the mode algebra of the operators, i.e. the Virasoro algebra.
As proof of the above statements would exceed the boundaries of my answer, I refer at this point to the literature. I can recommend both the book by Polchinski and two sets of lecture notes freely available on the internet by Maximilian Kreuzer and David Tong. 
A: A hint about OPEs:  They exist for CFTs in any dimension.  They play nicely with the holomorphic phenomena in 2d, but they don't require it.  They encode information about the singularities which occur when you multiply operator-valued distributions. 
Explicit examples in higher dimension are computationally challenging.
A: For the case in which one of the fields on the LHS is holomorphic (this is typically the case if you compute the OPE of the stress-energy tensor with some field), we write the OPE $\, \, \ \mathcal{O}_i(z) \mathcal{O}_j(w, \bar{w}) \sim \sum_k C_{ij}^k ( z-w) \mathcal{O}_k(w, \bar{w}) \, \, \,$ if and only if

$ \partial_{\bar{z}} \ \langle \, \mathcal{O}_i(z) \mathcal{O}_j(w, \bar{w})  
- \sum_k C_{ij}^k ( z-w) \mathcal{O}_k(w, \bar{w}) \rangle = 0  $.

This means that 

$\mathcal{O}_i(z) \mathcal{O}_j(w, \bar{w})= \sum_k C_{ij}^k ( z-w) \mathcal{O}_k(w, \bar{w}) + X$

where $X$ is something holomorphic and regular, i.e such that $ \partial_{\bar{z}} \langle X \rangle= 0$. The "regular" requirement is necessary because non-regular holomorphic functions are not truly holomorphic. Namely:

$\partial_{\bar{z}} \frac{1}{z} = \pi \delta (\mathbf{x})$, 

where $\mathbf{x}$ is the point in Euclidean coordinates. \
(There is an analogous definition of the OPE for the case of an antiholomorphic field $\mathcal{O}_i(\bar{z})$.)
In order to compute OPEs we usually use two tools: Ward identities (which give you correlation functions associated to symmetries) and Wick theorem (which allows you to write correlation functions of several fields in terms of 2-point functions).
