Operator product expansion (of free field theory) in functional quantization approach (Note: Those who are already familiar with the notion of OPE's might prefer to skip directly to the Question section below) 
Setup:
Consider, for simplicity, the quantum theory of the Euclidean K-G field.  That is, the "quantization" of the the classical field which obeys the equations of motion, 
\begin{align}
(\partial^2-m^2)\phi=0,
\end{align}
where $\partial^2\equiv \sum^4_{i=1} \partial^2/\partial (x^i)^2$ with $x^i$ denoting elements of 4-dimensional Cartesian coordinates.
Following the "canonical" quantization procedure for the free field a la Peskin and Schroeder, I find the 2-point (Schwinger?) function
\begin{align}
\langle 0| \phi(x_1)\phi(x_2)|0\rangle=\Delta_E(x_1-x_2), \tag{1} \label{2point}
\end{align}
where $\Delta_E$ denotes the Euclidean propagator/Green's function,
\begin{align}
\Delta_E (x_1-x_2)= \int \frac{d^4p}{(2\pi)^4}\frac{e^{i p\cdot (x_1-x_2)}}{p^2+m^2},
\end{align}
where $p,x_1,x_2$ are 4-vectors and "$\cdot$" denotes the usual Euclidean inner product.  
Now, the 2-point function \eqref{2point} clearly diverges in the coincidence limit $x_1\to x_2$, so it is customary to define the quadratic field operator via a "normal ordering" procedure.  This procedure is described in Peskin and Schroeder as "moving all the annihilation operators, $a_\vec{p}$, to the right of the creation operators, $a^\dagger_\vec{p}$" by hand, but this is essentially equivalent to subtracting off the vacuum contribution as, 
\begin{equation}
:\phi(x_1)\phi(x_2):\equiv \phi(x_1)\phi(x_2)-\Delta_E(x_1-x_2), \tag{2}\label{normal}
\end{equation}
so that the vacuum expectation of the normal-ordered operators vanish: $\langle 0|:\phi(x_1)\phi(x_2):0\rangle=0$, by definition.  Hence, one now has a well-defined quadratic operator,
\begin{align}
:\phi^2(x_2):\equiv \lim_{x_1\to x_2} \left\{\phi(x_1)\phi(x_2)-\Delta_E (x_1-x_2)\right\},
\end{align} 
in the sense that $\langle0|:\phi^2(x):|0\rangle$ is finite--and, in fact, zero for this choice of normal ordering.  (Of course, one could get a different vacuum expectation value via a different choice of normal ordering)
Definition of OPE coefficients: By an "operator product expansion" (OPE), I have in mind a short-distance asymptotic relation (see Scholarpedia article) which holds in any arbitrary (well-behaved) state,
\begin{align}
\langle \phi(x_1)\phi(x_2)\rangle \sim \sum_{A} C^{\phi\phi}_{A}(x_1,x_2) \langle \mathcal{O}^A(x_2) \rangle,
\end{align}
where $C^{\phi\phi}_A$ are the OPE coefficients of this expansion around $x_1\to x_2$ and the summation runs over all monomials of the fields and its derivatives.
For the free theory, it would appear such an expansion can be obtained via Taylor expanding the terms inside the normal ordering sign of \eqref{normal} as
\begin{align}
\phi(x_1)\phi(x_2)= \sum_{\alpha} \frac{(x_1-x_2)^\alpha}{\alpha!} :\partial_\alpha \phi(x_2)\phi(x_2):+\Delta_E (x_1-x_2), \label{expansion} \tag{3}
\end{align}
where $\alpha$ is a multi-index. 
Clearly, the preceding expansion holds trivially in the vacuum state.  Because the particle states of the free theory can be expressed in terms of the vacuum, I believe such an expansion will hold for an arbitrary particle state.
Reading off the nontrivial coefficients of the expansion \eqref{expansion}, I find
\begin{align}
C^{\phi\phi}_{I}(x_1,x_2)&=\Delta_E (x_1-x_2)\\
C^{\phi\phi}_{\partial_\alpha\phi\phi}(x_1,x_2)&=\frac{(x_1-x_2)^\alpha}{\alpha!}. \tag{4} \\
\end{align}
I expect to obtain more complicated OPE's in the free theory simply from applications of Wick's theorem.
Question: In the "canonical" approach, I seem to obtain the OPE coefficients for the free theory via normal ordering and Wick's theorem.  However, it's not clear to me how one defines/obtains the OPE coefficients working instead strictly in the path integral approach or even how to get started in this direction.
 A: The OPE is a statement about the asymptotic expansion when $x_1\longrightarrow x_2$ inside correlations. For this you need to introduce other fields just sitting there called spectator fields. It does not matter which formalism you use for calculating these correlations: canonical or path integral. For more mathematical precision about the OPE see my article https://arxiv.org/abs/1604.05259 If that's too much math look-up ref [114] I cited by Witten. 

Edit with a bit more detail:
The correlations for the free field are given by the Isserlis-Wick's theorem:
$$
\langle \phi(x_1)\cdots\phi(x_n)\rangle=\sum_{\mathcal{W}}\prod_{\{a,b\}\in\mathcal{W}} \langle\phi(x_a)\phi(x_b)\rangle
$$
where $\mathcal{W}$ runs over all set partitions of $\{1,2,\ldots,n\}$
which only contain blocks of size two.
One the other hand
$$
\langle :\phi(x_1)\phi(x_2):\ \phi(x_3)\cdots\phi(x_n)\rangle=
\sum_{\mathcal{W}\ {\rm nonsingular}}\prod_{\{a,b\}\in\mathcal{W}} \langle\phi(x_a)\phi(x_b)\rangle
$$
where by "nonsingular" I mean all pair partitions $\mathcal{W}$ that do not contain the pair $\{1,2\}$. For these contributions there is no divergence when you take the limit $x_1\rightarrow x_2$ with $x_2,\ldots,x_n$ distinct and fixed. The full expansion is indeed obtained by writing the Taylor expansion (in multiindex notation)
$$
\phi(x_1)=\sum_{\alpha\in\mathbb{N}_{0}^d}\frac{(x_1-x_2)^{\alpha}}{\alpha!}\ \partial^{\alpha}\phi(x_2)
$$
"inside" these expressions.
The full OPE is obtained from
$$
\langle \phi(x_1)\cdots\phi(x_n)\rangle=
\langle \phi(x_1)\phi(x_2)\rangle\ \langle\phi(x_3)\cdots\phi(x_n)\rangle+
\langle :\phi(x_1)\phi(x_2):\ \phi(x_3)\cdots\phi(x_n)\rangle
$$
namely
$$
\langle \phi(x_1)\cdots\phi(x_n)\rangle=\sum_{A}\mathcal{C}_{A}^{\phi\phi}(x_1,x_2)
\langle \mathcal{O}_{A}(x_2)\phi(x_3)\cdots\phi(x_n)\rangle
$$
with
$$
\mathcal{C}_{1}^{\phi\phi}(x_1,x_2)=\langle \phi(x_1)\phi(x_2)\rangle
$$
and
$$
\mathcal{C}_{:(\partial^{\alpha}\phi)\phi:}^{\phi\phi}(x_1,x_2)=\frac{(x_1-x_2)^{\alpha}}{\alpha!}
$$
