Suppose $x,y$ are two random variables, if $p(x|y)$ is normal distribution, and $p(y|x)$ is also normal distribution, can we say that $p(x,y)$ is also joint normal distribution?
I know that two normal random variable does not imply joint normal distribution as in： Is it possible to have a pair of Gaussian random variables for which the joint distribution is not Gaussian?
The answer is NO
The example you proposed in a comment do not work, due to the ,to me surprising, fact that $$ \int_{-\infty}^\infty \int_{-\infty}^\infty e^{-x^2 y^2} \; dx\; dy $$ do not converge, that is , it is $\infty$.
But there are many examples of bivariate densities with normal marginals, but that is not jointly normal. I will just give one such example. Define $$ f(x,y)= C e^{-[x^2+y^2 + 2xy(x+y+x y)]} $$ for $(x,y)\in \mathbb{R}^2$. Numerical integration gives $C\approx 1/4.290147838$. Since $x$ and $y$ play symmetrical roles it is enough to look at one of the conditionals. Let us first find the marginal of $Y$, which is (with the help of maple) $$ g(y) = C \frac{\sqrt{\pi}\, {\mathrm e}^{-\frac{y^{2} \left(y^{2}+2 y +1\right)}{2 y^{2}+2 y +1}}}{\sqrt{2 y^{2}+2 y +1}}$$ giving rise to the conditional $$ f(x \mid y) = \frac{{\mathrm e}^{-x^{2}-y^{2}-2 x y \left(x y +x +y \right)} \sqrt{2 y^{2}+2 y +1}}{\sqrt{\pi}\, {\mathrm e}^{-\frac{y^{2} \left(y^{2}+2 y +1\right)}{2 y^{2}+2 y +1}}} $$ which indeed is a normal density, with expectation and variance which is a function of $y$.
There is a complete monograph dedicated to this subject: Conditional Specification of Statistical Models by Barry C. Arnold, Enrique Castillo, and José María Sarabia.
