I am trying to construct the inverse covariance matrix of an AR(2) process of the form $X_t=\theta_1 X_{t-1}+\theta_2 X_{t-2} + \epsilon_t$ with i.i.d. $\epsilon_i$, $\mathbb{E} \epsilon_i =0$, $\mathbb{E}\epsilon_i\epsilon_j=\sigma^2 \delta_{i,j}$. We assume the process to be stationary.
I have already calculated the variance $Var(X_t)=\frac{(1-\theta_2)\sigma^2}{(1+\theta_2)((1-\theta_2)^2-\theta_1^2)}$ and the autocorrelation coefficients $\rho_1=\frac{\mathbb{E}X_t x_{t-1}}{Var(X_t)}=\frac{\theta_1}{1-\theta_2}$, $\rho_2=\frac{\mathbb{E}X_t x_{t-2}}{Var(X_t)}=\theta_2+\frac{\theta_1^2}{1-\theta_2}$ and $\rho_s=\theta_1 \rho_{s-1}+\theta_2 \rho_{s-2}$ for $s>2$.
How do I now proceed to construct the inverse covariance matrix? I tried to calculate it analytically but it did not work out.
EDIT: For fix $T$ I am looking for the inverse matrix of the covariance matrix $\Phi$ with entries $\Phi_{ij}=Cov(X_i,X_j)$, $i, j = 1,\dots,T$. The entries of this covariance matrix are already given by the autocorrelation coefficients above since the process is assumed to be stationary.
Any help would be very appreciated!
Ok, so this is not the full solution, but I will present how to get the explicit form for $\Phi$, maybe somebody figures out how to invert it. So, according to the above notation, $\Phi_{ij} = \rho_{|i-j|}$, where $\rho_0 = 1$ and $\rho_1$ and $\rho_2$ are specified in the question
We will now solve the recurrence relation $\rho_s = \theta_1 \rho_{s-1} + \theta_2 \rho_{s-2}$ with help of generating functions
Define $G(x) = \sum_{i=1}^{\infty}\rho_ix^i$, then use recurrence relation
$G(x) = \sum_{i=1}^{\infty}\rho_{i+1} x^i = \rho_1 + x \rho_2 + \sum_{i=2}^{\infty}x^i \theta_1\rho_{i} + \sum_{i=2}^{\infty}x^i \theta_2 \rho_{i-1} = \rho_1 + x(\rho_2-\theta_1\rho_1) + x \theta_1 G(x) + x^2 \theta_2 G(x) $
Let $c_1 = \rho_2-\theta_1\rho_1 = \theta_2$. Solving for $G(x)$, we get
$$G(x) = \frac{\rho_1 + x c_1}{1 - \theta_1x-\theta_2x^2} = \frac{\rho_1 + x c_1}{-\theta_2 (x-k_1)(x-k_2)} = \frac{\rho_1 + x c_1}{\theta_2 (k_1 - k_2)}(\frac{1}{x-k_2} - \frac{1}{x-k_1}) = \frac{\rho_1 + x c_1}{\theta_2 (k_1 - k_2)}(\frac{1}{k_1}\frac{1}{1-\frac{x}{k_1}} -\frac{1}{k_2}\frac{1}{1-\frac{x}{k_2}}) $$
Here we have factored the quadratic equation, resulting in solutions $k_{1,2} = \frac{1}{2 \theta_2}(-\theta_1 \pm \sqrt{\theta_1^2 + 4 \theta_2})$
Now we substitute the series $\frac{1}{1-x}=1+x+x^2+...$
$$G(x) = \frac{\rho_1 + x c_1}{\theta_2 (k_1 - k_2)} (\frac{1}{k_1}\sum_{i=0}^{\infty} \frac{x^i}{k_1^i} - \frac{1}{k_2}\sum_{i=0}^{\infty} \frac{x^i}{k_2^i}) = (\rho_1 + c_1 x) \sum_{i=0}^{\infty}a_ix^i $$ Where $$a_i = \frac{1}{\theta_2(k_1-k_2)}(\frac{1}{k_1^{i+1}} - \frac{1}{k_2^{i+1}})$$
Finally, $$\rho_{i+1} = \rho_1 a_i + c_1 a_{i-1} = \rho_1 a_i + \theta_2 a_{i-1}$$
So the next step is to try to simplify the explicit form of $\Phi_{ij}$ to see if it is somehow invertible. I'm sure there probably is a faster way to do this, if one somehow exploits some known symmetry of the matrix, but that's beyond me atm