Monge–Ampère with drift Let $I\subseteq \mathbb{R}$ be an interval. 
Let smooth $M(x,y):I\times(0,\infty) \to \mathbb{R}$ satisfies PDE:
$$
M_{xx}M_{yy}-M_{xy}^{2}+\frac{M_{y}M_{yy}}{y}=0.
$$
My question is to describe/characterize solutions of this PDE. 
Here are few nontrivial examples which solve the PDE:
\begin{align}
&1)\quad M(x,y)=x\ln x-\frac{1}{2}\frac{y^{2}}{x};\\
&2)\quad M(x,y)=x^{2}-y^{2};\\
&3)\quad M(x,y)=\sqrt{[\Phi'(\Phi^{-1}(x))]^{2}+y^{2}} \quad \text{where}\quad \Phi(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-\frac{t^{2}}{2}}dt.
\end{align}
Update 1:
I have few more questions regarding Robert Bryant's answer:
1) (Philosophical question) Which PDE's $F(M, M_x, M_y, M_{xx}, M_{xy}, \ldots)=0$ can be "encoded" as an exterior differential system on $\mathbb{R}^{m}$? 
2)  (Incorrect question) After having  exterior differential system, how do you come up with "right" change of variables? Does it mean that you tried all possible (and reasonable) ones  and this was the best one? Maybe this question is not correct because I did not specify what is the meaning of ``right change of variables''. In this case goal was backward heat equation. 
3) (Technical question) What happens if $p<0$? I think nothing changes all you need is requirement that $p\neq 0$. Am I Right?
4) (The most important question) Does this "local" knowledge actually allow me to write down at least one global nontrivial solution? For example, can you somehow trace back to the solution 3) that I mentioned above?  You can assume some boundary conditions. 
Update 2:
1) Accepted. 
2) Accepted. 
3) Accepted. P.S. I had in my mind to consider $-M(x,y)$ instead of $M(x,y)$. 
4) ACCEPTED! By the way my first solution corresponds to the case $u(p,t)=\frac{t}{2}-\frac{p^{2}}{4}$. My second solution corresponds to the case $u(p,t)=-e^{-t-p-1}$. 
 A: Perhaps the following observations will be of use to you:  First, for any (local) solution $M(x,y)$ of your equation on a simply-connected open domain $D\subset \mathbb{R}\times(0,\infty)$ in the $xy$-plane, consider its $1$-graph 
$$
(x,y,p,q)=\bigl(x,y,M_x(x,y),M_y(x,y)\bigr)
$$
in $xypq$-space.  This is a simply-connected surface $\Sigma$ in $4$-space on which $\Omega = \mathrm{d}x\wedge\mathrm{d}y$ is nonvanishing, but to which the two $2$-forms
$$
\Upsilon_1 = \mathrm{d}p\wedge\mathrm{d}x + \mathrm{d}q\wedge\mathrm{d}y
\qquad\text{and}\qquad
\Upsilon_2 = \bigl(y\,\mathrm{d}p+ q\,\mathrm{d}x\bigr)\wedge\mathrm{d}q
$$
pull back to be zero.  
Conversely, suppose given a simply connected surface $\Sigma$ in $xypq$-space (with $y>0$) on which $\Omega$ is nonvanishing but to which $\Upsilon_1$ and $\Upsilon_2$ pullback to be zero.  The $1$-form $p\,\mathrm{d}x + q\,\mathrm{d}y$ pulls back to $\Sigma$ to be closed (since $\Upsilon_1$ vanishes on $\Sigma$) and hence exact, and therefore there exists a function $m:\Sigma\to\mathbb{R}$ such that $\mathrm{d} m = p\,\mathrm{d}x + q\,\mathrm{d}y$ on $\Sigma$.  We then have (at least locally), $m = M(x,y)$ on $\Sigma$ and, by its definition, we have $p = M_x(x,y)$ and $q = M_y(x,y)$ on the surface.  Then the fact that $\Upsilon_2$ vanishes when pulled back to $\Sigma$ implies that $M(x,y)$ satisfies the desired equation.
Thus, we have encoded the given PDE as an exterior differential system on $\mathbb{R}^4$.  Note, though, that we can make a change of variables on the open set where $q>0$:  Set $y=qr$ and let $t=\tfrac12 q^2$.  Then, using these new coordinates on this domain, we have
$$
\Upsilon_1 = \mathrm{d}p\wedge\mathrm{d}x + \mathrm{d}t\wedge\mathrm{d}r
\qquad\text{and}\qquad
\Upsilon_2 = \bigl(r\,\mathrm{d}p+ \mathrm{d}x\bigr)\wedge\mathrm{d}t.
$$
Now, when we take an integral surface $\Sigma$ of these $2$-forms on which $\mathrm{d}p\wedge\mathrm{d}t$ is nonvanishing, it can be written locally as a graph of the form 
$$
(p,t,x,r) = \bigl(p,t,u_p(p,t),u_t(p,t)\bigr)
$$
(since $\Sigma$ is an integral of $\Upsilon_1$), where $u(p,t)$ satisfies $u_t + u_{pp} = 0$ (since $\Sigma$ is an integral of $\Upsilon_2$).  Thus, 'generically', your equation is equivalent to the backwards heat equation, up to a change of variables.   
(Note, by the way, that an integral surface of the $\Upsilon_i$ on which $\mathrm{d}p\wedge\mathrm{d}t$ vanishes identically must also have $\mathrm{d}t$ vanishing identically, since, otherwise, because $\mathrm{d}x\wedge\mathrm{d}t$ vanishes identically (since $\Upsilon_2$ vanishes), one would have $\mathrm{d}x$ and $\mathrm{d}p$ be multiples of $\mathrm{d}t$, which would then force $\mathrm{d}r\wedge\mathrm{d}t$ to vanish identically (since $\Upsilon_1$ vanishes), and hence each of $x,p,r$ would be functions of $t$ and this is impossible for a surface.  Thus, integral surfaces on which $\mathrm{d}p\wedge\mathrm{d}t$ vanishes identically are necessarily locally of the form
$$
(p,t,x,r) = \bigl(P(s),t_0,X(s),r\bigr)
$$
where $t_0$ is a constant and $r$ and $s$ are coordinates on the surface and $\bigl(P'(s),X'(s)\bigr)$ is never $(0,0)$.)
Thus, in some sense, your nonlinear Monge-Ampère equation is simply the (linear) backwards heat equation in disguise.  You should be able to use this to generate many nontrivial solutions and to solve the appropriate characteristic initial value problem.
Responses to updated question
1) In principle, any PDE can be encoded as an exterior differential system, often in more than one way.  There is a standard way to do it for Monge-Ampère equations, though, so that such an equation for one function of two variables can be written as an EDS on $\mathbb{R}^5$. For example, see the book Exterior Differential Systems by Bryant, et al for how this is usually done.
2) The first step is to determine the equivalence class of the EDS up to diffeomorphism (i.e., change of variables).  Not all Monge-Ampère equations are equivalent, but there are tests for when it is parabolic, when it is equivalent to a linear equation, etc.  In this particular case, I computed the conversation laws for the equation and realized that this space had infinite dimension, so I knew that it must be equivalent to a linear equation, in fact, the backwards heat equation.  For example, see my paper with Phillip Griffiths, Characteristic cohomology of differential systems II: Conservation laws for a class of parabolic equations, Duke Math. J. Volume 78, Number 3 (1995), 531-676.
3)  I think you mean $q<0$, not $p<0$.  In that case, yes, in fact, the mapping $(x,y,p,q)\mapsto(x,-y,p,-q)$ switches the domains $q>0$ and $q<0$ in $\mathbb{R}^4$ and preserves the EDS, so they are equivalent.
4) Yes.  Here are a couple of examples.  First,, if you follow the transformation through, you'll see that a solution $u(p,t)$ of $u_t+u_{pp}=0$ on some domain in $pt$-space corresponds to a solution of your original equation whose graph in $xyz$-space is the parametrized surface
$$
\begin{aligned}
x &= u_p(p,\tfrac12q^2)\\
y &=q u_t(p,\tfrac12q^2)\\
z &= p\,u_p(p,\tfrac12q^2) +q^2\,u_t(p,\tfrac12q^2) - u(p,\tfrac12q^2)
\end{aligned}
$$
You just need to choose your solution $u(p,t)$ on a $pt$-domain so that this surface is the graph of the form $z = M(x,y)$ over an appropriate domain and then $M$ will satisfy your
equation.  
For a simple example, take $u(p,t) = pt - \tfrac16 p^3$.  Running this through the above process yields in the solution
$$
M(x,y) = \tfrac23\bigl(\sqrt{x^2+y^2} + 2x\bigr)\bigl(\sqrt{x^2+y^2} - x\bigr)^{1/2}.
$$
As another example, take $u(p,t) = e^t\,\sin p$ on the open $pt$-strip $ \tfrac12\pi < p < \pi$ and $t>0$. Then the above surface in $xyz$-space is a graph $z = M(x,y)$ over the second quadrant in the $xy$-plane, so this gives a nontrivial solution $M(x,y)$ in the $xy$-domain where $x<0$ and $y>0$.  I believe that this solution is not the same as any of the solutions you have written down above.  In fact, I don't see how to get an 'explicit' formula for $M(x,y)$, i.e., to eliminate the parameter variables $p$ and $q$.
