Joint distribution of AR(1) model Let 
$$y_1 \sim   \mathcal{N}(\phi_0,\,\sigma^2),$$ $$y_t|y_{t-1},\,
\phi_0,\,\phi_1,\,\sigma^2 \sim   \mathcal{N}(\phi_0+\phi_1(y_{t-1} -
 \phi_0),\,\sigma^2),$$ for $t=2,3,\cdots,T$.
I want to find the joint distribution of $y_1,\,y_2,\,\cdots,\,y_T$. I think that it is multivariate Normal with mean a vector of $\phi_0$ and precision matrix 
$Q= \frac{1}{\sigma^2}
\begin{pmatrix}
 (1+\phi_1^2) & -\phi_1 & \cdots & 0 \\
 -\phi_1 & (1+\phi_1^2) & \cdots & 0 \\
 \vdots  & \vdots  & \ddots & \vdots  \\
 0 & 0 & \cdots & 1 
\end{pmatrix}.$
I would like to know if this is true and if there is a proof to read it.
 A: Let us write the joint density as
\begin{equation}
p(y_1,\ldots,y_T) = p(y_1)\,p(y_2\mid y_1) \, p(y_3 \mid y_2,y_1) \ldots \, p(y_T \mid y_{T-1},\ldots y_1).
\end{equation}
Furthermore, since the process is AR(1), the past values influence future values only via the latest value, i.e., we have the Markov property $p(y_t \mid y_1,\ldots,y_{t-1}) = p(y_t \mid y_{t-1})$. Substituting this in the factorization, we get 
\begin{equation}
p(y_1,\ldots,y_T) = p(y_1)\, \prod_{i=2}^T p(y_i \mid y_{i-1}).
\end{equation}
The marginal density of $y_1$ and the required conditional densities were given as assumptions. From now on, we shall ignore multiplicative constants (that are independent of $y$), since they are in the end be determined by the requirement that the joint density integrates to 1. 
\begin{equation}
\propto e^{-\frac{1}{2\sigma^2}(y_1 - \phi_0)^2} \times \prod_{i=2}^T e^{-\frac{1}{2\sigma^2}(y_i - \phi_0 - \phi_1\,(y_{i-1} - \phi_0))^2} = e^{-\frac{1}{2}\,E}
\end{equation}
where
\begin{equation}
E = \frac{1}{\sigma^2}(y_1 - \phi_0)^2 + \sum_{i=2}^T \frac{1}{\sigma^2}(y_i - \phi_0 - \phi_1(y_{i-1} - \phi_0))^2
\end{equation}
\begin{equation}
= \frac{1}{\sigma^2}(y_1 - \phi_0)^2 + \sum_{i=2}^T \frac{1}{\sigma^2}\left((y_i - \phi_0)^2 - 2\,(y_i - \phi_0)\,\phi_1\,(y_{i-1} - \phi_0) + \phi_1^2 (y_{i-1} - \phi_0)^2 \right)
\end{equation}
\begin{equation}
= \sum_{i=1}^{T-1}\frac{1}{\sigma^2}\,(1 + \phi_1^2)(y_i - \phi_0)^2  + \frac{1}{2\sigma^2} (y_T - \phi_0)^2 + \sum_{i=1}^{T-1}\frac{1}{\sigma^2}\,2\,(-\phi)\,(y_{i+1}-\phi_0)\,(y_i-\phi_0). 
\end{equation}
So the joint density is proportional to 
\begin{align}
\mathrm{exp}\bigg(-\frac{1}{2}\,\sum_{i=1}^{T-1}\frac{1}{\sigma^2}\,(1 + \phi_1^2)(y_i - \phi_0)^2  + \frac{1}{2\sigma^2} (y_T - \phi_0)^2 \\+ \sum_{i=1}^{T-1}2\,\frac{1}{\sigma^2}\,(-\phi)\,(y_{i+1}-\phi_0)\,(y_i-\phi_0)\bigg),
\end{align}
Observe that this exponent is a quadratic form of a vector consisting of variables $(y_i - \phi_0)$. Thus we conclude that the joint density is a multivariate normal with means $E(y_i) = \phi_0$ and the precision matrix can be read from the previous expression, since we have $E = (y - \phi_0\mathbf{1})\,\Sigma^{-1}\,(y-\phi_0\,\mathbf{1})$. Namely, 


*

*If $i=j$ and $i<T$,  $\Sigma^{-1}_{ij} = (1 + \phi_1^2) / \sigma^2$

*If $i=j=T$, $\Sigma^{-1}_{ij} = 1 / \sigma^2$

*If $|i-j|=1$, $\Sigma^{-1}_{ij} = -\phi_1 / \sigma^2$

*If $|i-j|>1$, $\Sigma^{-1}_{ij} = 0$,


which is indeed the form of the precision matrix that was claimed in the question.
