Convergence/Stability of SDE that depends on an ergodic process Let the following stochastic system be given
$$dX_t=-X_t \, dt+dW_t,$$
$$dY_t=Y_t(1-Y_t)X_t(dt+dV_t),$$
where $W_t$ and $V_t$ are independent Wiener processes, $X_0\sim\mathcal{N}(0,1/2)$ (the stationary measure), and $Y_0=1/2$.
Simulations suggest that $Y_t\rightarrow Y_\infty$ almost surely, where $Y_\infty$ is a random variable that takes the values $0$ and $1$ each with probability $1/2$.
How can this be proved?
Alternatively, convergence in probability would also work, but if possible I'm interested in the strongest possible result.
Also, can we relax the conditions on $X_t$ to be an arbitrary ergodic Markov process?

Currently, all the answers contain major gaps that need to be fixed.
Alternative ways of proving the statement are very welcome, in particular very principled or general methods of tackling this problem.
Please also note the more general version of this problem posted here on MathOverflow.

My original idea of a proof was the following:


*

*Show that the sample limit exists almost surely, i.e. show that $A=\{\omega: \lim_{t\rightarrow\infty}Y_t(\omega)\text{ exists}\}$ is measurable and that
$$\mathbb{P}(A)=1.$$

*Show that the limit is in $\{0,1\}$ almost surely if it exists
$$\mathbb{P}(Y_{\infty}=\lim_{t\rightarrow\infty}Y_t\in\{0,1\}|A)=1.$$
Because of 1), we have
$$\mathbb{P}(Y_{\infty}\in\{0,1\})=\mathbb{P}(Y_{\infty}=0)+\mathbb{P}(Y_{\infty}=1)=1.$$

*Make a symmetry argument for why $\mathbb{P}(Y_{\infty}=0)=\mathbb{P}(Y_{\infty}=1)=1/2$.

 A: Via PDE
We look at the Fokker-Plank or Komogorov forward equations. Let $p(t,x,y)$ be the probability density of the process $(X_t,Y_t)$. From the SDE, we know $p(t,x,y)=p_1(t,x)p_2(t,x,y)$ for some probability densities $p_1$ and $p_2$ and
\begin{align}
\frac{\partial}{\partial t}p_1(t,x) &= \frac{\partial}{\partial x}(xp_1(t,x))+\frac12\frac{\partial^2}{\partial x^2}p_1(t,x) \tag1 \\
\frac{\partial}{\partial t}p_2(t,x,y) &= -\frac{\partial}{\partial y}(xy(1-y)p_2(t,x,y))+\frac12\frac{\partial^2}{\partial y^2}\big((xy(1-y))^2p_2(t,x,y)\big) \tag2
\end{align}

Edit: As pointed out by @S.Surace PDE (2) dropped the terms
  involving partial derivatives of $x$. I am reconsidering this formulation.

It can be shown, not a trivial task (Maybe I will sketch the proof later.), that the dynamic solution converges to the stationary solution. It satisfies the above equations with both $\frac{\partial}{\partial t}p=0$.
Solving the first stationary equation of Equation (1), $p_1(\infty,x)\propto e^{-x^2}$. Substitute it into the stationary equation of Equation (2), and short hand $p_2(\infty,x,y)$ by $p_2$, we have
$$c(x) = -\frac1xy(1-y)p_2+\frac12\frac{\partial}{\partial y}\big((y(1-y))^2p_2\big)$$
or
$$\frac{\partial}{\partial y}p_2+2\frac{2(\frac12-y)-\frac1x}{y(1-y)}p_2-\frac{2c(x)}{(y(1-y))^2}=0.$$
The coefficient of $p$ has first order poles at $y=0$ and $y=1$, and the third term on the left hand side has a second order poles at $y=0$ and $y=1$. By the theory of ordinary differential equation, $p_2$ allows poles at $y\in\{0,1\}$. But we do need to ascertain there indeed are poles. The solution for $p_2$ is
$$p_2(x,y) = y^{-2+\frac2x}(1-y)^{-2-\frac2x}\left[2c(x)\int_{\frac12}^yds\, s^{-\frac2x}(1-s)^{\frac2x}+\frac1{16}p_2\Big(x,y=\frac12\Big)\right].$$
Consider the first term. Consider the neighbourhood of $y=0$. The analysis of the neighbourhood of $y=1$ will be similar. The highest order singularities is a first order pole when $x\not=2$. When $x=2$, the highest order singularity is $\frac{\ln y}{y}$. So the first term is never integrable at either end of $y$.
Consider the second term. Assume $p_2\big(x,y=\frac12\big)\not = 0.$ In the neighbourhood of $y=0$, it requires $x\in(0,2)$ for $p_2(x,y)$ to be integrable over $y$. In the neighbourhood of $y=1$, it requires $x\in(-2,0)$. So the second term is never integrable over the whole $y$ interval for $x\not=0$.
Therefore $p_2$ is not integrable. Thus the only stationary solution is $p_2(y)=0$. As a matter of fact one can integrate the above ODE as an algebraic function of the incomplete Beta function with exponents involving $x$.
Therefore the limiting probability distribution is $\delta\big(y(1-y)\big)$.
A: Inspired by the answer by @Hans (https://math.stackexchange.com/a/2554580/227280), we try to show that the stationary measure is not integrable.
We consider the Fokker-Planck equation of the coupled system:
$$\partial_t p(x,y,t)=\partial_x\Big(xp(x,y,t)\Big)+\frac12\partial_x^2p(x,y,t)-x\partial_y\Big(y(1-y)p(x,y,t)\Big)+\frac12x^2\partial_y^2\Big(y^2(1-y)^2p(x,y,t)\Big), \tag1$$
which when integrated over $x$ by parts (assuming that the density and its derivatives vanish for large $x$) can be written as 
$$\partial_t \int_{-\infty}^{\infty}p(x,y,t)dx=-\partial_y\Big(y(1-y)\int_{-\infty}^{\infty}xp(x,y,t)dx\Big)+\frac12\partial_y^2\Big(y^2(1-y)^2\int_{-\infty}^{\infty}x^2p(x,y,t)dx\Big). \tag2$$
Suppose that a stationary density $p(x,y)$ (sufficiently smooth) exists.
Then the stationary version of (2) reads, for $0<y<1$
$$y(1-y)\int_{-\infty}^{\infty}xp(x,y)dx-\frac12\partial_y\Big[y^2(1-y)^2\int_{-\infty}^{\infty}x^2p(x,y)dx\Big]=y(1-y)g(y)-\frac12\partial_y\Big[y^2(1-y)^2f(y)\Big]=0.$$
Solving this ODE for $f$ by integrating over $y$, we have
$$f(y)=\int_{-\infty}^{\infty}x^2p(x,y)dx=\frac{1}{y^2(1-y)^2}\Big[\frac{f(\frac12)}{16}+2\int_{\frac12}^yg(z)z(1-z)dz\Big], \quad 0<y<1.$$
Again assuming that $p(x,y)$ is well-behaved, we should also have that
$$\frac12=\int_{-\infty}^{\infty}x^2p(x)dx=\int_{\mathbb{R}\times(0,1)}x^2p(x,y)dxdy=\int_0^{1}f(y)dy.$$
It remains to be shown that the right-hand side diverges unless $f=g=0$. (As noted in the comments below, this has not been shown yet!)
We could then conclude that $p(x,y)=0$ for $0<y<1$ and that therefore the stationary measure (if it exists at all) must therefore be concentrated on the boundary of $\mathbb{R}\times[0,1]$.
A: Via SDE but is incomplete.
$$d\ln\frac Y{1-Y}=X\Big[\Big(1+\big(2Y-1\big)X\Big)dt+dV\Big].$$
We just need to show property (P) that the right-hand side when integrated has vanishing probability as $t\to\infty$ within any given arbitrarily large bound. The inverse transformation of $\ln\frac Y{1-Y}$ concentrates the measure outside the bound into intervals arbitrarily close to $\{0,1\}$. We know $V$ has Property (P) while $X$ does not vanish and approaches a stationary distribution. We need to make some estimate for the $dt$ term.
