Non-Sensible Estimates for MLE of AR Processes I am taking a course on Time Series Econometrics and I am solving a problem set that requires students to explicitly write maximum likelihood functions for, as an example, AR processes and estimate them by Maximum Likelihood. We are required to code both the exact likelihood and the conditional likelihood.
Right now, consider that one wishes to estimate an AR(2) process of the form

$y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \epsilon_t $, where $\epsilon_t \sim N(0, \sigma^2)$

So, we seek estimates of $(\phi_1, \phi_2)$ and $\sigma^2$. In order to code a function that takes as input the data $y_t$ and parameters $(\phi_1, \phi_2, \sigma^2)$ and outputs the exact log-likelihood, I followed Hamilton's book (Chapter 5). The routine I wrote seems to work fine and indeed finds the estimators. My question is on the theoretical side and follows.
To derive the expressions on Hamilton's book, one has to assume that the DGP being estimated is stationary. Otherwise, the formulas make no sense. But after getting the estimated coefficients, as a sanity check, I computed the roots of the AR polynomial and found that the estimated coefficients imply a non-stationary AR(2). Fearing a code mistake, I implemented the same routine using standard functions from the statsmodels package in Python found that: 1) estimates are similar and 2) the implied AR(2) is non-stationary as well.
I am very confused. On one hand, we are starting with the hypothesis that the original model is stationary. This enables us to derive neat formulas even for the exact likelihood. On the other hand, the estimated model using the exact likelihood derived under stationarity is non-stationary. What should I do? I must have understood something very wrong! Any ideas?
One last comment: I proposed another implementation for the exact likelihood in which my Python function checks, before any calculation, if the inputed values for $(\phi_1, \phi_2)$ imply a stationary AR(2) or not. It they don't, the function outputs $-\infty$. When I optimize this function, I get very different estimates. I recall my professor saying that non-stationary AR(p) models, for example, have an alternate stationary representation. Maybe I am finding this one when I do this "pseudo-likelihood", but I am not even sure this is reasonable. Any ideas on this as well? Thanks!
 A: The expression for the AR(2) log likelihood (eq. 5.3.8 in Hamilton) has a term for the joint density of $y_1$ and $y_2$.  This term indeed only makes sense if the process is stationary as it involves the stationary variance-covariance matrix of $y_1,y_2$.
If you run unconstrained optimisation on the expression you indeed risk ending up with meaningless estimates so you need to impose the constraints
$$
-1<\phi_2<\min(1-\phi_1,1+\phi_1)  \tag{1}
$$
when you do the optimization.  One way to do this is to work with some suitable one-to-one transformation of the parameters, e.g. the partial autocorrelations at lag 1 and 2,
$$
\phi_{11}=\frac{\phi_1}{1-\phi_2}
$$
and
$$
\phi_{22}=\phi_1,
$$
as stationarity is equivalent to these taking values between -1 and 1.  Or you could work with the atanh transforms of these, mapping the parameters to all of $\mathbb{R}^2$.  This is the transformation used by function arima in R when fitting by exact maximum likelihood.
Note that the term $\frac12\log\{(1+\phi_2)^2[(1-\phi_2)^2-\phi_1^2]\}$ in eq. 5.3.8 tends to $-\infty$ when the parameters approach the boundary defined by (1).  So the MLE will almost always be somewhere at the interior of the triangular region defined by (1) when fitting the model based on the exact likelihood.
