Application of Maximum Likelihood estimation (MLE) to the step of Feasible Generalized Least Square (FGLS) I have the following regression
$$y = X\beta +u$$
where $y$ and $u$ are $(n\times 1)$ and $X$ is a fixed $(n \times k)$ matrix with full column rank and $\beta$ is an unknown $(k\times 1)$ vector of parameters.
$E(u)=0$ and $var(u)=\sigma^2 V$ with $(n\times n) $ symmetric and positive definite non-diagonal matrix $V$.
Let $\hat {\beta_{OLS}}= (X'X)^{-1} X'y$
Suppose $$u_t =\phi u_{t-1} +v_t $$ where $v_t$ is a white noise process with variance $\sigma_v^2$ and $\phi$ is an unknown parameter with $|\phi|<1$.
Describe how to get the feasible generalized least square (FGLS) estimator for $\beta$.
In the FGLS estimation procedure, we should use the maximized likelihood estimator $\hat {\phi_{MLE}}$, which should be obtained through the conditional log-likelihood  function set conditional on the first observation.

What I did
I obtain the variance-covariance matrix as follows:

Since $\phi$ is unknown, the matrix $V$ is also unknown. Therefore, we need to use FGLS as follows:
Step 1: Estimate the generalized regression $y=X\beta + u$ by OLS and obtain $\hat {u}$.
Step 2: Construct the auxiliary regression $\hat {u_t}= \phi \hat {u_{t-1}}+ v_t$ in order to obtain $\hat {\phi_{MLE}}$ and $\hat {v_{MLE}}$.
For that, I write down a likelihood function
$$v_t \sim N(0, \sigma_v^2)$$
$\hat {u_t} \sim N(0, \sigma_v^2)$
$I(\theta, u)=- \frac{n}{2} log(2\pi) - \frac{n}{2} log(\sigma_v^2) - \frac{n}{2\sigma_v^2} (\hat {u_t}- \phi \hat {u_{t-1}})'(\hat {u_t} - \phi \hat {u_{t-1}}) $
where $\theta = [\phi \sigma_v^2]$
FOC:
$$\frac{\partial I(\theta, u)}{\partial \phi} = 0$$
$$\hat {\phi_{MLE}} = (\hat {u_{t-1}}' \hat {u_{t-1}})\hat {u_{t-1}}' \hat {u_{t}}$$
Step 3:
By using $\hat {\phi_{MLE}}$ and $\hat {v_{MLE}}$, I will obtain the matrix $V$ and calculate $\hat {\beta_{FGLS}} = (X'V^{-1}X)^{-1}X'V^{-1}y$

I couldn’t do the step 2. The question states the conditional log-likelihood function, however I cannot write down this log-likelihood function as conditional function to solve for $\hat {\phi_{MLE}}$ and $\hat {v_{MLE}}$.
Please help me how to calculate the step 2. Thank you for your all kinds of helps and suggestions.
 A: Note that $\hat{u}_t = \phi \hat{u}_{t-1} + v_t,\, v_t \sim \mathcal N\left(0, \sigma^2_{v}\right)$ implies $\hat{u}_t | \hat{u}_{t-1}\sim \mathcal N\left(\phi \hat{u}_{t-1}, \sigma^2_{v}\right)$ for $t = 2, \ldots, n$.
With the additional assumption $\hat{u}_1 \sim  \mathcal N\left(0, \sigma^2_{v}/\left(1 - \phi^2\right)\right)$ you get the full log-likelihood kernel
$$
l\left(\phi,\sigma^2_v;\hat{u}\right) = -\frac{n}{2} \cdot \ln\left(\sigma^2_v\right) + \frac{1}{2} \cdot \ln\left(1-\phi^2\right) -
\frac{1}{2\sigma^2_v} \left[\sum_{t=2}^n\left(\hat{u}_t-\phi\hat{u}_{t-1}\right)^2 + \left(1 - \phi^2\right)\hat{u}_1^2 \right].
$$
There is no closed-form solution for the maximum likelihood estimates $\hat{\phi}$ and $\hat{\sigma}^2_v$ in this case, but they can be determined iteratively.
Alternatively, you could treat $\hat{u}_1$ as deterministic/known and just maximize the conditional likelihood by maximizing the conditional log-likelihood kernel
$$
l_c\left(\phi,\sigma^2_v;\hat{u}\right) = -\frac{n-1}{2} \cdot \ln\left(\sigma^2_v\right) -
\frac{1}{2\sigma^2_v} \sum_{t=2}^n\left(\hat{u}_t-\phi\hat{u}_{t-1}\right)^2. 
$$
The corresponding score equations
$$
\begin{align}
\frac{\partial l_c}{\partial \phi}\left(\phi,\sigma^2_v;\hat{u}\right) =
\frac{1}{\sigma^2_v} \sum_{t=2}^n\left(\hat{u}_t-\phi\hat{u}_{t-1}\right) \cdot \hat{u}_{t-1} \overset{!}{=} 0 \\
\iff \sum_{t=2}^n\left(\hat{u}_t-\phi\hat{u}_{t-1}\right) \cdot \hat{u}_{t-1} = 0
\end{align}
$$
and
$$
\begin{align}
\frac{\partial l_c}{\partial \sigma^2_v}\left(\phi,\sigma^2_v;\hat{u}\right) =
-\frac{n-1}{2\sigma^2_v} + 
\frac{1}{2\left(\sigma^2_v\right)^2} \sum_{t=2}^n\left(\hat{u}_t-\phi\hat{u}_{t-1}\right)^2 \overset{!}{=} 0 \\
\iff -\left(n - 1\right) + 
\frac{1}{\sigma^2_v} \sum_{t=2}^n\left(\hat{u}_t-\phi\hat{u}_{t-1}\right)^2  = 0
\end{align}
$$
then yield the maximum conditional likelihood estimates
$$
\begin{align}
&\hat{\phi} = \frac{\sum_{t=2}^n\hat{u}_t\hat{u}_{t-1}}{\sum_{t=2}^n\hat{u}_{t-1}^2}, \\
&\hat{\sigma}^2_v = \frac{1}{n-1} \sum_{t=2}^n\left(\hat{u}_t-\hat{\phi}\hat{u}_{t-1}\right)^2, 
\end{align}
$$
from which you only need $\hat{\phi}$ to calculate $\hat{\beta}_\text{FGLS}$ since the matrix $V$ does not depend on $\sigma^2_v$.

Instead of the two-step FGLS procedure you could also use a (full) maximum likelihood approach with the original data that, under the assumption of a multivariate normal $u$, computes $\left(\hat{\beta}, \hat{\sigma}^2, \hat{\phi}\right)$ as maximizer of the likelihood
$$
\mathcal{L}\left(\beta, \sigma^2, \phi; y, X\right) = \left(2\pi\sigma^2\right)^{-n/2}|V|^{-1/2} \exp\left(-\frac{1}{2\sigma^2}\left(y-X\beta\right)^\top V^{-1} \left(y-X\beta\right)\right),
$$
where
$$
\begin{align}
&V = \frac{1}{1-\phi^2}
\begin{pmatrix}
1 & \phi & \phi^2 & \cdots & \phi^{n-1}\\
\phi & 1 & \phi & \cdots & \phi^{n-2}\\
\phi^2 & \phi & 1 & \cdots & \phi^{n-3}\\
\vdots & \vdots & \vdots & \ddots  & \vdots\\
\phi^{n-1} & \phi^{n-2} & \phi^{n-3}& \cdots & 1
\end{pmatrix},|V|=\left(1-\phi^2\right)^{-1},\\
&V^{-1} =
\begin{pmatrix}
1 & -\phi & 0 & \cdots & 0 & 0\\
-\phi & 1+\phi^2 & -\phi & \cdots & 0 & 0\\
0 & -\phi & 1+\phi^2 & \cdots & 0 & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\
0 & 0 & 0 & \cdots & 1+\phi^2 & -\phi\\
0 & 0 & 0 & \cdots & -\phi & 1
\end{pmatrix}.
\end{align}
$$
Usually, this would be done by iteratively computing the maximizer of the log-likelihood kernel
$$
l\left(\beta, \sigma^2, \phi; y, X\right) = -\frac{n}{2} \cdot \ln\left(\sigma^2\right) + \frac{1}{2} \cdot \ln\left(1-\phi^2\right) -
\frac{1}{2\sigma^2}\left(y-X\beta\right)^\top V^{-1} \left(y-X\beta\right).
$$
A: It would appear that Feasible Generalized Least Square (FGLS) should be used when the covariance matrix of the errors $cov(u) = \sigma^2V$ has a completely unknown form. However, in this problem, you have very specific prior knowledge of its form: it's the autocorrelation matrix of a 1st order autoregressive process, AR(1), which is a Toeplitz matrix. You even wrote this matrix explicitly.
Therefore, I feel the most appropriate and exact solution is to "directly" maximize the conditional log likelihood (conditioned on $\phi$), equivalently minimize the negative log likelihood conditioned on $\phi$, and then maximize over $\phi$:
For each specific value of $\phi$, the log likelihood of $u$ conditioned on $\phi$ is a multivariate Gaussian of $u$ ($u_{\phi} = y - X\beta_{\phi}$), but it's necessary to keep the $det(V_{\phi})$ term, as explained below. For a specific value of $\phi$, i.e., conditioned on $\phi$, this is maximized by the Generalized Least Squares (GLS): $\beta_{\phi} = (X'V^{-1}_{\phi}X)^{-1}X'V^{-1}_{\phi}y$.
But in order to maximize the overall unconditional (log) likelihood, it is necessary to compare likelihood values for different values of $\phi$, which requires to keep the $det(V_{\phi})$ term in the expression of the log likelihood.
Since the value of the unknown $\phi$ is in the bounded 1D range $-1 < \phi < 1$, it is a simple matter to maximize the unconditional likelihood over this range by a numerical 1D search (for example by a Fibonacci search for the maximum, etc.), or analytically by computing the $det(V_{\phi})$ as a function of $\phi$ and differentiating by $\phi$. (The analytical expression for $det(V_{\phi})$ as a function of $\phi$ can be obtained "easily", since $V$ has a very specific "easy" Toeplitz form.)
