# Confusing SARIMAX Parameter Estimates from Simulated ARMA Data

I've been working on simulating ARMA (Autoregressive Moving Average) time series data in Python and fitting ARIMA models using the SARIMAX class from the statsmodels package. However, I've encountered some puzzling results in the parameter estimates, more specifically when including an MA component in my data generating process I get $$\hat{c} = 0.002, \hat{\theta_{1}} = 0.261, \hat{\theta_2} = -0.069$$ while in fact $$c = 0.002, \theta_{1} = 0.3, \theta_2 = -0.1$$ (even though I have 80,000 observations).

I think it's fair to assume SARIMAX is implemented correctly but at the same time I don't see the error in my simulation. I was hoping maybe someone can spot my error.

Below is my implementation for generating ARMA data and fitting the models:

import numpy as np
from statsmodels.tsa.statespace.sarimax import SARIMAX

def generate_arma_data(intercept, ar_coefficients, ma_coefficients, std_eps, n_periods):

# For later use: ARMA(p, q)
p = len(ar_coefficients)
q = len(ma_coefficients)

# Generate normal errors centered around 0 with the specified standard error
eps = np.random.normal(loc=0, scale=std_eps, size=n_periods+max(p, q))

# Initialize
simulated_data = [0] * max(p, q)

# Iteratively generate the series
for t in range(max(p, q), len(eps)):
# For each AR coefficient we multiply it with the lagged data and sum to get the AR contribution
ar_component = sum(ar_coefficients[j] * simulated_data[t-j-1] for j in range(p))
# For each MA coefficient we multiply it with the lagged error and sum to get the MA contribution
ma_component = sum(ma_coefficients[j] * eps[t-j-1] for j in range(q))
# Finally add the intercept or constant and add the error of the current timestep
datapoint = intercept + ar_component + ma_component + eps[t]
simulated_data.append(datapoint)

return simulated_data[-n_periods::]

# Initialize
intercept = 0.002 # Constant in the ARMA specification
std_eps = 0.001 # Standard deviation of the errors
n_periods = 80000 # Number of observations in our sample
ma_coefs = [0.3, -0.1] # Coefficients for the MA terms

# Generate the data
data_ma_2 = generate_arma_data(intercept, [], ma_coefs, std_eps, n_periods)

# Fit the model
model = SARIMAX(data_ma_2, order=(0, 0, 2), seasonal_order=(0, 0, 0, 0), trend='c')
results_ma = model.fit()

# Display results
print(results_ma.summary()) # Yields terrible parameter estimates (c, ma.L1, ma.L2) = (0.002, 0.261, -0.069)

• Thanks for the comments I have cleaned up the code and made my confusion explicit. Commented Mar 19 at 17:59

If anyone runs into a similar issue: it seems to be the case that the variance of the errors is too low $$(\sigma_{\epsilon}^2 = (0.001)^2)$$ to properly estimate coefficients for the MA terms. Could be a number precision thing I'm not sure about that. Either way, increasing the variance of the errors fixes my issue.