This is not really an answer, just a sequence of comments that are all related, but are too long to put into a comment field.
First, some good news:
When $n=1$, there's always a (unique) solution for the desired $g$: Just note that the equation in this case reduces to
$$
g''(x) = \frac{\phi(x)}{\phi'(x)^2},
$$
and the hypotheses imply that the right hand side extends smoothly to $x=0$ because, by Whitney (or the Morse Lemma), we can write $\phi(x) = x^2\psi(x)$ where $\psi$ is smooth and positive on a neighborhood of $x=0$, and the smoothness of the right hand side follows from this. Of course, there will be a unique solution $g$ to the above ODE with $g(0)=g'(0)=0$.
For any $n$ and $\phi$ defined on a neighborhood of $0\in\mathbb{R}^n$ satisfying the given hypotheses, there exists a (unique) Taylor series that formally solves the given equation for $g$. Thus, if there exist smooth solutions in a neighborhood of $0\in\mathbb{R}^n$, then the difference of any two such vanishes to infinite order at the origin. Moreover, if $\phi$ is real-analytic, then the (unique) formal Taylor series solution for $g$ converges in a neighborhood of $0\in\mathbb{R}^n$ to the (unique) real-analytic solution. The proof of this is not very hard: You can always diagonalize the matrix of the quadratic terms in $\phi$ by an orthogonal rotation (which does not disturb the equation), and then a simple argument by induction shows how the formal Taylor series can be solved for uniquely order by order. Once one sees how this goes, one can then, with some work (that is a little messy, but not hard), show that the method of majorants applies to prove that if $\phi$ is real-analytic then the formal Taylor series for $g$ has a positive radius of convergence and hence solves the equation. Thus, when $\phi$ is real-analytic, there is a unique real-analytic solution.
Second, some not-so-good news:
Contrary to what was suggested in the comments, for general $\phi$, the equation does not reduce to a second order ODE along the integral curves of $\nabla\phi$. In fact, if the integral curves of $\nabla\phi$ are not straight lines (which is almost always the case), this is a degenerate PDE that has characteristics of the heat equation.
Consider the first nontrvial case, $n=2$. In a neighborhood of a point $p\not=0\in\mathbb{R}^2$ where the curvature at $p$ of the integral curve of $\nabla\phi$ passing through $p$ is nonzero, one can show that one can introduce new coordinates $(x,t)$ centered on $p$ in which the equation becomes
$$
g_{xx}(x,t) + a(x,t)\ g_{x}(x,t) + b(x,t)\ g_t(x,t) = \phi(x,t)
$$
where $a$ and $b$ are smooth functions and $b(0,0)<0$. (Essentially, in these coordinates, $\nabla\phi = \partial_x$.) Solving for $g_t$, one sees that this is a linear, inhomogeneous heat equation for $g$ in this neighborhood. However, I don't think that there will be much in the literature on the heat equation that will help you deal with the singularity or the cases in which the integral curves of $\nabla\phi$ have inflection points.
In higher dimensions, this kind of degenerate heat equation does not appear to be much studied at all, so I don't think you'll get much mileage out of the PDE literature there either.
About the only thing that this discussion might help you with if you want to prove the existence of a smooth solution is that you at least know (via Whitney) that there is a smooth function $g$ that solves the equation to infinite order at the origin. Whether you can somehow use this to prove that a smooth solution exists, I won't speculate.
There is one tiny ray of sunshine in the gloom: In the case that $\phi(x)=\tfrac{c}2\,|x|^2$ where $c>0$, the gradient lines of $\nabla\phi$ are straight lines, and the ODE argument does apply, showing that there is a unique smooth solution, namely $g = \tfrac1{4c}\,|x|^2 = \frac\phi{2c^2}$. (By the way, Kofi, you need to check your constants in your linear solution. Remember that, when $\phi = |x|^2$, you have $\nabla\phi(x) = 2x$, so $\nabla^2g = \tfrac14\,I$, so $g = \tfrac18\,|x|^2$, not $\tfrac12\,|x|^2$.)
Finally, here is a suggestion about a possible way to proceed: Suppose one turns the equation around and considers it as an equation for $\phi$ when one is given a smooth strictly convex function $g$ defined on a neighborhood of $0\in\mathbb{R}^n$ and having a nondegenerate minimum there. Let $h = h^{ij}$ be the (positive definite) inverse matrix of the Hessian of $g$, and consider the metric $ds^2 = h^{ij}\,\mathrm{d}x^i\mathrm{d}x^j$. Then your equation becomes the first-order PDE
$$
|\mathrm{d}\phi|^2 = \phi
$$
where the norm is computed with respect to the metric $ds^2$. Thus, this is a version of the eikonal equation for the metric $ds^2$, and it is well-known (by ODE techniques, essentially the method of Cauchy characteristics) that $\phi(x) = \tfrac14 d(x,0)^2$ is the only smooth solution of this equation in a neighborhood of $0\in\mathbb{R}^n$ that has a nondegenerate minimum there (and, conversely, that, for any strictly convex $g$, this is a smooth solution of the above equation). Here, $d(x,0)$ is the distance from $x$ to $0$ measured with respect to the Riemannian metric $ds^2$.
Now, the map that sends a smooth-strictly-convex-on-a-neighborhood-of-the-origin function $g$ to the squared distance function $\phi = \tfrac14\,d(x,0)^2$ for the metric $ds^2$ associated to $g$ as above will be differentiable when appropriate smooth (Frechet) structures on the domain and range are specified, and it might even be invertible (as it certainly is, by the argument above, if one restricts to the analytic category), so one might be able to solve the original problem by applying an appropriate implicit function theorem. This is really an analysis question, though, so good luck.
