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 $\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 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. And following this post, integrate with respect to $d\vec{x}$.

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 $\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 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. And following this post, integrate the joint density with respect to $d\vec{x}$.

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 $\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 on the index set of the matrix.

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. And following this post, integrate with respect to $d\vec{x}$.

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 $\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 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. And following this post, integrate with respect to $d\vec{x}$.