I have data concerning, say, ice cream sales and wish to predict future sales and also, to quantify the relationship between temperature and sales. The data has daily seasonality.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pmdarima as pm
import statsmodels.api as sm


def generate_synthetic_data(n_days=10, sd=5):
    day_time, sales, temperature = [], [], []
    for i in range(n_days): 
        day_time.extend([str(i)+'_0000', str(i)+'_0400', str(i)+'_0800', str(i)+'_1200', str(i)+'_1600', str(i)+'_2000'])
        #time.extend(['0000', '0400', '0800', '1200', '1600', '2000'])
        sales.extend([260, 322, 530, 547, 379, 234] + np.random.normal(0,sd*10,6))
        temperature.extend([15, 18, 25, 22, 21, 12] + np.random.normal(0,sd,6))
    df = pd.DataFrame({'day_time': day_time, 'sales': sales, 'temperature': temperature})
    return df

df = generate_synthetic_data(n_days)

A quick plot of the first 5 days shows what the data looks like: enter image description here

Clearly, sales and temperature are related. I fit a simple linear regression model.

slope, intercept, r_value, p_value, std_err = stats.linregress(df['temperature'], df['sales'])
# Slope is 11.8.

So I conclude that a unit increase in temperature results in an increase of 11.8 in ice cream sales. So far so good.

Now, I proceed to forecasting future ice cream sales. I split the data into train-test sets and fit a SARIMAX model, with temperature as an exogenous variable.

train_df = df.iloc[:294]
test_df = df.iloc[294:]
model = pm.auto_arima(y=train_df['sales'], X=train_df['temperature'].to_frame(), seasonal=True, m=6)

The output of this model is as follows:

enter image description here

The predictions on the test set look good, but this is not what I'm concerned about. The coefficient for temperature has now drastically changed (from 11.8 to 0.85). It wouldn't make much sense to change my interpretation now to a unit increase in temperature resulting in an additional 0.85 ice cream sales. Is this huge difference down to:

  1. Collinearity between sales and temperature?
  2. Misinterpreting ARIMAX coefficients, and I'd be better off with a linear regression with AR(I)MA errors? If this is the case, is there a package that implements this in Python? Or should I turn to R?
  3. Both of the above?

1 Answer 1


Reason 1 makes sense to me. Your first model has the following format:

$y_t = \beta_0 + \theta * x_t$

And if I am interpreting correctly your SARIMAX model format is:

$y_t = \beta_0 + sazon_t + a_1*y_{t-1} + \theta * x_t$

So having same $\theta$ values on both models your require the term $sazon_t + a_1*y_{t-1}$ to be estimated as 0 (zero) and that doesn't look like the case, as your series clearly has sazonality.

You could validate the sarimax coeficients understanding estimating a model without sazonality and check if coeficient to temperature approaches the previous model one.

Other things to be aware of:

  • pvalues indicate if the coefficient impact is being confidently estimated as different than zero by the model. The SARIMAX model seems to be estimating intercept, temperature and AR[1] as zero.
  • Using the temperature model as exogenous variable on the model requires that you have this variable value on prediction time. Just be sure this is the case to use this.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.