Consider the model $$ Y_i^{(\lambda)} = \alpha+\beta x_i + \varepsilon_i,\qquad \varepsilon_i,\ (i=1,\ldots,n) \sim \mathrm{i.i.d.}\ N(0,\sigma^2) $$ where $$ \begin{align} y_i^{(\lambda)} & = \begin{cases} \dfrac{y_i^\lambda-1}{\lambda(\operatorname{GM}(y))^{\lambda -1}} , &\text{if } \lambda \neq 0 \\[12pt] \operatorname{GM}(y)\log{y_i} , &\text{if } \lambda = 0 \end{cases} \\[12pt] \text{and } & \operatorname{GM}(y) = (y_1\cdots y_n)^{1/n}\text{ is the geometric mean.} \end{align} $$
As my estimate of $\lambda$, I use the value that minimizes the sum of square of residuals when $\alpha$ and $\beta$ have likewise been estimated by least squares.
In finding a confidence region for $\alpha$ and $\beta$, how should one take into account the uncertainty in the estimate of the Box–Cox parameter $\lambda$?