ARMA Model constant estimation in statsmodels library

Question

I have a question concerning the statsmodels library, I am just trying to play a bit with that library in order to understand how it work.

from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.tsa.ar_model import ar_select_order
from matplotlib import pyplot as plt

X = ArmaProcess([1, -0.75, 0.15, -0.4, 0.42])
sample_ = X.generate_sample(nsample=800, scale=0.1) # create dummy sample
select_order = ar_select_order(sample_, 10, old_names=False, trend='c')
select_order.model.fit().summary() # fit the model

I get the following results (signs of estimated coefficient are the opposite sign of the dummy model used to generate the data)

but when I try to look at the predicted value from the model I clearly do not get what I am expected

fit = select_order.model.fit()
plt.plot(fit.predict(1, 900))

aka I can see well the auto regressive feature for the first 800 points (the one simulated) but not for the points from 800 to 900, I am probably doing something wrong like not indicating the index etc, but if someone has idea or indication I would definitely be grateful. Please note I am a real novice with this package. In Advance thank you very much

Answer 1

As I already explained in the comment: I think your confusion arises from the difference between sampling/simulating and predicting. For example, doing a 2-step prediction for an $AR(p)$-process, $$X_{n+1}= \sum_{h \geq 1} \theta_h X_{n+1-h} + \epsilon_{n+1} \\ X_{n+2}= \sum_{h \geq 1} \theta_h X_{n+2-h} + \epsilon_{n+2} = \sum_{h \geq 2} \theta_h X_{n+1-h} + \theta_1 \left( \sum_{h \geq 1} \theta_h X_{n+1-h} + \epsilon_{n+1} \right) + \epsilon_{n+2}$$

What is the best linear prediction? Since we are dealing with linear time series models, answering this essentially amounts to setting all the future $\epsilon_t$ equal to zero, i.e. $$\widehat{X_{n+1}} = \sum_{h \geq 1} \theta_h X_{n+1-h} \\ \widehat{X_{n+2}}= \sum_{h \geq 2} \theta_h X_{n+1-h} + \theta_1 \left( \sum_{h \geq 1} \theta_h X_{n+1-h} \right) $$ and so on. Eventually, this will become zero (which is also expected, because with the ACF $\rho(h) = \operatorname{cov}[X_n, X_{n+h}] \to 0$ for $h \to \infty$, this means that $X_{n+h}$ should become independent of the time series up to time $n$, and hence your best guess would be the mean [i.e. a constant]).

Unfortunately I am not particularly familiar with the statsmodels API (I find the forecast package in R much nicer to use, statsmodels seems pretty messy). There seems to be a simulate() method, but apparently only for the classes under statsmodels.tsa.statespace (maybe someone can confirm this or prove me wrong). Anyway, I tried to code this by hand, building on your code, just to illustrate the difference:

values = fit.fittedvalues
n = values.shape[0]
h = 100    # number of forecast steps
values = np.concatenate((values, np.zeros(h)), axis=0)  # extend time series with zeroes
nlags = len(fit.ar_lags)

# now predict stepwise
for i in range(0, h):
    lagged_values = np.concatenate((np.array([1.]), values[(n+i-nlags):(n+i)]), axis=0)
    values[n+i] = -(fit.params @ lagged_values) + np.random.normal(0, np.sqrt(fit.sigma2), 1)
    
plt.plot(values)
plt.plot(fit.predict(1, n+h))

(orange: samples, blue: prediction)

PS: I share your confusion that the coefficients seem to be the negatives of the ones used for generating the samples. There are two different conventions in the literature, and apparently this is inconsistent accross statsmodels. I am really not fond of this package...

[ Just noticed: The [1.] should then probably be a [-1.]... Doesn't really matter since the intercept is so small. ;) ]