Why does ARIMA not perform well?

Question

I created a simple AR(1) process with a constant=1 and coefficient=0.5:

def ar_1(k,b,n=1000):
  s1=[0]
  for _ in range(n-1):
    s1.append(k + b*s1[-1] + 0.01*np.random.normal() )
  return s1

S = pd.DataFrame(ar_1(1, 0.5))

And then ran an ARIMA model on the series:

armodel = ARIMA(S,order=(1,0,0))
armodel_fit = armodel.fit()
armodel_fit.summary()

and got the following summary:

whereas running a stasmodels linear regression on the same series gives perfect results. The following is the code and the model summary:

y = S.shift(-1).fillna(method='ffill')
X = S
X = sm.add_constant(X)
lin_reg = sm.OLS(y,X)
m = lin_reg.fit()
m.summary()

As you can see, linear regression gives perfectly good results whereas the ARIMA doesn't perform well on the simplest autoregressive process. So why is it that the ARIMA model is giving these results? Am I not interpreting the results correctly? How does the ARIMA model actually work behind the scenes?

Answer 1

TLDR

The reason is because ARIMA class does regression with AR(1) errors when a constant is present, not the AR(1) model that you expect and created the series for. ARIMA class estimates AR(1) as you expect only when the constant is zero, i.e. unconditional mean is zero. I mean statsmodels v0.12.1.

Theory

The AR(1) that OP generated the series for is: $$x_t=c+\phi x_{t-1}+\varepsilon_t$$ The model that is being estimated by the code OP invoked is a different one, and is called regression with AR(1) errors. However, being Python developers, the authors of statsmodels package didn't care for conventions, and still call it ARIMA. This is what they're estimating: $$x_t=\mu+u_t$$ where $u_t=\varphi u_t+\varepsilon_t$

Here's how I found out about it. This is from statsmodels.tsa.arima.model.ARIMA doc: This model incorporates both exogenous regressors and trend components through “regression with ARIMA errors”. and the trend parameter: Default is ‘c’ for models without integration, and no trend for models with integration.

Proof of red herring

The models surely look similar. In fact, we have $E[x_t]=\mu=\frac c {1-\phi}$. Hence, if the statsmodels ARIMA could fit correctly to OP's data, it would get the following: $\hat\phi\approx 0.5$ and $\hat\mu\approx 2$. Needless to say, a properly estimated AR(1) process should have rendered $\hat c\approx 1$ and $\hat\varphi\approx 0.5$, i.e. what OP expected to see.

However, as OP has shown it got $\mu\approx 1.3$ and $\phi\approx 1$. So, this threw some people off the trail to the correct answer to OP's question. They thought this was the issue. They also found an excuse for Python's failure: the initial point was "too far" from the unconditional mean yada-yada.

Of course, they were right to point to a problem. Here's chart of the series and the initial point is far from the rest, indeed.

So, let's see how it's resolved. Let's fit the model after cutting out first 50 observations.

Here's the output:

SARIMAX Results                                
==============================================================================
Dep. Variable:                      y   No. Observations:                  950
Model:                 ARIMA(1, 0, 0)   Log Likelihood                3047.631
Date:                Wed, 13 Jan 2021   AIC                          -6089.262
Time:                        23:23:05   BIC                          -6074.693
Sample:                             0   HQIC                         -6083.711
                                - 950                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.9989      0.001   3267.676      0.000       1.998       2.000
ar.L1          0.4816      0.029     16.677      0.000       0.425       0.538
sigma2      9.561e-05   4.41e-06     21.673      0.000     8.7e-05       0.000
===================================================================================

Look at the estimated parameters. They match egression with AR(1) process, but don't match AR(1) process that OP created: the AR coefficient $\approx 0.5$ but the constant is not 1, it is $\approx 2$, i.e. close to what $\mu$ should be not $c$.

Can we overcome estimation issue?

It turns out that if OP used the specific fitting method, then ARIMA class would find the correct parameters even without cutting the head of the series. Here's example:

armodel = ARIMA(S,order=(1,0,0))
armodel_fit = armodel.fit(method='hannan_rissanen')
print(armodel_fit.summary())

Which generates the output:

SARIMAX Results                                
==============================================================================
Dep. Variable:                      y   No. Observations:                 1000
Model:                 ARIMA(1, 0, 0)   Log Likelihood              -11077.953
Date:                Thu, 14 Jan 2021   AIC                          22161.906
Time:                        02:57:38   BIC                          22176.629
Sample:                             0   HQIC                         22167.502
                               - 1000                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          1.9932      0.001   3132.867      0.000       1.992       1.994
ar.L1          0.4909      0.010     48.699      0.000       0.471       0.511
sigma2         0.0001   1.38e-06     76.655      0.000       0.000       0.000

No, these are not AR(1) parameters $c,\phi$, of course, but these are the best regression with AR(1) parameters $\mu,\varphi$ that fit the series. Again, this doesn't solve the problem that ARIMA is not really ARIMA, but it does solve the estimation issue.

AR(1) done properly in Python

Now that we understand the problem, let us consider solutions. The most straightforward solution is to not use misleading and weak ARIMA class, but instead use a dedicated class for AR(p) models. Here's how to estimate AR(p) with Python:

from statsmodels.tsa.ar_model import AutoReg
 res = AutoReg(ar_1(1,0.5), lags = [1]).fit()
 res.summary()

The output as expected:

AutoReg Model Results
Dep. Variable:  y   No. Observations:   1000
Model:  AutoReg(1)  Log Likelihood  3175.129
Method: Conditional MLE S.D. of innovations 0.010
Date:   Mon, 11 Jan 2021    AIC -9.188
Time:   17:47:33    BIC -9.174
Sample: 1   HQIC    -9.183
1000        
coef    std err z   P>|z|   [0.025  0.975]
intercept   1.0134  0.009   117.039 0.000   0.996   1.030
y.L1    0.4935  0.004   113.903 0.000   0.485   0.502

Note, we didn't have to cut the head out of the series. It simply estimated the process without any fuss.

ARIMA class in statsmodels works only without constant

Can ARIMA class do anything right? Here's example of how ARIMA class correctly estimates AR(1) with zero mean, i.e. in my notation $c=0$:

armodel = ARIMA(ar_1(0,0.5),order=(1,0,0))
armodel_fit = armodel.fit()
armodel_fit.summary()

SARIMAX Results
Dep. Variable:  y   No. Observations:   1000
Model:  ARIMA(1, 0, 0)  Log Likelihood  3205.138
Date:   Mon, 11 Jan 2021    AIC -6404.276
Time:   18:42:41    BIC -6389.553
Sample: 0   HQIC    -6398.680
- 1000      
Covariance Type:    opg     
coef    std err z   P>|z|   [0.025  0.975]
const   6.197e-05   0.001   0.104   0.917   -0.001  0.001
ar.L1   0.4786  0.028   17.043  0.000   0.424   0.534
sigma2  9.626e-05   4.45e-06    21.638  0.000   8.75e-05    0.000

What about SARIMAX class in statsmodels?

Let's add this to confusion: although ARIMA class uses SARIMAX class under the hood, it defines the trend parameter for the constant differently!

In SARIMAX the constant with a trend parameter is defined as you'd expect in AR(1) model, and not like in regression with AR(1) errors process. Yet, although the process is correct, the estimation routine still sucks:

from statsmodels.tsa.statespace.sarimax import SARIMAX
model = SARIMAX(S, order=(1,0,0), trend='c')
res = model.fit()
res.summary()

Output:

SARIMAX Results
Dep. Variable:  y   No. Observations:   1000
Model:  SARIMAX(1, 0, 0)    Log Likelihood  1827.315
Date:   Tue, 12 Jan 2021    AIC -3648.629
Time:   03:15:15    BIC -3633.906
Sample: 0   HQIC    -3643.033
- 1000      
Covariance Type:    opg     
coef    std err z   P>|z|   [0.025  0.975]
intercept   0.0526  0.003   19.215  0.000   0.047   0.058
ar.L1   0.9721  0.003   334.121 0.000   0.966   0.978
sigma2  0.0015  9.16e-06    162.577 0.000   0.001   0.002
Ljung-Box (L1) (Q): 9.08    Jarque-Bera (JB):   5992619.42
Prob(Q):    0.00    Prob(JB):   0.00
Heteroskedasticity (H): 0.03    Skew:   14.59
Prob(H) (two-sided):    0.00    Kurtosis:   381.12

It fails for the reason correctly pointed by @cfulton and it shouldn't have, if it was written properly.

Here's how SARIMAX fits properly with the head of the series cut off:

SARIMAX Results
Dep. Variable:  y   No. Observations:   950
Model:  SARIMAX(1, 0, 0)    Log Likelihood  3047.631
Date:   Thu, 14 Jan 2021    AIC -6089.262
Time:   00:12:00    BIC -6074.693
Sample: 0   HQIC    -6083.711
- 950       
Covariance Type:    opg     
coef    std err z   P>|z|   [0.025  0.975]
intercept   1.0362  0.058   17.943  0.000   0.923   1.149
ar.L1   0.4816  0.029   16.674  0.000   0.425   0.538
sigma2  9.563e-05   4.41e-06    21.668  0.000   8.7e-05 0.000

Notice, how unlike ARIMA class this one estimates the constant to be around 1, i.e. close to AR(1) constant $c$, and not the unconditional mean $\mu$ as in regression with ARIMA.

Conclusion

This is all because Python statsmodels ARIMA related packages are a mess. They can't even follow their own made up conventions consistently.

Next, I'm planning to show you how, when implemented properly, the state space representation of AR(1) process should estimate OP's series, without the need to cut its head out. Both ARIMA ans SARIMAX classes use statsmodels' state space model utilities. It's a lame excuse to blame the initial point for a failure of estimation.

Answer 2

What is going wrong with ARIMA?

The reason is, as @chris-haug suggested, that the large distance between the initial value and the mean of the process at the beginning of the series make it appear unlikely (in the sense of the likelihood function) that the series is stationary (in the sense of all the values being drawn from the same conditional distribution).

It does not have to do with whether the model is parameterized with the constant in mean form or in intercept form (see the end of this answer for more details on that).

Technically, the reason is that ARIMA model assumes that, in stationary models, the first observation is drawn from the unconditional distribution implied by the parameters. The unconditional distribution implied by your generating parameters is $N(c / (1 - \phi), \sigma^2 / (1 - \phi^2))$. In your model, $c = 1, \phi = 0.5,$ and $\sigma^2 = 0.0001$, so the unconditional distribution is $N(2, 0.00013)$.

Recall that the rule of thumb for Gaussian distributions is that 99.7% of the mass of the distribution falls within 3 standard deviations of the mean. Here, the standard deviation is $\sqrt{0.00013} = 0.0115$ so here that means that 99.7% of the mass of this distribution lies in the range $[2 - 0.0345, 2 + 0.0345] = [1.9655, 2.0345]$. In summary, it is incredibly unlikely that the value 0 would arise from this unconditional distribution.

So, when the model tries out your parameters, it finds that they are extremely unlikely to fit the data. We can check this using the SARIMAX model:

model = sm.tsa.SARIMAX(S, order=(1,0,0), trend='c')
print(model.loglikeobs([1., 0.5, 0.01**2]))

yields:

[-14996.45760938      3.5751061       2.9770201       2.65990702 ...

The first observation, S[0] = 0 is just incredibly unlikely.

Now, you might think that given our discussion above, the second observation, S[1] = 1.004714 is still extremely unlikely, so why is its likelihood so much higher? The answer is that conditional on the first value being zero, then observing the next value to be 1.004714 is actually right in line with the parameters of the model. This is because, given these parameters, the model would have predicted $y_2 = 1 + 0.5 y_1 = 1$.

Why is OLS doing better?

This leads directly into the answer of why running OLS does a better job. That's because OLS works exactly by conditioning on the first observation. So OLS is not evaluating how good or bad the first observation is at all, it's simply saying: conditional on the first observation's value, what are the parameters that minimize the sum of squared errors for the remaining observations.

Can we understand what ARIMA is doing here?

Basically, it appears to the model that the first few observations were not drawn from the same conditional distribution as the remaining observations. Time series that have this property are called non-stationary, and one characteristic of them is that they have a unit root - which is to say, they have $\phi = 1$.

Fitting the model we constructed earlier, we can see that the model is trying to fit the data by estimating non-stationary parameters (actually, these models are constrained to be stationary, so instead it is estimating parameters that are almost non-stationary):

results = model.fit()
print(results.summary())

yields:

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                      0   No. Observations:                 1000
Model:               SARIMAX(1, 0, 0)   Log Likelihood                1827.942
Date:                Mon, 11 Jan 2021   AIC                          -3649.884
Time:                        13:44:30   BIC                          -3635.161
Sample:                             0   HQIC                         -3644.288
                               - 1000                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept      0.0354      0.003     12.908      0.000       0.030       0.041
ar.L1          0.9810      0.003    338.304      0.000       0.975       0.987
sigma2         0.0015   9.04e-06    166.245      0.000       0.001       0.002
===================================================================================
Ljung-Box (L1) (Q):                  25.71   Jarque-Bera (JB):           7100021.57
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               0.03   Skew:                            16.69
Prob(H) (two-sided):                  0.00   Kurtosis:                       414.44
===================================================================================

So it has estimated $\phi$ to be close to 1, as expected.

Something that's interesting to notice is that even though the estimated intercept isn't what you were expecting, the model has still done a pretty good job estimating the mean of the series, since $\mu = c / (1 - \phi) = 0.0354 / (1 - 0.9810) = 1.86$. The model is trying to do the best it can, given that the data appears to be approximately non-stationary.

How to move forward with ARIMA

In your case, the answer is just to initialize the simulated data closer to mean, just as @Chris Haug suggested.

More generally, if your dataset was not simulated, you would likely think of the first few observations as outliers. Then, depending on the context, you would either want to model them (for example by using a fat-tailed distribution rather than a Gaussian distribution) or, if they weren't central to the analysis, you might drop them.

Here, if we use your code but drop the first 6 observations, we get results that are in line with what you were expecting:

model_subset = sm.tsa.SARIMAX(S[6:], order=(1,0,0), trend='c')
results_subset = model_subset.fit()
print(results_subset.summary())

yields:

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                      0   No. Observations:                  994
Model:               SARIMAX(1, 0, 0)   Log Likelihood                3193.052
Date:                Mon, 11 Jan 2021   AIC                          -6380.103
Time:                        14:31:47   BIC                          -6365.398
Sample:                             0   HQIC                         -6374.512
                                - 994                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept      1.1044      0.055     19.952      0.000       0.996       1.213
ar.L1          0.4479      0.028     16.180      0.000       0.394       0.502
sigma2       9.48e-05   4.23e-06     22.417      0.000    8.65e-05       0.000
===================================================================================
Ljung-Box (L1) (Q):                   0.02   Jarque-Bera (JB):                 0.25
Prob(Q):                              0.89   Prob(JB):                         0.88
Heteroskedasticity (H):               1.04   Skew:                            -0.04
Prob(H) (two-sided):                  0.72   Kurtosis:                         3.02
===================================================================================

What about regression with ARIMA errors?

There is very little practical difference between using regression with ARIMA errors form (i.e. the constant term as a mean value) versus an ARIMA model with the constant as an intercept. This post by Rob Hyndman notes that "there is not much to choose between the models in terms of forecasting ability, but the additional ease of interpretation in [regression with ARMA errors] makes it attractive.".

For example, we can run regression with AR errors on the stationary subset of the data as follows:

model_subset_reg = sm.tsa.SARIMAX(S[6:], order=(1,0,0), exog=np.ones_like(S[6:]))
results_subset_reg = model_subset_reg.fit()
print(results_subset_reg.summary())

which yields:

                               SARIMAX Results                                
==============================================================================
Dep. Variable:                      0   No. Observations:                  994
Model:               SARIMAX(1, 0, 0)   Log Likelihood                3193.054
Date:                Mon, 11 Jan 2021   AIC                          -6380.109
Time:                        14:35:47   BIC                          -6365.403
Sample:                             0   HQIC                         -6374.518
                                - 994                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          2.0002      0.001   3576.525      0.000       1.999       2.001
ar.L1          0.4500      0.028     16.258      0.000       0.396       0.504
sigma2       9.48e-05   4.23e-06     22.417      0.000    8.65e-05       0.000
===================================================================================
Ljung-Box (L1) (Q):                   0.04   Jarque-Bera (JB):                 0.25
Prob(Q):                              0.83   Prob(JB):                         0.88
Heteroskedasticity (H):               1.04   Skew:                            -0.04
Prob(H) (two-sided):                  0.72   Kurtosis:                         3.02
===================================================================================

This is giving basically the same results as before, and notice that we can recover the intercept directly from the formula $c = (1 - \phi) \mu = (1 - 0.4500) * 2.0002 = 1.10011$, which is essentially identical to the result we got when estimating the intercept directly.