How to find the transformed density and log likelihood for this family of distributions? Let's consider the family of transformations given by
$$g_a(Y)=\begin{cases}
\frac{e^{aY}-1}{a} & \text{ for } a\neq 0 \\
Y & \text{ for } a=0 
\end{cases}$$
for $Y\in\mathbb{R}$. Analogous to the estimation of the Box-Cox parameter $\lambda$, the parameter $a\in\mathbb{R}$ can be estimated using a profile likelihood approach. Let $Y_1,...,Y_n\in\mathbb{R}$ be independent
responses together with corresponding predictors $\mathbf{x}_1,...,\mathbf{x}_n\in\mathbb{R}^n$. We assume for $a$ that there exist $\textbf{b}\in\mathbb{R}^p$ and $\sigma^2>0$ such that $g_a(Y_i)\sim N(\textbf{x}_i^T \textbf{b},\sigma^2)$ for $i=1,...,n$.
Here is my questions:

*

*Can we derive such a density function $f_Y$ of the untransformed observations $Y_i$. If so, would it be
$$f_{Y_i}(\textbf{x})=\frac{1}{\sqrt{2\pi \sigma^2}}e^{-\frac{1}{2}\cdot \frac{(\textbf{x}-\textbf{x}_i^T\textbf{b})^2}{\sigma^2}}$$


*How would we find the log-likelihood function $\ell(a,b,\sigma^2;y_1,...,y_n)$? I know that $\sum_{i=1}^n \log(f_{Y_i}(y_i))$.
Thanks in advance.

Only one calculation is needed: we have to differentiate $g_a.$  Everything else is just substitutions and taking logarithms.
For convenience, let $Z=g_a(Y)$ and write $\mu = \mathbf x_i^\prime \beta.$  You say $Z$ has a Normal distribution, which means it has a density function given by
$$f(z) = \frac{1}{\sigma\sqrt{2\pi}}\, \exp\left(-\frac{1}{2\sigma^2}(z-\mu)^2\right)$$
The transformation $y\to g_a(y) = z$ is everywhere increasing, differentiable, and one-to-one.  For $a\ne 0$ its differential is
$$\mathrm d z = \mathrm d\left(\frac{e^ay-1}{a}\right) = e^{ay}\mathrm dy.$$
Consequently, $Z$ has a continuous distribution supported on the image of $g_a$ with probability element

$$h(y)\mathrm d y = f(z(y))\mathrm dz = \frac{1}{\sigma\sqrt{2\pi}}\, \exp\left(-\frac{1}{2\sigma^2}((e^{ay}-1)/a-\mu)^2\right)\, e^{ay}\mathrm dy .$$

The logarithm of the density $h$ is
$$\log h(y) = -\log(\sigma) - \frac{1}{2}\log (2\pi) - \frac{1}{2\sigma^2}\left(\frac{e^{ay}-1}{a}-\mu\right)^2 + ay.$$
Consequently, the log likelihood of $n$ independent observations according to your model is the sum of the logs of their individual likelihoods,

$$\Lambda(a,\beta,\sigma; \mathbf x, \mathbf y)) = - n\log\sigma - \frac{n}{2}\log(2\pi) - \frac{1}{2\sigma^2}\sum_{i=1}^n \left(\frac{e^{ay_i}-1}{a}-\mathbf x_i^\prime \beta \right)^2 + a\sum_{i=1}^n y_i.$$

When $a=0$, taking the limit of the forgoing as $a\to 0$ reduces each $(e^{ay_i}-1)/a$ to $y_i$ and the $+ay_i$ terms vanish, reproducing the usual likelihood for this Normal regression model.
