Given the simple linear regression model $$ y_i = \beta_0 + \beta_1 x_i + \epsilon_i$$ where $\beta_0$ and $\beta_1$ are fixed paramters, $x_i$ are nonrandom variables and the errors $\epsilon_i$ are gaussian distributed according to $\epsilon_i \sim N(0,\sigma^2)$.
I know that in this setting the $y_i$ are also gaussian distributed and it holds: $y_i \sim N(\beta_0 + \beta_1 x_i, \sigma^2)$.
It is straigt forward to show that $E[y]=\beta_0 + \beta_1 x_i$ and that $Var[y]=\sigma^2$, but what I do not understand is, how I could show that the $y_i$ are indeed gaussian distributed (and not distributed to some other distribution with the mean $E[y]=\beta_0 + \beta_1 x_i$ and $Var[y]=\sigma^2$).
I mean it makes sense that the responses are also gaussian when all the other terms are nonrandom expect for the gaussian noise but what I try to do is to derive this mathematically, i.e., I tried to plug in a gaussian for $\epsilon_i$ and plug this into the uppermost equation but that doesnt add up to a gaussian with the mean as calculated above.
Is it possible to show mathematically that $y_i$ are indeed gaussian?
