# Sufficient statistic for Gaussian $AR(1)$

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.

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})). }$$

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)$$