Sufficient statistic for Gaussian $AR(1)$

Question

Question Does the Gaussian $AR(1)$ model, with a fixed sample size $T$, have nontrivial sufficient statistics?

The model is given by $$ y_t = \rho y_{t-1}, \, t = 1, \cdots, T, \; \epsilon_i \stackrel{i.i.d.}{\sim} \mathcal{N}(0, \sigma^2), $$ parametrized by $(\rho, \sigma^2)$. Let $y = (y_1, \cdots, y_T)'$. Conditional on $y_0$, the likelihood function is \begin{align*} L(y) &= \frac{1}{ ( \sqrt{2 \pi \sigma^2} )^n } e^{-\frac{1}{2\sigma^2} ( \sum_{t = 1}^T (y_t - \rho y_{t-1})^2)} \\\\ &= \frac{1}{ ( \sqrt{2 \pi \sigma^2} )^n } e^{-\frac{1}{2\sigma^2} y' P(\rho) y}, \end{align*} where $P(\rho)$ is the $T \times T$ tridiagonal Toeplitz matrix given by $$ P(\rho) = \begin{bmatrix} 1 + \rho^2 & -\rho & 0 & 0 & & & \\ -\rho & 1 + \rho^2 & -\rho & 0 & & & \\ 0 &-\rho & 1 + \rho^2 & -\rho & & & \\ & & & & \ddots& & \\ & & & & -\rho & 1 + \rho^2 & -\rho\\ & & & & 0 & -\rho & 1\\ \end{bmatrix}. $$

My guess is that, if a nontrivial sufficient statistic exists, it would involve factorizing $P(\rho)$.

Comment

This situation is a little different with the case of the Gaussian linear model $$ Y = X \beta + \epsilon, \;\; \epsilon \stackrel{d}{\sim} \mathcal{N}(0, \sigma^2 I_T) $$ parametrized by $\beta \in \mathbb{R}^p$ and $\sigma^2 > 0$. The design matrix $X, T \times p$ is considered fixed. The likelihood function is \begin{align*} L(y) &= \frac{1}{ ( \sqrt{2 \pi \sigma^2} )^n } e^{-\frac{1}{2\sigma^2} ( Y - X \beta)'( Y - X \beta) } \\\\ &= \frac{1}{ ( \sqrt{2 \pi \sigma^2} )^n } e^{-\frac{1}{2\sigma^2} \left[ ( Y - X \hat{\beta})'( Y - X \hat{\beta}) + (\hat{\beta} - \beta)' X' X (\hat{\beta} - \beta) \right]}. \end{align*} This makes $(\hat{\beta}, s^2)$ sufficient (and minimal), where $\hat{\beta}$ is the OLS estimate $\hat{\beta}$ and $s^2 = \frac{1}{T-1} ( Y - X \hat{\beta})'( Y - X \hat{\beta})$. But in the $AR(1)$ case, it doesn't make sense to consider the covariates $y_{t-1}$ as fixed.

Answer 1

Turns out the answer is simpler than expected.

Consider the conditional form of the likelihood function is (still conditioning on $y_0$) \begin{align*} L(y) &= \frac{1}{ ( \sqrt{2 \pi \sigma^2} )^n } e^{-\frac{1}{2\sigma^2} ( \sum_{t = 1}^T (y_t - \rho y_{t-1})^2)}. \\\\ \end{align*}

As in the case of linear model $$ \sum_{t = 1}^T (y_t - \rho y_{t-1})^2 = \sum_{t = 1}^T (y_t - \hat{\rho} y_{t-1})^2 + (\hat{\rho} - \rho)^2 \sum_{t = 1}^T y_{t-1}^2 $$ where $\hat{\rho} = \frac{\sum_{t = 1}^T y_{t-1} y_t}{\sum_{t = 1}^T y_{t-1}^2}$ is the OLS estimate, and $\sum_{t = 1}^T (y_t - \hat{\rho} y_{t-1})^2$ is the residual sum of squares.

Evidently, $$ T_1 = (\sum_{t = 1}^T (y_t - \hat{\rho} y_{t-1})^2,\, \hat{\rho},\, \sum_{t = 1}^T y_{t-1}^2) $$ is a minimal sufficient statistic.

The only difference with the linear model case is the additional term $\sum_{t = 1}^T y_{t-1}^2 = X'X$. This is not surprising. Inspecting the calculation for the linear model, we see that taking the design matrix $X$ as fixed is unnecessary. One only needs to take $X'X$ as fixed. In the time series case, the statement "taking $\sum_{t = 1}^T y_{t-1}^2$ as fixed" doesn't make sense. Instead it appears as part of a sufficient statistic.

Equivalently, \begin{align*} \sum_{t = 1}^T (y_t - \rho y_{t-1})^2 &= \sum_{t = 1}^T y_t^2 - 2 \rho \sum_{t = 1}^T y_{t-1} y_t + \rho^2 \sum_{t = 1}^T y_{t-1}^2 \\\\ &= (1 + \rho^2) \sum_{t = 1}^T y_{t-1}^2 + y_T^2 - 2 \rho \sum_{t = 1}^T y_{t-1} y_t \end{align*} means $$ T_2 = (\sum_{t = 1}^T y_{t-1}^2,\, y_T^2,\, \sum_{t = 1}^T y_{t-1} y_t) $$ is another minimal sufficient statistic. Geometrically, it is clear that $T_2$ contains the same information as $$ T_1 = \mbox{ (sum of OLS square residuals, OLS estimate, squared norm of $(y_{t-1})$). } $$

Answer 2

Here is a direct proof.

Suppose the model is of the form $$Y_1=\varepsilon_1\quad ,\quad Y_t=\rho Y_{t-1}+\varepsilon_t \quad,\,t=2,3,\ldots,n$$

where $\varepsilon_1,\ldots,\varepsilon_n$ are i.i.d $N(0,\sigma^2)$.

In other words, we have $$Y_1=\varepsilon_1\quad,\quad Y_t=\sum_{j=0}^{t-1} \rho^j \varepsilon_{t-j} \quad,\,t=2,3,\ldots,n$$

That is, $Y=(Y_1,Y_2,\ldots,Y_n)'$ is a linear transformation of $\varepsilon=(\varepsilon_1,\varepsilon_2,\ldots,\varepsilon_n)'$:

$$ \begin{bmatrix} Y_1 \\ Y_2 \\ Y_3 \\ \vdots \\ Y_n \end{bmatrix} =\underbrace{\begin{bmatrix} 1 & 0 &0 &\cdots & 0 \\ \rho & 1 & 0 & \cdots & 0 \\ \rho^2 & \rho & 1 &\cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & 0 \\ \rho^{n-1} & \rho^{n-2} & \rho^{n-3} & \cdots & 1 \end{bmatrix}}_{L_{\rho}} \begin{bmatrix}\varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \vdots \\ \varepsilon_n \end{bmatrix} $$

Now, $$\varepsilon \sim N_n(0,\sigma^2 I_n)\implies Y=L_{\rho}\,\varepsilon \sim N_n(0,\sigma^2 L_{\rho}{L_{\rho}}')$$

Since $\det L_{\rho}=1$, the matrix $L_{\rho}$ is nonsingular for any $\rho$.

So the density of $Y$ at $y=(y_1,\ldots,y_n)\in \mathbb R^n$ is

$$f_{\rho,\sigma^2}(y)\propto \exp\left\{-\frac1{2\sigma^2} y' (L_{\rho}{L_{\rho}}')^{-1} y'\right\}$$

It is not difficult to verify that

$$(L_\rho {L_\rho}')^{-1}=(L_\rho^{-1})' L_\rho^{-1}= \begin{bmatrix} 1+\rho^2 & -\rho & 0 & 0 &\cdots& 0 \\-\rho & 1+\rho^2 & -\rho & 0 &\cdots & 0 \\ 0 & -\rho & 1+\rho^2 & -\rho & \cdots &0 \\ & & \ddots & \ddots & \ddots & \\ 0 & 0 & \cdots & -\rho & 1+\rho^2 & -\rho \\ 0 & 0 & \cdots & 0 & -\rho & 1 \end{bmatrix} $$

Therefore,

\begin{align} f_{\rho,\sigma^2}(y)&\propto \exp\left\{-\frac1{2\sigma^2}\left[(1+\rho^2)\sum_{t=1}^{n-1} y_t^2-2\rho\sum_{t=1}^{n-1} y_t y_{t+1}+y_n^2\right]\right\} \\&=\exp\left\{-\frac{(1+\rho^2)}{2\sigma^2}\sum_{t=1}^{n-1} y_t^2 + \frac{\rho}{\sigma^2}\sum_{t=1}^{n-1} y_t y_{t+1}-\frac1{2\sigma^2}y_n^2\right\} \end{align}

Indeed a minimal sufficient statistic for $(\rho,\sigma^2)$ is $$T(Y_1,\ldots,Y_n)=\left(\sum_{t=1}^{n-1} Y_t^2,\sum_{t=1}^{n-1} Y_t Y_{t+1},Y_n^2\right)$$