This isn't a complete answer, but I will give you a starting point to try to tackle part 2.

First notice that $\mathbf{x}_0^T (X^TX)^{-1}\mathbf{x}_0$ is a [change of variables](https://en.wikipedia.org/wiki/Probability_density_function#Function_of_random_variables_and_change_of_variables_in_the_probability_density_function) $\nu:\mathbb{R}^{n \times p} \mapsto \mathbb{R}$ which you can even more simply think of as mapping $V : \mathbb{R}^k \mapsto \mathbb{R}$ via a [pairing function](https://en.wikipedia.org/wiki/Pairing_function) on the index set of the matrix $(X^TX)^{-1}$.

This vector-to-scalar transformation will have a joint density $f_{Y, X}(y, \vec{x}) = F_{X}(\vec{x}) \delta (y - V(\vec{x}))$ where $\delta$ is the [Dirac delta distribution](https://en.wikipedia.org/wiki/Dirac_delta_function). And following this [post](https://stats.stackexchange.com/questions/48304/how-to-find-marginal-distribution-from-joint-distribution-with-multi-variable-de), integrate the joint density with respect to $d\vec{x}$.