Why Python does not seem to correctly report fitted values in statmodels ARIMA when differencing is involved?

Question

When fitting an ARIMA model using the statsmodels python implementation I see the following behaviour, python does not seem to correctly provide the values for the differenced lags. I am comparing the results with the ones obtained using the R ARIMA implementation. Scenario:

• Seasonal non-stationary data which requires
• Single differencing
• Seasonal first differencing of period 4
• We end up with a ARIMA(0,1,1)(0,1,1,4) model

Please note the issue occurs when using the python ARIMA implementation when importing from statsmodels.tsa.arima.model import ARIMA I've seen webtutorials where it seems to function correctly but they seem to be using the previous implementation from statsmodels.tsa.arima_model import ARIMA, this implementation seems to currently be deprecated https://github.com/statsmodels/statsmodels/issues/3884

I am using statsmodels.version = 0.13.5

Python model code below:

import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
df=pd.read_pickle('ffp2\datasets\euretail.pk')
df.plot()
model=ARIMA(df, order=(0,1,1), seasonal_order=(0,1,1,4)).fit()

However the lags where the differencing occurs they present incorrect fitted values, see below:

pd.concat((df,model.fittedvalues, model.resid), axis=1).head(15)

            value      py_fit   py_resid      r_fit   r_resid
1996-03-31  89.13    0.000000  89.130000  89.078541  0.051459 <---- (Issue Here)
1996-06-30  89.52   89.130003   0.389997  89.496685  0.023315
1996-09-30  89.88   89.520000   0.360000  89.864432  0.015568
1996-12-31  90.12   89.879998   0.240002  90.143352 -0.023352
1997-03-31  89.19  134.685000 -45.495000  89.380354 -0.190354 <---- (Issue Here)
1997-06-30  89.78   89.579993   0.200007  89.621984  0.158016
1997-09-30  90.03   90.193547  -0.163547  90.164354 -0.134354
1997-12-31  90.38   90.222837   0.157163  90.251028  0.128972
1998-03-31  90.27   89.463010   0.806990  89.602346  0.667654
1998-06-30  90.77   90.951737  -0.181737  90.937687 -0.167687
1998-09-30  91.85   91.019609   0.830391  91.077857  0.772143
1998-12-31  92.51   92.389250   0.120750  92.397782  0.112218
1999-03-31  92.21   92.044928   0.165072  92.056141  0.153859
1999-06-30  92.52   92.752252  -0.232252  92.744481 -0.224481
1999-09-30  93.62   93.066856   0.553144  93.083940  0.536060

Python fitted model below (for reference)

Fitted model results below:

                                    SARIMAX Results                                    
=======================================================================================
Dep. Variable:                           value   No. Observations:                   64
Model:             ARIMA(0, 1, 1)x(0, 1, 1, 4)   Log Likelihood                 -34.642
Date:                         Thu, 01 Dec 2022   AIC                             75.285
Time:                                 14:09:56   BIC                             81.517
Sample:                             03-31-1996   HQIC                            77.718
                                  - 12-31-2011                                         
Covariance Type:                           opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ma.L1          0.2903      0.155      1.872      0.061      -0.014       0.594
ma.S.L4       -0.6912      0.132     -5.250      0.000      -0.949      -0.433
sigma2         0.1810      0.034      5.316      0.000       0.114       0.248
===================================================================================
Ljung-Box (L1) (Q):                   0.25   Jarque-Bera (JB):                 1.91
Prob(Q):                              0.62   Prob(JB):                         0.38
Heteroskedasticity (H):               0.76   Skew:                            -0.22
Prob(H) (two-sided):                  0.54   Kurtosis:                         3.77
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

R model code for reference below

fit <- Arima(euretail, order=c(0,1,1), seasonal=c(0,1,1))
fitted(fit)
resid(fit)

Dataset below:

      Qtr1   Qtr2   Qtr3   Qtr4
1996  89.13  89.52  89.88  90.12
1997  89.19  89.78  90.03  90.38
1998  90.27  90.77  91.85  92.51
1999  92.21  92.52  93.62  94.15
2000  94.69  95.34  96.04  96.30
2001  94.83  95.14  95.86  95.83
2002  95.73  96.36  96.89  97.01
2003  96.66  97.76  97.83  97.76
2004  98.17  98.55  99.31  99.44
2005  99.43  99.84 100.32 100.40
2006  99.88 100.19 100.75 101.01
2007 100.84 101.34 101.94 102.10
2008 101.56 101.48 101.13 100.34
2009  98.93  98.31  97.67  97.44
2010  96.53  96.56  96.51  96.70
2011  95.88  95.84  95.79  95.94

Answer 1

There are two issues to consider here:

First, due to the differencing operator, the residuals for the first d + sD observations should be ignored because the predictions remain associated with a diffuse prior, which means that although the predicted mean is defined (and this is what is reported), the associated variance is infinite. Here d=1, D=1, and s=4, so that would be the first 5 observations. To avoid confusion, Statsmodels should probably set those to be NaN, since the values are not generally useful.

Second, I don't know exactly what R is returning for predictions and residuals here, but the difference between the Statsmodels results is not due to a bug in Statsmodels. For example, you can use Stata or Eviews to compute the predictions and residuals and they will match Statsmodels.

It is easy to see that R is doing something non-standard here by looking at its predictions. Because the predictions are one-step-ahead forecasts, the prediction of the very first observation cannot be be based on any data, and the model doesn't include a constant term. This is why the prediction from Statsmodels is equal to 0. So the question is, what is R's prediction of 89.078541 based on?

I guess that R is doing something additional to try and return some kind of predictions / residuals that are more useful than the diffuse predictions, but this is not documented anywhere that I can see.