Going from conditional distributions to the joint distribution function Assume that we have a random vector $\mathbf{Y}=(Y_1,...,Y_n)$ with unknown joint distribution. We know however that the conditional distribution of an entry $Y_i$ given the rest of the vector, without the entry $Y_i$ ($Y_{-i}$) is normal:
$$ Y_i|Y_{-i}\sim N(\mu+\beta Y_kY_j,\sigma^2) $$
with $k,j\neq i$ standing in some relation with $i$. For example $k=i-1, j=k+1$ and $\mu$ some constant. Is it possible to find the joint distribution function? I.e. if for example we would have
$$ Y_i|Y_{-i}\sim N(\mu+\beta Y_k,\sigma^2) $$
then the joint probability function would be a multivariate normal (with some resitriction on the parameter space (i.e. $\beta<1$), but what about the case above? I do understand that the joint distribution cannot be normal but isn't it possible that the joint distribution is just another multivariate distribution?
A short Gibbs sampling code in R shows that it is at least plausible that such a distribution exists as the sampler does not diverge (depending on the parameter space):
Y<-1
Z<-1
col<-c()
for (i in 1:10000){
X<-rnorm(1,1+0.1*Y*Z,1)
Y<-rnorm(1,1+0.1*X*Z,1)
Z<-rnorm(1,1+0.1*X*Y,1)
col<-cbind(col,c(X,Y,Z))
}

scatterplot3d::scatterplot3d(x=col[1,],y=col[2,],z=col[3,])

 A: To determine whether or not conditional distributions are compatible with one another, an approach is to apply the Hammersley-Clifford Theorem:

Under the condition that the supports of the full conditionals $g_{j}$
  are the entire sets $\mathcal{Y}_j$, the joint distribution $g$
  satisfies $$ g(y_1,\ldots,y_p) \propto \prod_{j=1}^p \;
  {g_{\ell_j}(y_{\ell_j}|y_{\ell_1},
  \ldots,y_{\ell_{j-1}},y_{\ell_{j+1}}^\prime,\ldots,y_{\ell_p}^\prime)
  \over g_{\ell_j}(y_{\ell_j}^\prime|y_{\ell_1},\ldots,y_{\ell_{j-1}},
  y_{\ell_{j+1}}^\prime,\ldots,y_{\ell_p}^\prime)} $$ for every
  permutation $\ell$ on $\{1,2,\ldots,p\}$ and every $y,y' \in
  \mathcal{Y}=\prod_{j=1}^p \mathcal{Y}_j$.

In the current setting, it however is impossible to find a joint distribution: the products $y_i y_j y_k$ appearing in the joint mean that the joint cannot be Normal. But since all conditionals are Normal, there is no other joint distribution possible.
Take as an example
\begin{align*}X|Y,Z&\sim\mathcal{N}(\beta YZ,1)\\
Y|X,Z&\sim\mathcal{N}(\beta XZ,1)\\
Z|Y,X&\sim\mathcal{N}(\beta XY,1)
\end{align*}
and apply Hammersley-Clifford's formula with $x'=y'=z'=0$. The joint does not depend on the choice of $(x',y',z')$ and this choice makes a lot of terms vanish. The log joint density should be (up to an additive constant)
$$-\{x^2+y^2+(z-\beta xy)^2-\beta^2 x^2 y^2\}/2$$
since all other terms cancel. Which means that the joint density is proportional to
$$\exp(-\{x^2+y^2+(z-\beta xy)^2-\beta^2 x^2 y^2\}/2)$$
If I integrate this density in $z$, the term $(z-\beta xy)^2$ vanishes and integrate to a constant in $(x,y)$. I am thus left with
$$\exp(-\{x^2+y^2-\beta^2 x^2 y^2\}/2)$$
which does not integrate in $x$ conditional on $y$ when $y^2>\beta^{-2}$. Therefore there cannot exist a joint probability density.
