Agreeing with the answers above, I am also under the impression that path integrals are a mathematical tool whose reach extends beyond the real of quantum mechanics. However, it seems that, at least naively, path integrals lead always to an operator formulation of the type found in quantum mechanics. It won't be unitary, but is still very useful for the study of the system, as exemplified by the use of the transfer matrix method in classical spin chains. Moreover, it can be potentially made unitary, possibly after a Wick rotation.
The statement above is well known, see for example https://www2.perimeterinstitute.ca/videos/lecture-3-boundary-conditions-and-extended-defects. However, let me try to expand on the material covered in the video (I do believe the stuff below should be scattered throughout the literature but I have failed to find a reference that contains it, if anyone finds it, let me know!). The basic idea is the following. To every codimension 1 surface $\Sigma$ we have boundary conditions we can put on our path integrals. We can then consider the Hilbert space $\mathcal{H}_\Sigma$ "spanned" by these boundary conditions. By these we mean that the Hilbert space is some space of functionals $\Psi(\varphi)=\langle\varphi|\Psi\rangle$ of the boundary conditions $\varphi$. We can say that this space is spanned by "vectors of the form $|\varphi\rangle$ and that they form an orthonormal basis, in the sense that the space of boundary conditions has a measure $D\varphi$ such that
$$\langle\varphi'|\varphi\rangle=\delta(\varphi'-\varphi),$$
where this delta function with respect to the measure and
$$1=\int D\varphi\,|\varphi\rangle\langle\varphi|.$$
In particular, operators on this space $\mathcal{O}$ are completely determined by its matrix entries
$$\langle\varphi'|\mathcal{O}|\varphi\rangle$$
since by inserting a resolution of the identity, we have
$$\langle\varphi|\mathcal{O}|\Psi\rangle=\int D\varphi'\,\langle\varphi|\mathcal{O}|\varphi'\rangle\langle\varphi'|\Psi\rangle.$$
Now we can consider a path integral on a "spacetime" $M$ with an incoming boundary $\Sigma_i$ and an outgoing boundary $\Sigma_f$. Then the path integral defines a "time evolution" operator $U(M):\mathcal{H}_{\Sigma_i}\rightarrow\mathcal{H}_{\Sigma_f}$ given by
$$\langle\varphi_f|U(M)|\varphi_i\rangle=\int_{\{\phi|\phi|_{\Sigma_i}=\varphi_i,\phi|_{\Sigma_f}=\varphi_f\}}\mathcal{D}\phi\,e^{-\frac{1}{\hbar}S(\phi)}=:\langle 1\rangle,$$
We have introduced this correlation function notation to avoid continuously writing path integrals but one shouldn't forget that these correlations depend on the spacetime $M$ and the boundary conditions. We can use this operator on a thin cylinder $M=\Sigma\times[0,\epsilon]$ to perform "quantization". Given some observable $\mathcal{o}$ on this surface, we can construct its operator counterpart $\mathcal{O}$ as follows. In order to define the state $\mathcal{O}|\phi\rangle$ we will demand that for all final boundary conditions $\phi_f$ and small $\epsilon$ we have
$$\langle\varphi'|U(M)\mathcal{O}|\varphi\rangle=\int_{\{\phi|\phi|_{\Sigma\times\{0\}}=\varphi,\phi|_{\Sigma\times\{\epsilon\}}=\varphi'\}}\mathcal{D}\phi\,e^{-\frac{1}{\hbar}S(\phi)}\mathcal{o}(\phi)=\langle\mathcal{o}\rangle.$$
We can in fact obtain an explicit expression for the matrix elements of the operator if we assume that in the limit $\epsilon\rightarrow 0$ we have $U(M)\rightarrow 1$, yielding the "quantization" formula
$$\langle\varphi'|\mathcal{O}|\varphi\rangle=\lim_{\epsilon\rightarrow 0}\int_{\{\phi|\phi|_{\Sigma\times\{0\}}=\varphi,\phi|_{\Sigma\times\{\epsilon\}}=\varphi'\}}\mathcal{D}\phi\,e^{-\frac{1}{\hbar}S(\phi)}\mathcal{o}(\phi)=\lim_{\epsilon\rightarrow 0}\langle\mathcal{o}\rangle.$$
One should convince oneself that this does work in a simple example. Following the video above, we can use the 1d scalar field with action
$$S(\phi)=\int\text{d}t\frac{1}{2}\dot{\phi}^2$$. In 1d the only connected codimension 1-surface is a point. Thus, the Hilbert space is the space $\mathcal{H}=L^2(\mathbb{R})$ wave functions $\Psi(\varphi)$ of the values $\varphi$ our scalar field $\phi$ can take at said point. We will be interested in performing path integrals over $M=[0,\epsilon]$.
There are two fundamental tools we will use to compute the associated path integrals. One is that the standard path integral computing the time evolution operator can be explicitly done using, for example, zeta-function regularization
$$\langle1\rangle=\int_{\{\phi:[0,\epsilon]\rightarrow\mathbb{R}|\phi(0)=\varphi_i,\phi(\epsilon)=\varphi_f\}}\mathcal{D}\phi\,e^{-\frac{1}{\hbar}S(\phi)}=\frac{1}{\sqrt{2\pi\hbar\epsilon}}e^{-\frac{1}{2\hbar}\frac{(\phi_f-\phi_i)^2}{\epsilon}}=:K(\phi_f,\epsilon;\phi_i,0).$$
In the limit as $\epsilon\rightarrow 0$ this indeed tends to $\delta(\phi_f-\phi_i)$, as we assumed in our general explanation above. The second tool is the Schwinger-Dyson equation,
$$\left\langle\frac{\delta S}{\delta\phi(t)}F(\phi)\right\rangle=\hbar\left\langle\frac{\delta F}{\delta\phi(t)}\right\rangle,$$
for any observable $F$. In our specific case the equations of motion are
$$\frac{\delta S}{\delta\phi(t)}=-\ddot{\phi}(t),$$
so that the Schwinger-Dyson equation reduces to
$$\langle\ddot{\phi}(t)F(\phi)\rangle=-\hbar\left\langle\frac{\delta F}{\delta\phi(t)}\right\rangle,$$
In particular, taking $F=1$, we obtain that
$$0=\langle\ddot{\phi}(t)\rangle=\frac{d^2}{dt^2}\langle \phi(t)\rangle,$$
i.e. the equations of motion hold as long as there are no other operator insertions. Moreover, it is clear from the boundary conditions in the path integral that
$\langle \phi(0)\rangle=\varphi\langle 1\rangle$ and $\langle\phi(\epsilon)=\varphi'\langle 1\rangle$. We conclude that
$$\langle\phi(t)\rangle=\left(\varphi+\frac{\varphi'-\varphi}{\epsilon}t\right)\langle 1\rangle$$
With these tools in hand we can compute the operator counterparts of different observables. Let us start with the operator corresponding to $\phi(0)$, which yields
$$\langle\varphi'|\Phi|\varphi\rangle=\lim_{\epsilon\rightarrow 0}\langle\phi(0)\rangle=\varphi\lim_{\epsilon\rightarrow 0}\langle 1\rangle=\varphi\delta(\varphi'-\varphi).$$
In other words
$$\langle\varphi|\Phi|\Psi\rangle=\varphi\langle\varphi|\Psi\rangle,$$
i.e. it is simply the multiplication operator.
Let us now move the the operator corresponding to $\dot{\phi}(0)$
$$\langle\varphi'|\Pi|\varphi\rangle=\lim_{\epsilon\rightarrow 0}\langle\dot{\phi}(0)\rangle=\lim_{\epsilon\rightarrow 0}\frac{\varphi'-\varphi}{\epsilon}\langle 1\rangle.$$ To compute this limit we note that
$$\frac{d}{d\varphi'}\langle 1\rangle=-\frac{1}{\hbar}\frac{\varphi'-\varphi}{\epsilon}\frac{1}{\sqrt{2\pi\hbar\epsilon}}e^{-\frac{1}{2\hbar}\frac{(\varphi'-\varphi)^2}{\epsilon}},$$
so that
$$\langle\varphi'|\Pi|\varphi\rangle=\lim_{\epsilon\rightarrow 0}\langle\dot{\phi}(0)\rangle=\lim_{\epsilon\rightarrow 0}-\hbar\frac{d}{d\varphi'}\langle 1\rangle=-\hbar\frac{d}{d\varphi'}\delta(\varphi'-\varphi).$$
We conclude that
$$\langle\varphi|\Pi|\Psi\rangle=-\hbar\frac{d}{d\varphi}\langle\varphi|\Psi\rangle,$$
i.e. $\Pi$ is the derivative operator.
We are of course recovering the theory of a free particle, as expected. I just decided to call it $\phi$ instead of $x$ or $q$ to highlight the fact that a Gaussian classical statistical field theory in 1d recovers this (Wick rotated) operator formulation. We can go ahead with these methods and compute other interesting things. For example, we can define the Hamiltonian
$$H=-\hbar\frac{d}{d\epsilon}U(M)|_{\epsilon=0}.$$ By taking now two derivatives with respect to $\varphi'$, one can then show that
$$\langle\varphi|H|\Psi\rangle=-\frac{1}{2}\frac{d^2}{d\varphi^2}\langle\varphi|\Psi\rangle,$$
as expected. We can also for example see that the operator corresponding to $\ddot{\phi}(0)$ vanishes due to the equations of motion.
Some final remarks are in order:
- It is interesting to repeat this in a massive field theory with the additional term $\frac{1}{2}m\phi^2$. One should recover the usual operator formulation of the harmonic oscillator. Moreover, in this case, the operator corresponding to $\ddot{\phi}$ is $m\Phi$ by the equations of motion.
- This procedure clarifies the relationship between canonical quantization and path integrals. This is particularly clear if one thinks of phase space not as the set of initial conditions (positions and momenta), but rather takes the covariant perspective that it is the set of solutions to the equations of motion. We have relations like $\ddot{\phi}=0$ in the free particle case relating the restrictions of observables to this phase space. This gets implemented at the quantum level by the Schwinger-Dyson equation and our quantization formula, explaining why the starting point for quantization is the set of observables in phase space. The symplectic structure then gets related to the commutators by the mechanism explained in the video above.
- One might think that this quantization scheme achieves the goal of assigning to every classical observable a quantum operator. Sadly, I don't think this is the case. At least in $d>1$ expressions like $\langle\phi(0)^2\rangle$ are quite simply ill-defined and need to be regularized, which shows that quantizing $\phi(0)^2$ is not straightforward and requires normal ordering.
