Matrix representation of the OLS of an AR(1) process, Is there any precise way to express the OLS estimator of the centred error terms $\{u_t\}
_{t=1}^{n}$ that follows an AR(1) process? In other words, for
\begin{equation}
u_t=\rho u_{t-1}+\varepsilon_t,\quad \varepsilon_t\sim N(0,\sigma^2)
\end{equation}
is there a matrix representation for
\begin{equation}
\hat{\rho}=\frac{(1/n)\sum\limits_{t=1}^{n}u_tu_{t-1}}{(1/n)\sum\limits_{t=1}^{n}u_{t-1}^2}
\end{equation}
? I suspect there should be. However, I seem to fail to find it in Hamilton or other sources or derive an elegant expression myself.
Much appreciated in advance

To facilitate our analysis, we will use the following $(n-1) \times n$ matrices:
$$\mathbf{M}_0 \equiv \begin{bmatrix}
1 & 0 & 0 & \cdots & 0 & 0 \\
0 & 1 & 0 & \cdots & 0 & 0 \\
0 & 0 & 1 & \cdots & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 1 & 0 \\
\end{bmatrix}
\quad \quad \quad 
\mathbf{M}_1 \equiv \begin{bmatrix}
0 & 1 & 0 & \cdots & 0 & 0 \\
0 & 0 & 1 & \cdots & 0 & 0 \\
0 & 0 & 0 & \cdots & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 0 & 1 \\
\end{bmatrix},$$
and the following $n \times n$ matrices:
$$\begin{align}
\mathbf{G}_0 
&\equiv \mathbf{M}_0^\text{T} \mathbf{M}_1
= \begin{bmatrix}
0 & 1 & 0 & \cdots & 0 & 0 \\
0 & 0 & 1 & \cdots & 0 & 0 \\
0 & 0 & 0 & \cdots & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 0 & 1 \\
0 & 0 & 0 & \cdots & 0 & 0 \\
\end{bmatrix} \\[20pt]
\mathbf{G}_1 
&\equiv \mathbf{M}_0^\text{T} \mathbf{M}_0
= \begin{bmatrix}
1 & 0 & 0 & \cdots & 0 & 0 \\
0 & 1 & 0 & \cdots & 0 & 0 \\
0 & 0 & 1 & \cdots & 0 & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & \cdots & 1 & 0 \\
0 & 0 & 0 & \cdots & 0 & 0 \\
\end{bmatrix}.
\end{align}$$
Given the observable time-series vector $\mathbf{u} = (u_1,...,u_n)$ we can then write the model in matrix form as:
$$\mathbf{M}_1 \mathbf{u} = \rho \mathbf{M}_0 \mathbf{u} + \sigma \boldsymbol{\varepsilon}
\quad \quad \quad \quad \quad 
\boldsymbol{\varepsilon} \sim \text{N}(\mathbf{0}, \mathbf{I}).$$
The OLS estimator for the parameter $\rho$ is:
$$\begin{align}
\hat{\rho}_\text{OLS} 
&= (\mathbf{u}^\text{T} \mathbf{M}_0^\text{T} \mathbf{M}_0 \mathbf{u} )^{-1} (\mathbf{u}^\text{T} \mathbf{M}_0^\text{T} \mathbf{M}_1 \mathbf{u} ) \\[12pt]
&= (\mathbf{u}^\text{T} \mathbf{G}_1 \mathbf{u} )^{-1} (\mathbf{u}^\text{T} \mathbf{G}_0 \mathbf{u} ) \\[12pt]
&= \frac{\mathbf{u}^\text{T} \mathbf{G}_0 \mathbf{u}}{\mathbf{u}^\text{T} \mathbf{G}_1 \mathbf{u}} \\[12pt]
&= \frac{\sum_{i=1}^{n-1} u_i u_{i+1}}{\sum_{i=1}^{n-1} u_i^2 }.
\end{align}$$
Note that the OLS estimator for an auto-regressive process is not equivalent to the MLE, since the log-likelihood contains a log-determinant term that is a function of the auto-regression parameter.  The MLE can be obtained via iterative methods if desired.
