Estimate of the intercept is off in a simulated AR(1) model

Question

I've been working with a SARIMAX model for forecasting and found myself struggling to accurately interpret its long-term forecasts. To better understand the underlying mechanics and perhaps pinpoint where my confusion originates, I decided to revisit the fundamentals and examine the behavior of a simpler AR(1) model with a constant term.

I generate data according to the following AR(1) process:

data[i] = c + phi * data[i-1] + eps[i]

where c is a constant, phi is the autoregressive coefficient, and eps[i] is a white noise error term. After generating the data, I first regress the generated series on just a constant and then with its lagged value as an explanatory variable. When regressing on just the constant, the estimate I get is very close to c, as expected. However, things get tricky when I add the lagged data as an explanatory variable.

Despite the mean of the generated data being consistent with the theoretical expectation of $c/(1−phi)$, the constant estimate from regressing $X$ on both a constant and $X_{t-1}$ significantly deviates from the original c used in the data generation process. The coefficient for $X_{t-1}$ seems reasonable, but the constant term is different from c in the DGP.

I'm puzzled by this discrepancy. I initially thought my data generation process (DGP) might be incorrect, but the mean of the data aligns with theoretical expectations, suggesting the DGP is correct.

Questions:

1. Is my expectation to approximate c from my DGP when including the lagged variable as an explanatory variable incorrect?

2. Are there any common pitfalls or considerations in estimating AR(1) model parameters that I might be overlooking?

Any insights or suggestions on what might be causing these unexpected estimates and how to address them would be greatly appreciated!

Answer 1

Here is some python code showing everything working okay:

import numpy as np
from sklearn.linear_model import LinearRegression

eps = np.random.randn(1000)
data = list(np.random.randn(1))

for error in eps:
    data += [3 + 0.5*data[-1] + error]

lr.fit(np.array(data).reshape(-1,1)[:-1,:], np.array(data).reshape(-1,1)[1:,:])

I get lr.coef_ = 0.49958832 and lr.intercept_ = 2.93268167 which seems pretty close?

I am guessing you just had a coding error.