Lack of Maslov index in the path integral formalism Introduction
Consider Feynman's famous path integral formula
\begin{equation}
K(x_a,x_b) = \int \mathcal{D}[x(t)] \exp \left[ \frac{i}{\hbar} \int_{t_a}^{t_b} dt \, \mathcal{L}(x(t),\dot{x}(t),t) \right] \, ,
\end{equation}
where the $\mathcal{D}[x(t)]$ represents the mathematically problematic "sum over paths" and $x_a \equiv x(t_a)$, $x_b \equiv x(t_b)$. The semiclassical limit of the above expression, which can be defined as the limit $\hbar \to 0$, can be approximated by the saddle point method and results in
\begin{equation}
K_{SC}(x_a,x_b) = \left( \frac{1}{2 \pi i \hbar} \right)^{-\frac{n}{2}}  \sum_{branches} \bigg| \det \frac{\partial x_b}{\partial p_a} \bigg|^{-\frac{n}{2}}\exp \left[\frac{i}{\hbar} S(x_a,x_b,t) \right] \, ,
\end{equation} 
where $SC$ stands for semiclassical, $S$ is the classical action, $n$ is half the phase space dimension, $p_a$ is the initial momentum and the sum is over all "branches" patched by points where $\dot{x}(t)=0$ (these points are examples of what we call caustics). Now, numerous studies conducted by physicists in the past have established that this formula cannot be right for all times, even for quadratic Hamiltonians, where the semiclassical approximation should be exact. The rigorous reason for that was first exposed in (1) and has to do with cohomology theory: the patching of Lagrangian submanifolds between caustics is accompanied by a jump in the determinant squared mapping, which is intrinsically related to the Morse (or Maslov) index acting on the Lagrangian Grassmannian.
Now, since the semiclassical propagator can be derived from more rigorous formalisms than the path integral, the proper form of $K_{SC}$ has been obtained by Gutzwiller (2) and rigorously justified by many other mathematicians (see (1) again). It is
\begin{equation}
K_{SC}(x_a,x_b) = \left( \frac{1}{2 \pi i \hbar} \right)^{-\frac{n}{2}} \sum_{branches} \bigg| \det \frac{\partial x_b}{\partial p_a} \bigg|^{-\frac{n}{2}}\exp \left[\frac{i}{\hbar} S(x_a,x_b,t) - \frac{i \pi \nu}{2} \right] \, ,
\end{equation} 
where $\nu$ is the Maslov index, related to the cohomology class of the trajectories being considered. This is a celebrated result called Gutzwiller's formula.
Question
If we consider points $x_a$ and $x_b$ where $p$ is a single valued function of $x$, then the determinant in $K_{SC}$ has no jumps, we are away from caustics and the Maslov index is zero. Now, if $p$ is multivalued we have to account for the branches, which I think cannot be directly derived from Feynman's path integral. This means that for each and every system where the classical trajectories have a non-trivial cohomology group, the path integral is wrong. Now, the most ridiculous example of that would be the harmonic oscillator, whose exact quantum propagator (identical to its semiclassical approximation, because the Hamiltonian is quadratic) can be seen here. Now, if I am correct in my reasoning, then this propagator is wrong. I checked Feynman & Hibbs and there's a footnote that says that this propagator is indeed wrong for long times, that is, after the momentum has a caustic. They refer to an article that pretty much guesses the right answer (3). 
My question is: can the proper Gutzwiller formula be derived from Feynman's path integral? If this is not the case, then each and every propagator calculated from the path integral formula will be wrong for Hamiltonians with non-linear potentials.
P.S.: I have the feeling that this is related to the fact that even if physicists say that "the path integral sums over all paths", it does not sum over paths with different cohomology classes. 

References:
(1) Characteristic Class Entering in Quantization Conditions, by Vladimir Arnol'd
(2) Chaos in Classical and Quantum Mechanics, by Martin Gutzwiller
(3) Propagator for the simple harmonic oscillator, by Thorber and Taylor, American Journal of Physics
 A: *

*Apart from a full proof of the Gutzwiller's formula in the context of the Feynman path integral (FPI), then OP is essentially asking
the following:

Does the FPI know about the metaplectic correction/Maslov index and caustics? 

Answer: Yes, the FPI does contain these semiclassical phase factors. 

*In practice, let $$\Delta t_M~:=~t_f-t_i~>~0.\tag{1}$$ 
To render the oscillatory FPI convergent, insert the Feynman's $i\epsilon$ prescription 
$$\Delta t_M\quad\to\quad\Delta t_M-i\epsilon.\tag{2}$$ 
Or equivalently, Wick-rotate $$\Delta t_E~:=~i\Delta t_M,\tag{3}$$ where ${\rm Re}(\Delta t_E)>0$. Here the letters $M$ and $E$ stands for Minkowski and Euclid, respectively.

*Let us for simplicity assume there are no instantons. The WKB/stationary phase approximation yields for $\hbar \to 0$:
$$\langle q_f, t_f | q_i,t_i \rangle 
~=~ \int_{q(t_i)=q_i}^{q(t_f)=q_f} \! {\cal D}q~\exp\left(\frac{i}{\hbar} S\right) 
~\stackrel{\text{WKB}}{\sim}~ 
\frac{\exp\left(\frac{i}{\hbar} S_{\rm cl}\right)}{\sqrt{{\rm Det}H}}.\tag{4}$$

*The Gelfand-Yaglom formula yields that the functional determinant is given by
$$\frac{{\rm Det}H}{{\rm Det}H^{(0)}}~=~\frac{\prod_{n\in\mathbb{N}}\left(1- \left(\frac{\Delta t_M}{T_n}\right)^2\right)}{\prod_{n\in\mathbb{N}}\left(1- \left(\frac{\Delta t_M}{T_n^{(0)}}\right)^2\right)}, \tag{5}$$
where $T_n$ denotes the times for caustics. They appear squared because of time-reflection symmetry. The superscript $(0)$ denotes the corresponding free theory with $T_n^{(0)}=\infty$ and 
$$ {\rm Det}H^{(0)}~=~ \frac{2\pi i\hbar\Delta t_M}{m} .\tag{6}$$
Hence
$$ {\rm Det}H~\stackrel{(2)+(3)}{=}~ \frac{2\pi i\hbar\Delta t_M}{m}
\prod_{n\in\mathbb{N}}\left(1- \frac{\Delta t_M}{T_n}\right)\left(1+ \frac{\Delta t_M}{T_n}\right) .\tag{7}$$ 

*The Maslov index comes as follows (keeping in mind that this is an analytic continuation from Euclidean time):


*

*We get a factor 
$$\exp\left(-\frac{i\pi}{4}\right)\tag{8}$$ 
in the FPI (4) from the ${\rm Det}H^{(0)}$ determinant (6).

*Note that for $0<\Delta t_M < T_1$, the infinite product (7) is positive, and we haven't passed any caustics yet. However for each new caustic $T_n$ (that $\Delta t_M$ passes), the functional determinant (7) develops a new negative eigenvalue, and we get another factor 
$$\exp\left(-\frac{i\pi}{2}\right)\tag{9}$$
in the FPI (4). In other words, the Morse index of the pertinent Hessian becomes relevant. 


*Example. Consider the quantum harmonic oscillator (QHO) and with characteristic frequency
$$ \frac{2\pi}{T}~=~\omega~=~\sqrt{\frac{k}{m}},\tag{10}$$
and caustics
$$ T_n ~=~n\frac{T}{2} \tag{11}$$
for every half-period.
Altogether, the caustic-crossing induces the well-known metaplectic correction phase factor
$$\exp\left(-\frac{i\pi}{2}\left(\frac{1}{2}+\left[\frac{2\Delta t_M}{T}\right]\right)\right),\tag{12}$$
where $\left[x\right]$ denotes the integer part of $x$.   
References:


*

*R.P. Feynman & A.R. Hibbs, Quantum Mechanics and Path Integrals, 1965; Section 3-11.

*J. Polchinski, String Theory Vol. 1, 1998; Appendix A.
