Rigorously why there should be an operator product expansion in conformal field theory?
If we look at ordinary holomorphic functions, then the OPE is nothing but the Laurent series of the product of the two functions.
If you expect that the OPE is to be formalized as a rigorous convergence of operators in some space or something like that, then I have to disappoint you: The formalized version of operator product expansions are vertex operator algebras, in which the variables manipulated have entirely formal character and all the power series are formal.
The latter works because the OPEs follow from the structure of the Virasoro algebra and the transformation of the operators in the OPE under the Virasoro generators. When we view the OPE as a statement about holomorphic functions, this must be proven through use of the residue theorem. When we view the OPE as a statement in a vertex operator algebra, this way of computing the coefficients in the OPE through the behaviour under the Virasoro generators becomes the definition of the OPE.
In this way the construction of the OPE becomes an algebraic procedure that tells you how to associate a vertex operator algebra to a given "algebra of modes" such as the Virasoro algebra. It is then a rigorous algebraic statement that representations of the Virasoro algebra induce representations of its associated vertex operator algebra and therefore it is meaningful to use this vertex operator algebra in conformal field theory, since the space of state carries a representation of the Virasoro algebra and therefore also of this vertex operator algebra.
If you do not phrase your axiomatic CFT in the language of these algebras, then the existence of the OPE is an axiom itself, i.e. it is assumed, not shown, to exist.
For an introduction to a rigorous look at conformal field theories treating the Virasoro algebra as the fundamental starting point, see Schottenloher's book "A Mathematical Introduction To Conformal Field Theory".
Disclaimer: I am not giving an answer which is as mathematically rigorous as @ACuriousMind's but I believe it satisfies a physicist's definition of rigour. Moreover, it also is generalizable to higher dimensions.
Think of Operator Product Expansion as a corollary of State-Operator correspondence in any CFT in arbitrary d-dimension. In general, if you think of (backward-)evolution under dilatation operator, then you can evolve any arbitrary state defined on a unit (d-1)-sphere by 'shrinking' it on smaller and smaller spheres till you have a state defined on a sphere of small radius, $\epsilon\to0$. Now this state is local enough to be called an 'operator', or a linear sum of operators, $\mathcal O_k$s. Notice that there is no obstruction to this procedure as long as one doesn't encounter any other operator insertion while radially evolving towards the origin. TBH, this is where I loose my mathematical rigor because I am not being precise about the circumstances under which I can write a localized state as an operator but for most physicists I know this picture works. So once you are ready to accept this shrunk state as a linear sum over operators, think of a situation where you start with an initial state defined on (d-1)-sphere (of unit radius, say) by inserting two operators, $\mathcal O_i$ and $\mathcal O_j$, at arbitrary locations inside the sphere, say at $\vec{x}_1$ and $\vec{x}_2$ $\left(|\vec{x}_1|,|\vec{x}_1|<1\right)$. Now if you think about applying the radial evolution procedure I described above to this state then you can write it as sum of operator insertions at origin and this statement is `exact' as long as you include all the operators that need to be included. This procedure is accurate as long as you don't have a third operator inserted within the unit sphere. This whole process can be made slightly more rigorous and this rigor gives rise to the OPE coefficients $C^k_{ij}$ as has been shown in the works of Dolan and Osborn. So now you have written an OPE expansion as follows, $$\mathcal O_i(\vec x_1) \mathcal O_j(\vec x_2) = \sum_k \left.C^k_{ij} (\vec x_1, \vec x_2,\vec x,\vec\partial_x)\mathcal O_k(\vec x)\right|_{\vec x=0}~.$$ So you don't need necessarily write an OPE around one of the original insertions $\vec x_1, \vec x_2$, and not even origin $\vec x=0$ as I have described above. You can use translation and rotation to move around and write the OPE around arbitrary point as long as you don't have any obstruction due to any other operator insertions.
Edit: see David Simmons-Duffin lecture note's section 6.1-6.3 for the kind of discussion I have presented.