Transformation of random variables and joint distributions

I don't know how to get Mathematica to get the joint distribution explicitly for a general value of $n$ but here is how one can easily see the pattern to figure out the general solution.

First the distribution of the mean:

marginalDistribution = TransformedDistribution[Sum[y[i], {i, n}]/n, 
  Table[y[i] \[Distributed] NormalDistribution[0, \[Sigma]], {i, n}], 
  Assumptions -> \[Sigma] > 0]
{#, marginalDistribution/.n->#} &/@Range[2,10]

$$ \begin{array}{cc} 2 & \text{NormalDistribution}\left[0,\frac{\sigma }{\sqrt{2}}\right] \\ 3 & \text{NormalDistribution}\left[0,\frac{\sigma }{\sqrt{3}}\right] \\ 4 & \text{NormalDistribution}\left[0,\frac{\sigma }{2}\right] \\ 5 & \text{NormalDistribution}\left[0,\frac{\sigma }{\sqrt{5}}\right] \\ 6 & \text{NormalDistribution}\left[0,\frac{\sigma }{\sqrt{6}}\right] \\ 7 & \text{NormalDistribution}\left[0,\frac{\sigma }{\sqrt{7}}\right] \\ 8 & \text{NormalDistribution}\left[0,\frac{\sigma }{2 \sqrt{2}}\right] \\ 9 & \text{NormalDistribution}\left[0,\frac{\sigma }{3}\right] \\ 10 & \text{NormalDistribution}\left[0,\frac{\sigma }{\sqrt{10}}\right] \\ \end{array} $$

So we see that the marginal distribution of $\bar{y}$ is

NormalDistribution[0, σ/Sqrt[n]]

The joint distribution of $\bar{y}$ and, say, $y_1$ is given by

jointDistribution = TransformedDistribution[{y[1], Sum[y[i], {i, n}]/n}, 
  Table[y[i] \[Distributed] NormalDistribution[0, \[Sigma]], {i, n}]]
{#, jointDistribution /. n -> #} & /@ Range[2, 10] // TableForm

$$ \begin{array}{cc} 2 & \text{MultinormalDistribution}\left[\{0,0\},\left( \begin{array}{cc} \sigma ^2 & \frac{\sigma ^2}{2} \\ \frac{\sigma ^2}{2} & \frac{\sigma ^2}{2} \\ \end{array} \right)\right] \\ 3 & \text{MultinormalDistribution}\left[\{0,0\},\left( \begin{array}{cc} \sigma ^2 & \frac{\sigma ^2}{3} \\ \frac{\sigma ^2}{3} & \frac{\sigma ^2}{3} \\ \end{array} \right)\right] \\ 4 & \text{MultinormalDistribution}\left[\{0,0\},\left( \begin{array}{cc} \sigma ^2 & \frac{\sigma ^2}{4} \\ \frac{\sigma ^2}{4} & \frac{\sigma ^2}{4} \\ \end{array} \right)\right] \\ 5 & \text{MultinormalDistribution}\left[\{0,0\},\left( \begin{array}{cc} \sigma ^2 & \frac{\sigma ^2}{5} \\ \frac{\sigma ^2}{5} & \frac{\sigma ^2}{5} \\ \end{array} \right)\right] \\ 6 & \text{MultinormalDistribution}\left[\{0,0\},\left( \begin{array}{cc} \sigma ^2 & \frac{\sigma ^2}{6} \\ \frac{\sigma ^2}{6} & \frac{\sigma ^2}{6} \\ \end{array} \right)\right] \\ 7 & \text{MultinormalDistribution}\left[\{0,0\},\left( \begin{array}{cc} \sigma ^2 & \frac{\sigma ^2}{7} \\ \frac{\sigma ^2}{7} & \frac{\sigma ^2}{7} \\ \end{array} \right)\right] \\ 8 & \text{MultinormalDistribution}\left[\{0,0\},\left( \begin{array}{cc} \sigma ^2 & \frac{\sigma ^2}{8} \\ \frac{\sigma ^2}{8} & \frac{\sigma ^2}{8} \\ \end{array} \right)\right] \\ 9 & \text{MultinormalDistribution}\left[\{0,0\},\left( \begin{array}{cc} \sigma ^2 & \frac{\sigma ^2}{9} \\ \frac{\sigma ^2}{9} & \frac{\sigma ^2}{9} \\ \end{array} \right)\right] \\ 10 & \text{MultinormalDistribution}\left[\{0,0\},\left( \begin{array}{cc} \sigma ^2 & \frac{\sigma ^2}{10} \\ \frac{\sigma ^2}{10} & \frac{\sigma ^2}{10} \\ \end{array} \right)\right] \\ \end{array} $$

So the general distribution is a multivariate normal

MultinormalDistribution[{0, 0}, {{σ^2, σ^2/n}, {σ^2/n, σ^2/n}}]

The general form of the joint density function can then be found with

FullSimplify[PDF[MultinormalDistribution[{0, 0}, {{σ^2, σ^2/n}, {σ^2/n, σ^2/n}}], {y, ybar}],
  Assumptions -> {σ > 0, n > 1}]

$$\frac{n e^{-\frac{n \left(n \text{ybar}^2+y^2-2 y \text{ybar}\right)}{2 (n-1) \sigma ^2}}}{2 \pi \sqrt{n-1} \sigma ^2}$$