What is the exact defenition of ARMA model in statsmodels and how to predict the next value using the summary of the fitted model manually?

Question

I've fitted a time series (Y) on the ARMA(2,1) model using statsmodels in python. let's leave alone that the selected order is not the best for my time series (Y).

The summary of the fitted model can be seen below. The AR and MA coefficients are shown as ar.L1.Y, ar.L2.Y, and ma.L1.Y. The constant value is indicated as const. The standard deviation of the white noise is also given by S.D. of innovations.

Therefore, we may say that our ARMA model is defined as follows:

$Y_{t}= C+ A_{1}Y_{t-1}+A_{2}Y_{t-2}+B_{1}e_{t-1}+e_{t}$

However, I am not certain about it, since I could not find the ARMA model definition in statsmodels package. I know there is no unique definition for the ARMA model through different software packages. For instance, the MA coefficients are defined by a negative sign in some packages and textbooks.

Questions:

1- Is my interpretation of ARMA coefficients and parameters with respect to the summary correct? (Paragraph 2)

2- What is statsmodels definition for ARMA's equation?

3- Lets say we have fitted the model using $(Y_{0},..., Y_{10})$. How can we predict $Y_{11}$?

Well, my guess is to do as follows. But not sure, especially about the value of $e_{10}$.

$Y_{11}= C+ A_{1}Y_{10}+A_{2}Y_{9}+B_{1}e_{10}+e_{11}$, in which $e_{11}$ is a generated random number based on $N(0,\sigma)$. However, what is $e_{10}$ and how can I calculate it?

                                   ARMA Model Results                                  
=======================================================================================
Dep. Variable:                           Y       No. Observations:                  359
Model:                              ARMA(2, 1)   Log Likelihood                -127.666
Method:                                css-mle   S.D. of innovations              0.344
Date:                         Thu, 12 Mar 2020   AIC                            265.331
Time:                                 09:35:15   BIC                            284.748
Sample:                             01-01-2017   HQIC                           273.052
                                  - 01-04-2017                                         
================================================================================================
                              coef       std err      z        P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------------------
const                        4.4386      0.911      4.873      0.000       2.653       6.224
ar.L1.Y                      1.0940      1.398      0.783      0.434      -1.645       3.833
ar.L2.Y                     -0.1096      1.373     -0.080      0.936      -2.801       2.582
ma.L1.Y                     -0.1028      1.398     -0.074      0.941      -2.843       2.637
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1            1.0179           +0.0000j            1.0179            0.0000
AR.2            8.9667           +0.0000j            8.9667            0.0000
MA.1            9.7242           +0.0000j            9.7242            0.0000
-----------------------------------------------------------------------------

Answer 1

1. Your interpretation is correct
2. It is essentially the equation that you have defined above but you must move the constant. It is more like $Y_t = Y'_t+C$ where $Y'_t = A_1Y'_{t−1}+A_2Y'_{t−2}+B_1e_{t−1}+e_t$
3. In your prediction of $Y_{11}$ you would set $e_{11}$ to its mean value of 0 rather than drawing a random number for it from the Gaussian distribution. $e_{10}$ can only be estimated and is generally estimated by iteration when fitting the model itself. You can read more about actually fitting the model here. The estimates themselves are stored in the fitted statsmodels object .resid property.

Here is a minimal example to play with

import matplotlib.pyplot as plt
import numpy as np
from statsmodels.tsa.arima_process import arma_generate_sample
from statsmodels.tsa.arima_model import ARMA

np.random.seed(878015)

# create simulated ARMA(1,1) dataset
ar_coefs = [1, -0.5] # ar1 = 0.5 -> function used makes us pass in negative 
ma_coefs = [1, 0.2] # ma1 = 0.2
sigma = 0.2 # scale of e_t
nsample = 1000

x = np.arange(nsample)
y = arma_generate_sample(ar_coefs, ma_coefs, nsample=nsample, sigma=sigma, )

plt.figure()
plt.plot(x,y)
plt.ylabel(r'$y_t$')
plt.xlabel(r'$t$')
plt.show()

# fit ARMA model to generated data
model = ARMA(y, order=(1,1))
results = model.fit()
print(results.summary())

#                              ARMA Model Results                              
#==============================================================================
#Dep. Variable:                      y   No. Observations:                 1000
#Model:                     ARMA(1, 1)   Log Likelihood                 201.309
#Method:                       css-mle   S.D. of innovations              0.198
#Date:                Fri, 20 Mar 2020   AIC                           -394.617
#Time:                        12:14:42   BIC                           -374.986
#Sample:                             0   HQIC                          -387.156
#                                                                              
#==============================================================================
#                 coef    std err          z      P>|z|      [0.025      0.975]
#------------------------------------------------------------------------------
#const          0.0088      0.014      0.631      0.528      -0.018       0.036
#ar.L1.y        0.4563      0.046      9.989      0.000       0.367       0.546
#ma.L1.y        0.2102      0.050      4.219      0.000       0.113       0.308
#                                    Roots                                    
#=============================================================================
#                  Real          Imaginary           Modulus         Frequency
#-----------------------------------------------------------------------------
#AR.1            2.1915           +0.0000j            2.1915            0.0000
#MA.1           -4.7565           +0.0000j            4.7565            0.5000
#-----------------------------------------------------------------------------
#5.551115123125783e-17

# collect fitted parameters
fit_c = results.params[0] # = 0.0088
fit_ar1 = results.arparams[0] # = 0.4563
fit_ma1 = results.maparams[0] # = 0.2102

# make predictions for last 10 time steps
n=10
y_hat_manual = (y[-n-1:-1]-fit_c)*fit_ar1 + results.resid[-n-1:-1]*fit_ma1 + fit_c

# make same predictions using statsmodels and compare 
y_hat_sm = results.predict(start=-n, dynamic=False)

pred_resid = y_hat_manual - y_hat_sm
print(abs(pred_resid).max())
# 5.551115123125783e-17