Why is $f_{Y|X}(y|x) = f_\varepsilon(y - g(x))$ for the regression model $Y = g(X) + \varepsilon$? Suppose we have the model
$$
Y = g(X) + \varepsilon,
$$
where the errors are zero-mean and independent of $X$.
I read that the conditional probability density function $f_{Y|X}(y|x)$ of $Y$ can be written as
$$
f_{Y|X}(y|x) = f_\varepsilon(y - g(x)),
$$
where $f_\varepsilon$ is the density of $\varepsilon$. How can this be proven?
Here is my attempt:
$$
\begin{align}
f_{Y|X}(y|x) = f_{Y|X}(Y=y|X=x) &= f_{Y|X}(g(X)+\varepsilon=y|X=x) \\
&= f_{g(X)+\varepsilon|X}(g(x)+\varepsilon=y|X=x) \\
&= f_{g(X)+\varepsilon|X}(\varepsilon = y - g(x)|X=x) \\
\end{align}
$$
Since $\varepsilon$ is independent of $X$ I then drop the $'X=x'$ in the function argument to get
$$
\begin{align}
f_{Y|X}(y|x) &= f_{g(X)+\varepsilon|X}(\varepsilon = y - g(x)).
\end{align}
$$
The problem then is how to get the correct subscript on the density, eg. how to get $f_\varepsilon$ instead of $f_{g(X)+\varepsilon|X}$?
 A: As mentioned in the comments, it follows pretty directly. But if you want to see it by going through a bunch of explicit steps, here's one way to do it.
The conditional distribution of $Y$ given $X$ can be obtained by integrating $\epsilon$ out of the joint distribution of $Y$ and $\epsilon$, given $X$:
$$p(y \mid x) = \int_{-\infty}^\infty p_{Y,\epsilon \vert X}(y, \epsilon \mid x) d\epsilon$$
Using the chain rule of probability, the integrand can be factorized as follows:
$$= \int_{-\infty}^\infty p_{Y \vert X,\epsilon}(y \mid x, \epsilon) p_{\epsilon \vert X}(\epsilon \mid x) d\epsilon$$
$\epsilon$ and $X$ are independent, so the conditional distribution of $\epsilon$ given $X$ is equal to the marginal distribution of $\epsilon$:
$$= \int_{-\infty}^\infty p_{Y \vert X,\epsilon}(y \mid x, \epsilon) p_\epsilon(\epsilon) d\epsilon$$
$Y$ is given by a deterministic function of $X$ and $\epsilon$; if we knew the values of $X$ and $\epsilon$, we would know the value of $Y$ with certainty. We can therefore write the conditional distribution of $Y$ given $X$ and $\epsilon$ as a Dirac delta function $\delta(g(x) + \epsilon - y)$. This places all probability mass at the point where $y = g(x) + \epsilon$, and zero mass at all other points.
$$= \int_{-\infty}^\infty \delta(g(x) + \epsilon - y) p_\epsilon(\epsilon) d\epsilon$$
According to the rules of integration with delta functions, $\int \delta(z-c) g(z) dz = g(c)$. So, our integral simplifies to:
$$= p_\epsilon(y - g(x))$$
