Workshop 4, Stats for AI - Bloque 2

Stationarity, cointegration and ARIMA/SARIMA forecasting models

Author

Alberto Dorantes D., Ph.D.

Published

October 1, 2026

Abstract

In this workshop (4 sessions) we learn what makes a time series stationary, how to transform and test a series for stationarity (Dickey-Fuller and KPSS tests), why regressions between non-stationary series can be spurious, and how to test for cointegration. Then we learn the Box-Jenkins methodology to build ARIMA/SARIMA models: we recognize AR and MA patterns in simulated data, model the IGAE of Mexico, evaluate the forecasts out of sample against simple benchmarks and a machine-learning model, and add an external explanatory variable (ARIMAX).

1 Directions for this Workshop

This workshop covers four class sessions (Parts A to D). Together with Workshop 3, it is the evidence for the topic of time series. The grading criteria are the same as in Workshop 3: code replication (20%), your personal notes and brief explanations of each step (30%), and the solution of the challenges (50%).

You have to submit your .ipynb file through Canvas.

2 Learning objectives

Part A (session 1): Stationarity

Define a stationary series and explain why the AR(1) model is stationary only when |\phi_1|<1.
Compute the long-run mean and the half-life of an AR(1) process.
Transform a business series (logs, first difference, seasonal difference) and interpret the result as a % change.
Run and interpret the Augmented Dickey-Fuller (ADF) and KPSS tests, choosing the right specification.

Part B (session 2): Spurious regression and cointegration

Show by simulation why regressions between non-stationary series can be spurious.
Test for cointegration with the Engle-Granger method, and interpret an error-correction model.
Connect these ideas with common pitfalls of machine-learning models trained on time series.

Part C (session 3): ARIMA models, the building blocks

Write and interpret AR(p), MA(q) and ARMA(p,q) models, and the ARIMA/SARIMA notation.
Read the ACF and PACF plots and recognize AR, MA and ARMA patterns.
Find a stationary version of a real series and choose the differencing parameters d and D.

Part D (session 4): Building and evaluating a forecasting model

Calibrate a SARIMA model using residual diagnostics (ACF/PACF and Ljung-Box) and information criteria (AIC).
Interpret the coefficients of the final model.
Evaluate forecasts out of sample against a naive benchmark and a machine-learning model.
Add an external variable (ARIMAX) and understand what it requires for forecasting.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import requests
from io import BytesIO, StringIO
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.graphics.tsaplots import plot_acf

headers = {'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/120.0 Safari/537.36'}

3 PART A. Stationarity

3.1 From the random walk to the AR(1) model

In Workshop 3 we modeled the log of the S&P 500 as a random walk with drift:

Y_t = \phi_0 + Y_{t-1} + \varepsilon_t

We saw that the uncertainty about Y_t keeps growing with time (SD=\sigma_\varepsilon\sqrt{t}): the series has no level it returns to. Now let’s give a free coefficient to Y_{t-1}:

Y_t = \phi_0 + \phi_1 Y_{t-1} + \varepsilon_t

This is the autoregressive model of order 1, or AR(1). The random walk is the special case \phi_1=1.

3.2 What is a stationary series?

A series is (weakly) stationary if its statistical behavior does not depend on the period we look at. Formally, three conditions must hold:

The mean is the same for every period: E[Y_t]=\mu for all t.
The variance is the same for every period: Var(Y_t)=\sigma_Y^2 for all t.
The covariance between Y_t and Y_{t-k} depends only on the distance k (the lag), not on t.

Intuitively, if you cut a stationary series into pieces, every piece looks statistically similar. A series with a trend, or whose variance grows over time like the random walk, is non-stationary.

3.3 The AR(1) is stationary only if |\phi_1|<1

If the AR(1) is stationary, its mean is the same in t and t-1. Taking expected values of the equation:

\mu = \phi_0 + \phi_1\mu \quad\Rightarrow\quad \mu=\frac{\phi_0}{1-\phi_1}

This long-run mean exists only if \phi_1\neq1. With \phi_1=1 (random walk) there is no mean to return to. Similarly, the variance of a stationary AR(1) is

Var(Y_t)=\frac{\sigma_\varepsilon^2}{1-\phi_1^2}

which is positive and finite only if |\phi_1|<1. Then:

|\phi_1|<1: stationary. The series is mean-reverting: after a shock it returns to \mu.
\phi_1=1: unit root (random walk). Non-stationary; shocks are permanent.
|\phi_1|>1: explosive. Rarely useful for business data.

How fast does the series return to its mean? A shock today still has a fraction \phi_1^k of its effect k periods later. The half-life is the number of periods until half of the effect of a shock has disappeared:

\phi_1^{k}=0.5 \quad\Rightarrow\quad k_{1/2}=\frac{\ln(0.5)}{\ln(|\phi_1|)}

3.4 In-class exercise: simulating AR(1) processes

The following function simulates an AR(1). We pass the same shocks to every simulation so that the only difference between them is \phi_1:

def simulate_ar1(phi0, phi1, shocks, y0=0.0):
    """Simulate Y_t = phi0 + phi1*Y_{t-1} + e_t using the given shocks."""
    y = np.zeros(len(shocks))
    y[0] = y0
    for t in range(1, len(shocks)):
        y[t] = phi0 + phi1 * y[t-1] + shocks[t]
    return y

n = 600
rng = np.random.default_rng(10)
shocks = rng.normal(0, 0.01, n)    # sigma of the shock = 0.01
phi0 = 0.01
phi1_values = [0.10, 0.90, 0.98, 1.00]

fig, axes = plt.subplots(4, 1, figsize=(12, 10), sharex=True)
sims = {}
for ax, p1 in zip(axes, phi1_values):
    sims[p1] = simulate_ar1(phi0, p1, shocks)
    ax.plot(sims[p1])
    if p1 < 1:
        ax.axhline(phi0 / (1 - p1), color='red', linestyle='--', label=r'long-run mean $\phi_0/(1-\phi_1)$')
        ax.legend(loc='lower right')
    ax.set_title(f'phi1 = {p1}')
    ax.grid(True)
plt.tight_layout()
plt.show()

Let’s compare the theory with the simulated series. For the sample statistics we use only the second half of each series, so that the starting value Y_0=0 does not distort them:

rows = []
for p1 in phi1_values:
    y = sims[p1][n // 2:]
    rows.append({
        'phi1': p1,
        'theoretical mean': phi0 / (1 - p1) if p1 < 1 else np.nan,
        'sample mean (2nd half)': y.mean(),
        'lag-1 autocorrelation': pd.Series(y).autocorr(1),
        'half-life (periods)': np.log(0.5) / np.log(p1) if p1 < 1 else np.inf,
    })
ar1_table = pd.DataFrame(rows).round(4)
ar1_table

   phi1  theoretical mean  ...  lag-1 autocorrelation  half-life (periods)
0  0.10            0.0111  ...                 0.0202               0.3010
1  0.90            0.1000  ...                 0.8303               6.5788
2  0.98            0.5000  ...                 0.9217              34.3096
3  1.00               NaN  ...                 0.9999                  inf

[4 rows x 5 columns]

Answer in class:

What happens to the series when \phi_1 is small (0.10)? And when it is close to 1?
For \phi_1=0.98, how many periods does it take for half of a shock to disappear? What does that imply for a forecast 12 periods ahead?
The series with \phi_1=0.98 and the one with \phi_1=1 look similar over short windows. Why will it be hard for a statistical test to tell them apart?

3.5 Autocorrelation: a first look

The autocorrelation at lag k is the correlation between Y_t and Y_{t-k}. For a stationary AR(1), the theory says \rho_k=\phi_1^k: it decays geometrically. For a random walk, the sample autocorrelations decay very slowly, which is a visual sign of a unit root. The autocorrelation function (ACF) plots \rho_k for many lags:

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(sims[0.90], lags=40, zero=False, ax=axes[0], title='ACF, phi1 = 0.90 (stationary)')
plot_acf(sims[1.00], lags=40, zero=False, ax=axes[1], title='ACF, phi1 = 1 (random walk)')
plt.tight_layout()
plt.show()

The blue shaded area is the 95% confidence band: bars outside it are significantly different from zero. We will use the ACF intensively in Part C.

3.6 Real data: sales of consumer products

We work again with the monthly dataset of volume sales (in tons) of the brands of one consumer-product category:

url = "http://www.apradie.com/ec2004/salesbrands3.xlsx"
response = requests.get(url, headers=headers)
response.raise_for_status()
sales_brand = pd.read_excel(BytesIO(response.content), sheet_name='Resultado')

sales_brand['date'] = pd.to_datetime(sales_brand['date'])
sales_brand = sales_brand.set_index('date').sort_index().asfreq('MS')
brand_cols = [c for c in sales_brand.columns if c.startswith('qbrand')]
print(brand_cols)

['qbrand1tons', 'qbrand2tons', 'qbrand3tons', 'qbrand4tons', 'qbrand5tons', 'qbrand6tons']

Let’s plot brands 2 and 6:

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
for ax, col in zip(axes, ['qbrand2tons', 'qbrand6tons']):
    ax.plot(sales_brand.index, sales_brand[col])
    ax.set_title(f'Volume sales: {col}')
    ax.grid(True)
plt.tight_layout()
plt.show()

Both series show trend and seasonality, so they are not stationary: their mean changes over time. Most business variables (sales, revenue, costs, prices) behave like this.

3.7 Transformations that make a series stationary

3.7.1 Logs and differences

The usual recipe for business series is:

Take the natural log. This turns proportional changes into additive ones and stabilizes the variance when it grows with the level.
Take a difference of the log.

Why is the difference of logs a % change? If R_t is the simple % change, then Y_t=Y_{t-1}(1+R_t) and

\ln(Y_t)-\ln(Y_{t-1})=\ln\left(\frac{Y_t}{Y_{t-1}}\right)=\ln(1+R_t)\approx R_t \quad \text{(for small } R_t\text{)}

The log difference is the continuously compounded % change. For small changes it is almost equal to the simple % change. For large changes it is not, and we convert back with R_t=e^{\Delta \ln Y_t}-1. For example, a log difference of -0.25 is a simple change of e^{-0.25}-1=-22.1\%, not -25\%.

With monthly data there are two useful differences:

First difference of the log: \ln Y_t-\ln Y_{t-1} = monthly % change.
Seasonal difference of the log: \ln Y_t-\ln Y_{t-12} = annual % change, month by month. Since it compares each month with the same month of the previous year, it also removes the seasonality.

ln_sales = np.log(sales_brand[brand_cols])
msales_growth = ln_sales.diff(1)    # monthly % change
ysales_growth = ln_sales.diff(12)   # annual % change

fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)
for col in ['qbrand2tons', 'qbrand6tons']:
    axes[0].plot(msales_growth.index, msales_growth[col], label=col)
    axes[1].plot(ysales_growth.index, ysales_growth[col], label=col)
axes[0].set_title('Monthly % change (first difference of log)')
axes[1].set_title('Annual % change (seasonal difference of log)')
for ax in axes:
    ax.axhline(0, color='black', linewidth=0.8)
    ax.legend()
    ax.grid(True)
plt.tight_layout()
plt.show()

The monthly % change no longer has a trend, but it shows large, regular peaks and valleys: the seasonality is still there. The annual % change has no trend and no seasonal pattern. Whether these series are stationary is hard to judge by eye. We need a statistical test.

3.8 The Augmented Dickey-Fuller (ADF) test

The Dickey-Fuller test (Dickey and Fuller 1979) starts from the AR(1) model and subtracts Y_{t-1} from both sides:

\Delta Y_t = \phi_0 + (\phi_1-1)Y_{t-1} + \varepsilon_t = \phi_0 + \gamma Y_{t-1} + \varepsilon_t

H_0: \gamma=0 (that is, \phi_1=1): the series has a unit root, so it is non-stationary.
H_A: \gamma<0 (that is, \phi_1<1): the series is stationary.

That is why it is called a unit-root test. The Augmented version adds lags of \Delta Y_t to the regression to absorb short-term autocorrelation. adfuller chooses how many lags to add with the AIC criterion.

Under H_0 the test statistic does not follow a t distribution, so we use the special Dickey-Fuller critical values (adfuller computes the p-value for us).

3.8.1 Which version of the test?

The regression argument decides which deterministic terms enter the test regression. The choice matters:

`regression`	Test regression includes	Use it when the series…
`'n'`	nothing	fluctuates around zero (rare)
`'c'`	a constant	fluctuates around a non-zero mean, no trend. Typical for growth rates / returns
`'ct'`	a constant and a linear trend	has a visible trend. Typical for log levels (then H_A = stationary around a trend)

A helper function to report the test:

from statsmodels.tsa.stattools import adfuller, kpss

def adf_test(x, regression='c'):
    """Run the ADF test and return the key results in a dictionary."""
    res = adfuller(x.dropna(), regression=regression, autolag='AIC')
    return {'ADF stat': res[0], 'p-value': res[1], 'lags used': res[2], 'obs': res[3]}

Now we apply it to three versions of brand 2: the log level, the monthly % change and the annual % change:

q2 = 'qbrand2tons'
adf_q2 = pd.DataFrame({
    'log level (ct)': adf_test(ln_sales[q2], 'ct'),
    'monthly % change (c)': adf_test(msales_growth[q2], 'c'),
    'annual % change (c)': adf_test(ysales_growth[q2], 'c'),
}).T.round(4)
adf_q2

                      ADF stat  p-value  lags used   obs
log level (ct)         -0.3013   0.9895       11.0  42.0
monthly % change (c)   -2.7994   0.0584       11.0  41.0
annual % change (c)    -2.4884   0.1183        3.0  38.0


p_level = adf_q2.loc['log level (ct)', 'p-value']
p_annual = adf_q2.loc['annual % change (c)', 'p-value']
verdict_level = "reject" if p_level < 0.05 else "cannot reject"
verdict_annual = "reject" if p_annual < 0.05 else "cannot reject"

How to read the table:

If the p-value is less than 0.05, we reject H_0 (unit root) at the 5% level and treat the series as stationary.
If the p-value is greater than 0.05, we cannot reject the unit root and treat the series as non-stationary.

For brand 2, with the log level the p-value is 0.9895, so we cannot reject the unit root. With the annual % change the p-value is 0.1183, so we cannot reject the unit root.

Two cautions about the ADF test

Low power near the unit root. Remember the in-class exercise: a series with \phi_1=0.98 looks almost like a random walk. With few observations the ADF test often cannot reject H_0 even if the series is stationary. “Cannot reject” is not “proven non-stationary”.
Seasonality. A strong seasonal pattern in the monthly % change can distort the test. That is one more reason to prefer the annual % change for monthly business data.

3.8.2 A cross-check: the KPSS test

The KPSS test (Kwiatkowski et al. 1992) reverses the hypotheses: H_0 = the series is stationary. Using both tests gives a more reliable conclusion:

ADF (H_0: unit root)	KPSS (H_0: stationary)	Conclusion
reject	do not reject	stationary
do not reject	reject	non-stationary
reject	reject	ambiguous: possibly a structural break or long memory
do not reject	do not reject	ambiguous: not enough data to decide

kpss_stat, kpss_p, _, _ = kpss(ysales_growth[q2].dropna(), regression='c', nlags='auto')
print(f"KPSS for the annual % change of brand 2: stat = {kpss_stat:.3f}, p-value = {kpss_p:.3f}")

KPSS for the annual % change of brand 2: stat = 0.199, p-value = 0.100

Note: kpss reports p-values only between 0.01 and 0.10 (it prints a warning when the true value is outside that range).

4 CHALLENGE 1 (Part A)

Work in your Notebook. Report your numbers and interpret them in your own words.

Build one table with the ADF p-values of all brands in the dataset (brand_cols) for three versions of each series: log level ('ct'), monthly % change ('c') and annual % change ('c'). Tip: loop over brand_cols and use the adf_test function.
Add a column with the KPSS p-value of the annual % change.
For each brand, decide which version you would use for modeling and justify it in one sentence. Do all brands lead to the same decision? Identify at least one brand where the tests are ambiguous and explain what could cause it (look at its plot).
Using the in-class AR(1) simulation, compute the half-life of a shock for \phi_1=0.5, 0.9 and 0.99 with monthly data. Express it in months and explain what it would mean for the sales of a brand.

5 PART B. Spurious regression and cointegration

5.1 Spurious regression: a simulation

Let’s simulate two random walks that are completely independent: they are built with different random shocks, so there is no relationship between them. Then we regress one on the other:

rng = np.random.default_rng(7)
n = 200
x_rw = np.cumsum(rng.normal(0, 1, n))
y_rw = np.cumsum(rng.normal(0, 1, n))
sim_df = pd.DataFrame({'y': y_rw, 'x': x_rw})

plt.figure(figsize=(12, 4))

<Figure size 1200x400 with 0 Axes>

plt.plot(sim_df['y'], label='y (random walk 1)')

[<matplotlib.lines.Line2D object at 0x000001D6AB23D010>]

plt.plot(sim_df['x'], label='x (random walk 2)')

[<matplotlib.lines.Line2D object at 0x000001D6AC36CEC0>]

plt.title('Two independent random walks')

Text(0.5, 1.0, 'Two independent random walks')

plt.legend()

<matplotlib.legend.Legend object at 0x000001D6AB23D2B0>

plt.grid(True)
plt.show()

spur = smf.ols('y ~ x', data=sim_df).fit()
print(spur.summary().tables[1])

==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -5.8696      0.686     -8.552      0.000      -7.223      -4.516
x              0.2813      0.034      8.285      0.000       0.214       0.348
==============================================================================

print(f"R-squared = {spur.rsquared:.3f}")

R-squared = 0.257

Very often (try other seeds) the regression finds a “significant” relationship, with a high R-squared, between two series that have nothing to do with each other. This is a spurious regression (Granger and Newbold 1974).

How often does this happen? We repeat the experiment 500 times for different sample sizes, and count how many times the p-value of the slope is below 0.05. If the test worked properly, that should happen about 5% of the time, because the true slope is zero:

def spurious_rate(n, reps=500, differences=False, seed=0):
    """Share of regressions with p<0.05 between two INDEPENDENT random walks."""
    rng = np.random.default_rng(seed)
    count = 0
    for _ in range(reps):
        x = np.cumsum(rng.normal(size=n))
        y = np.cumsum(rng.normal(size=n))
        if differences:
            x, y = np.diff(x), np.diff(y)
        p = sm.OLS(y, sm.add_constant(x)).fit().pvalues[1]
        count += p < 0.05
    return count / reps

spurious_table = pd.DataFrame({
    'n': [50, 200, 1000],
    'p<0.05 in LEVELS': [spurious_rate(k) for k in [50, 200, 1000]],
    'p<0.05 in DIFFERENCES': [spurious_rate(k, differences=True) for k in [50, 200, 1000]],
})
spurious_table

      n  p<0.05 in LEVELS  p<0.05 in DIFFERENCES
0    50             0.700                  0.072
1   200             0.818                  0.050
2  1000             0.922                  0.050


rate_200 = float(spurious_table.loc[spurious_table['n'] == 200, 'p<0.05 in LEVELS'].iloc[0])

With n=200, the regression in levels declares a “significant” relationship in about 82% of the cases instead of 5%. The problem gets worse with more data. In differences (which are stationary) the rate is close to the expected 5%.

Why? If y and x are independent random walks, the regression error y_t-\beta_0-\beta_1x_t is itself a random walk. It is non-stationary and highly autocorrelated, which breaks the OLS assumptions. The t-statistics no longer follow a t distribution and become larger as n grows. A classic warning sign: an R^2 higher than the Durbin-Watson statistic.

from statsmodels.stats.stattools import durbin_watson
print(f"R-squared = {spur.rsquared:.3f},  Durbin-Watson = {durbin_watson(spur.resid):.3f}")

R-squared = 0.257,  Durbin-Watson = 0.061

A Durbin-Watson statistic close to 0 means that the residuals are strongly positively autocorrelated.

5.2 Spurious regression with real data

We download from the World Bank two annual indicators for Mexico that should not be directly related: the infant mortality rate (deaths per 1,000 live births) and an index of the value of merchandise exports. We use the World Bank API directly:

def get_wb(indicator, country='MEX', start=1980, end=2024):
    """Download one World Bank indicator for one country as a pandas Series."""
    url = (f"https://api.worldbank.org/v2/country/{country}/indicator/{indicator}"
           f"?format=json&date={start}:{end}&per_page=500")
    js = requests.get(url, timeout=60).json()
    s = pd.Series({int(d['date']): d['value'] for d in js[1]}, dtype=float)
    return s.sort_index()

wb = pd.DataFrame({
    'infant_mortality': get_wb('SP.DYN.IMRT.IN'),
    'exports': get_wb('TX.VAL.MRCH.XD.WD'),
}).dropna()
wb.tail()

      infant_mortality  exports
2020              13.4    109.6
2021              13.1    130.1
2022              12.8    151.8
2023              12.4    155.8
2024              12.1    162.2

fig, ax1 = plt.subplots(figsize=(12, 4))
ax1.plot(wb.index, wb['infant_mortality'], color='tab:blue', label='Infant mortality (left)')

[<matplotlib.lines.Line2D object at 0x000001D6AB202A50>]

ax2 = ax1.twinx()
ax2.plot(wb.index, wb['exports'], color='tab:orange', label='Merchandise exports index (right)')

[<matplotlib.lines.Line2D object at 0x000001D6AB23CEC0>]

ax1.set_title('Mexico: infant mortality and merchandise exports')

Text(0.5, 1.0, 'Mexico: infant mortality and merchandise exports')

fig.legend(loc='upper center')

<matplotlib.legend.Legend object at 0x000001D6AB23CC20>

ax1.grid(True)
plt.show()

We regress the log of exports on the log of infant mortality, first in levels and then in differences (annual % changes):

wb_log = np.log(wb)
reg_levels = smf.ols('exports ~ infant_mortality', data=wb_log).fit()
reg_diffs = smf.ols('exports ~ infant_mortality', data=wb_log.diff().dropna()).fit()

wb_results = pd.DataFrame({
    'slope': [reg_levels.params['infant_mortality'], reg_diffs.params['infant_mortality']],
    'p-value': [reg_levels.pvalues['infant_mortality'], reg_diffs.pvalues['infant_mortality']],
    'R-squared': [reg_levels.rsquared, reg_diffs.rsquared],
    'Durbin-Watson': [durbin_watson(reg_levels.resid), durbin_watson(reg_diffs.resid)],
}, index=['levels (logs)', 'differences (% changes)']).round(4)
wb_results

                          slope  p-value  R-squared  Durbin-Watson
levels (logs)           -2.0169   0.0000     0.9212         1.5875
differences (% changes)  3.4032   0.6791     0.0103         2.2280

Discuss in class: Did the regression in levels find a “significant” relationship? Does it make sense that lower infant mortality causes higher exports? What happens to the relationship in differences?

Both series move with the long-term development of the country: one goes down and the other goes up over the decades. A regression in levels picks up these common trends, not a direct relationship. But be careful: not every regression between non-stationary series is spurious. The exception is cointegration.

5.3 Cointegration

5.3.1 The idea: a drunk and her dog

Two series are cointegrated if each one is non-stationary (has a unit root), but there is a linear combination of them that is stationary (Engle and Granger 1987):

Y_t \text{ and } X_t \text{ are I(1)}, \quad\text{but}\quad u_t = Y_t-\beta_0-\beta_1X_t \text{ is stationary (I(0))}

A classic image (Murray 1994): a drunk walking home with her dog. Each one wanders like a random walk, so you cannot predict where either will be in one hour. But they never get too far from each other: when they separate, the dog runs back to her (or she calls the dog). The distance between them is stationary. There is a long-run equilibrium, and deviations from it are temporary.

When two series are cointegrated, the regression in levels is not spurious: it estimates the long-run equilibrium relationship.

5.3.2 A simulated cointegrated pair

rng = np.random.default_rng(3)
n = 300
x_c = np.cumsum(rng.normal(0, 1, n))                       # random walk
u_c = simulate_ar1(0, 0.5, rng.normal(0, 1, n))           # stationary deviation
y_c = 1 + 2 * x_c + u_c                                    # cointegrated with x
coint_df = pd.DataFrame({'y': y_c, 'x': x_c})

plt.figure(figsize=(12, 4))

<Figure size 1200x400 with 0 Axes>

plt.plot(coint_df['y'], label='y')

[<matplotlib.lines.Line2D object at 0x000001D6AB0C8C20>]

plt.plot(2 * coint_df['x'], label='2x')

[<matplotlib.lines.Line2D object at 0x000001D6AB0C8D70>]

plt.title('A cointegrated pair: both wander, but y - 2x stays close to a constant')

Text(0.5, 1.0, 'A cointegrated pair: both wander, but y - 2x stays close to a constant')

plt.legend()

<matplotlib.legend.Legend object at 0x000001D6AB0C8980>

plt.grid(True)
plt.show()

5.3.3 The Engle-Granger test

The steps are:

Check that both series are I(1): ADF does not reject in levels, but rejects in first differences.
Run the regression in levels: Y_t=\beta_0+\beta_1X_t+u_t.
Test whether the residuals \hat u_t are stationary with a unit-root test.

There is a subtlety in step 2. The residuals come from a regression that chose \beta to make them as small and as stationary-looking as possible. Because of that, the usual ADF critical values reject too often. We must use the Engle-Granger critical values, which are more demanding. The function coint from statsmodels runs steps 1 and 2 with the correct critical values:

from statsmodels.tsa.stattools import coint

# Wrong way (for comparison): plain ADF on the residuals
resid_c = smf.ols('y ~ x', data=coint_df).fit().resid
p_wrong = adfuller(resid_c, regression='n')[1]

# Right way: Engle-Granger with the proper critical values
eg_stat, p_eg, eg_crit = coint(coint_df['y'], coint_df['x'], trend='c')
print(f"Plain ADF on residuals (wrong critical values): p = {p_wrong:.4f}")

Plain ADF on residuals (wrong critical values): p = 0.0000

print(f"Engle-Granger test (coint):                     p = {p_eg:.4f}")

Engle-Granger test (coint):                     p = 0.0000

H_0 of the Engle-Granger test: no cointegration. If the p-value is below 0.05 we conclude that the series are cointegrated. For this simulated pair, cointegration is detected clearly. For the independent random walks of the spurious example it should not be:

p_eg_spur = coint(sim_df['y'], sim_df['x'], trend='c')[1]
print(f"Engle-Granger p-value, independent random walks: {p_eg_spur:.4f}")

Engle-Granger p-value, independent random walks: 0.4547

What about infant mortality and exports?

p_eg_wb = coint(wb_log['exports'], wb_log['infant_mortality'], trend='c')[1]
wb_verdict = "are cointegrated, so the regression in levels could be interpreted as a long-run relationship" if p_eg_wb < 0.05 else "are NOT cointegrated, so the regression in levels is spurious"
print(f"Engle-Granger p-value: {p_eg_wb:.4f}")

Engle-Granger p-value: 0.9859

The p-value is 0.9859: according to the test, the two series are NOT cointegrated, so the regression in levels is spurious. With only 20 annual observations the test has low power, so we combine it with the economic logic.

5.4 A real cointegrated pair: WTI vs. Brent oil prices

WTI (West Texas Intermediate, US) and Brent (North Sea) are two benchmark crude oils. They are close substitutes, so if one becomes much cheaper than the other, buyers switch and arbitrage tends to close the gap. That suggests cointegration. We download from Yahoo Finance the daily prices of the front-month futures contracts of WTI (CL=F) and Brent (BZ=F), from 2008 to the end of 2025, and compute monthly averages.

The helper function get_yahoo downloads one ticker, tries again if Yahoo does not respond, and saves a local copy in the folder yahoo_data/, so the next time you run the code it does not need the internet:

import os
import time
import yfinance as yf

def get_yahoo(ticker, start='2000-01-01', end='2026-01-01', cache_dir='yahoo_data', retries=3):
    """Download daily closing prices from Yahoo Finance as a pandas Series.

    - `end` is NOT included (yfinance convention).
    - If Yahoo returns no data, it waits and tries again.
    - The first successful download is saved in the folder `yahoo_data/`; later renders read
      that file and do not need the internet. Delete the folder to download fresh data.
    """
    os.makedirs(cache_dir, exist_ok=True)
    name = ticker.replace('^', '').replace('=', '_')
    cache_file = os.path.join(cache_dir, f"{name}_{start}_{end}.csv")
    if os.path.exists(cache_file):
        return pd.read_csv(cache_file, index_col=0, parse_dates=True).iloc[:, 0].rename(ticker)
    for attempt in range(1, retries + 1):
        try:
            s = yf.download(ticker, start=start, end=end, auto_adjust=True, progress=False)['Close']
            s = s.squeeze('columns').dropna() if isinstance(s, pd.DataFrame) else s.dropna()
            if len(s) > 0:
                s = s.rename(ticker)
                s.index.name = 'date'
                s.to_csv(cache_file)
                return s
            print(f"Attempt {attempt}/{retries}: Yahoo returned no data for {ticker}")
        except Exception as e:
            print(f"Attempt {attempt}/{retries} for {ticker} failed: {type(e).__name__}")
        time.sleep(5 * attempt)
    raise RuntimeError(f"Could not download {ticker} from Yahoo Finance. Try again later.")

oil_daily = pd.DataFrame({'wti': get_yahoo('CL=F', start='2008-01-01', end='2026-01-01'),
                          'brent': get_yahoo('BZ=F', start='2008-01-01', end='2026-01-01')})

We use monthly averages. Using averages also avoids the day of April 20, 2020 when the WTI futures price closed negative (about -37 dollars per barrel), which has no logarithm.

oil = oil_daily.resample('MS').mean().loc[:'2025-12-01'].dropna()
ln_oil = np.log(oil)

plt.figure(figsize=(12, 4))

<Figure size 1200x400 with 0 Axes>

plt.plot(ln_oil.index, ln_oil['wti'], label='log WTI')

[<matplotlib.lines.Line2D object at 0x000001D6ABA956A0>]

plt.plot(ln_oil.index, ln_oil['brent'], label='log Brent')

[<matplotlib.lines.Line2D object at 0x000001D6ABA957F0>]

plt.title('Monthly average oil prices (logs)')

Text(0.5, 1.0, 'Monthly average oil prices (logs)')

plt.legend()

<matplotlib.legend.Legend object at 0x000001D6ABA95400>

plt.grid(True)
plt.show()

Step 0: are both series I(1)?

oil_i1 = pd.DataFrame({
    'p-value ADF level (c)': [adf_test(ln_oil[c], 'c')['p-value'] for c in ['wti', 'brent']],
    'p-value ADF 1st diff (c)': [adf_test(ln_oil[c].diff(), 'c')['p-value'] for c in ['wti', 'brent']],
}, index=['log WTI', 'log Brent']).round(4)
oil_i1

           p-value ADF level (c)  p-value ADF 1st diff (c)
log WTI                   0.0224                       0.0
log Brent                 0.0428                       0.0

We expect high p-values in levels (unit root not rejected) and very low p-values in differences.

Step 1: the long-run regression:

oil_lr = smf.ols('brent ~ wti', data=ln_oil).fit()
beta_oil = float(oil_lr.params['wti'])
print(oil_lr.summary().tables[1])

==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      0.1615      0.054      2.976      0.003       0.055       0.268
wti            0.9789      0.013     76.613      0.000       0.954       1.004
==============================================================================

The estimated long-run elasticity is \hat\beta_1= 0.979: a 1% higher WTI price goes with about 0.98% higher Brent price in the long run.

Step 2: is the deviation from equilibrium (the spread) stationary?

spread = oil_lr.resid

plt.figure(figsize=(12, 4))

<Figure size 1200x400 with 0 Axes>

plt.plot(spread.index, spread)

[<matplotlib.lines.Line2D object at 0x000001D6ABB72CF0>]

plt.axhline(0, color='black', linewidth=0.8)

<matplotlib.lines.Line2D object at 0x000001D6ABB72A50>

plt.axvspan(pd.Timestamp('2011-01-01'), pd.Timestamp('2014-12-01'), color='orange', alpha=0.2,
            label='2011-2014: US pipeline bottleneck at Cushing')

<matplotlib.patches.Rectangle object at 0x000001D6ABB7FED0>

plt.title('Deviation from the long-run relationship (residuals)')

Text(0.5, 1.0, 'Deviation from the long-run relationship (residuals)')

plt.legend()

<matplotlib.legend.Legend object at 0x000001D6ABB72BA0>

plt.grid(True)
plt.show()

p_eg_oil = coint(ln_oil['brent'], ln_oil['wti'], trend='c')[1]
oil_verdict = "reject" if p_eg_oil < 0.05 else "cannot reject"
print(f"Engle-Granger p-value: {p_eg_oil:.4f}")

Engle-Granger p-value: 0.0143

The Engle-Granger p-value is 0.0143, so we reject the null hypothesis of no cointegration at the 5% level.

Look at the residual plot. Between 2011 and 2014 the spread moved far from zero for a long time: US oil production grew fast and could not reach the coast because of pipeline bottlenecks, so WTI became much cheaper than Brent. This kind of structural break makes the deviation look less stationary and reduces the power of the test. Statistical tests do not replace the knowledge of the business.

5.4.1 Error correction: how the gap closes

If the series are cointegrated, deviations from equilibrium should predict future changes that close the gap. We check this with an error-correction model (ECM):

\Delta \ln Brent_t = \alpha + \gamma\,\hat u_{t-1} + \delta\,\Delta \ln WTI_t + e_t

If \gamma<0, when Brent is above its equilibrium with WTI (\hat u_{t-1}>0), Brent tends to fall in the next period, so the gap closes.

ecm_df = pd.DataFrame({'d_brent': ln_oil['brent'].diff(),
                       'd_wti': ln_oil['wti'].diff(),
                       'u_lag': spread.shift(1)}).dropna()
ecm = smf.ols('d_brent ~ u_lag + d_wti', data=ecm_df).fit()
gamma = float(ecm.params['u_lag'])
gamma_p = float(ecm.pvalues['u_lag'])
print(ecm.summary().tables[1])

==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     -0.0001      0.002     -0.054      0.957      -0.005       0.004
u_lag         -0.1742      0.038     -4.558      0.000      -0.250      -0.099
d_wti          0.7870      0.022     36.603      0.000       0.745       0.829
==============================================================================

The error-correction coefficient is \hat\gamma= -0.174 (p-value = 0). A negative and significant \gamma means that each month about 17% of the deviation from equilibrium is corrected.

5.5 Why this matters for AI and machine learning

These ideas apply to any model trained on time series, not only to linear regression:

Spurious features. Two trending variables are always highly correlated. A model that uses one to predict the other can look excellent in-sample and fail out of sample. Correlation screening for feature selection is especially dangerous with trending data.
Distribution shift. A non-stationary series changes its mean (and often its variance) over time. The data in the test period comes from a different distribution than the training data, which violates the basic assumption of most ML models.
Tree-based models cannot extrapolate. A random forest or gradient-boosting model predicts values inside the range it saw in training. If sales keep growing, the model will predict a flat line at the training maximum. Modeling growth rates (stationary) and converting back to levels solves this.
Chronological splits. Train/test splits and cross-validation must respect time order; random shuffling leaks future information into training.
Cointegration = useful structure. When series share a long-run equilibrium (prices of substitutes, spot and futures prices), the deviation from equilibrium is a powerful, stationary feature. This is the basis of pairs trading in finance.

6 CHALLENGE 2 (Part B)

Spurious regression experiment. Extend spurious_rate to compute the rejection rate (p<0.05) for n = 25, 50, 100, 250, 500 and 1000 in three cases:
1. two independent random walks, in levels;
2. the same random walks, in first differences;
3. two independent stationary AR(1) series with \phi_1=0.5 (use simulate_ar1), in levels.
Plot the three rejection rates against n. Explain in your own words why case (a) gets worse with n, and why case (c) is much better than (a) even though the series are autocorrelated.
Cointegration with real data. Download from Yahoo Finance (with get_yahoo) the daily prices of the RBOB gasoline futures (RB=F, dollars per gallon) and the WTI futures (CL=F, dollars per barrel) from 2006-01-01 to 2025-12-31. Convert both to monthly averages and take logs.
1. Show that both series are I(1).
2. Run the Engle-Granger test with coint. Are they cointegrated?
3. Estimate the error-correction model for gasoline. Interpret \hat\gamma: how many months does it take for half of a deviation to disappear? (Hint: half-life =\ln(0.5)/\ln(1+\hat\gamma).)
A pair that should not be cointegrated. Repeat the Engle-Granger test for log WTI vs. the 10-year US Treasury yield (^TNX in Yahoo Finance, in %, monthly average, 2006–2025). Run the regression in levels and in differences, and decide whether the regression in levels is spurious. Justify with your test results.

What I will look for: correct code and results; correct choice of test specifications; interpretations that use your numbers; and clear explanations in your own words.

7 PART C. ARIMA models: the building blocks

We import the additional tools for this part:

import os
import warnings
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.tsa.stattools import pacf
from statsmodels.graphics.tsaplots import plot_pacf
from statsmodels.stats.diagnostic import acorr_ljungbox

warnings.filterwarnings('ignore')   # hide convergence warnings while fitting many models

7.1 The Box-Jenkins methodology

The ARIMA/SARIMA family of models was popularized by Box and Jenkins (1970). The idea is to explain a time series Y_t using only its own past: its past values and its past random shocks. The methodology is an iterative cycle:

Identification: find a stationary version of the series and look at its autocorrelations to propose a model.
Estimation: estimate the coefficients.
Diagnostics: check whether the residuals look like white noise (no autocorrelation left). If not, go back to step 1 and adjust the model.
Forecasting: use the final model to predict, and evaluate the forecasts on data the model did not see.

7.2 AR(p): the series depends on its own past values

In Part A we studied the AR(1) model. An autoregressive model of order p includes p past values (lags) of the series:

Y_t=\phi_0+\phi_1Y_{t-1}+\phi_2Y_{t-2}+\dots+\phi_pY_{t-p}+\varepsilon_t

A shock \varepsilon_t enters Y_t; Y_t enters Y_{t+1}; Y_{t+1} enters Y_{t+2}; and so on. A shock is never completely forgotten, but in a stationary AR model its effect fades geometrically (in an AR(1), a fraction \phi_1^k remains after k periods). We say that an AR process has an infinite but fading memory.

7.3 MA(q): the series depends on past shocks

A moving-average model of order q explains Y_t with the current shock and the last q shocks:

Y_t=\theta_0+\varepsilon_t+\theta_1\varepsilon_{t-1}+\dots+\theta_q\varepsilon_{t-q}

Here a shock affects the series for exactly q periods and then it is completely forgotten: an MA process has a finite memory of q periods. Example: an unexpected promotion boosts sales this month and partially next month, and then its effect disappears.

7.4 ARMA(p,q), ARIMA(p,d,q) and SARIMA

ARMA(p,q) combines both types of terms:

Y_t=\phi_0+\sum_{k=1}^{p}\phi_kY_{t-k}+\sum_{k=1}^{q}\theta_k\varepsilon_{t-k}+\varepsilon_t

ARMA models require a stationary series. If we have to difference the series d times to make it stationary, the model is called ARIMA(p,d,q), where the “I” stands for integrated.

When the data are monthly or quarterly, we often also need seasonal terms: lags of 12 months (or 4 quarters) and their multiples. The full notation is

\text{SARIMA}(p,d,q)(P,D,Q)_s

Parameter	Meaning
p	number of AR terms (lags 1, 2, …, p)
d	number of first differences (Y_t-Y_{t-1})
q	number of MA terms (shocks of t-1, …, t-q)
P	number of seasonal AR terms (lags s, 2s, …)
D	number of seasonal differences (Y_t-Y_{t-s})
Q	number of seasonal MA terms (shocks of t-s, t-2s, …)
s	number of periods in one year (12 for monthly, 4 for quarterly)

If we add external explanatory variables (X’s), the model becomes an ARIMAX / SARIMAX model. That is why the Python class is called SARIMAX.

7.5 Identification tools: ACF and PACF

7.5.1 The autocorrelation function (ACF)

The autocorrelation at lag k is the correlation between Y_t and Y_{t-k}. Since a stationary series has the same variance in every period, it simplifies to

\rho_k=\frac{Cov(Y_t,Y_{t-k})}{SD(Y_t)\,SD(Y_{t-k})}=\frac{\gamma(k)}{\sigma_Y^2}

where \gamma(k) is the autocovariance at lag k.

7.5.2 The partial autocorrelation function (PACF)

The partial autocorrelation at lag k measures the relationship between Y_t and Y_{t-k} after removing the effect of the intermediate lags Y_{t-1},\dots,Y_{t-k+1}. One way to compute it is the coefficient of Y_{t-k} in the regression

Y_t=\alpha+\beta_1Y_{t-1}+\dots+\beta_kY_{t-k}+e_t \qquad\Rightarrow\qquad PACF(k)=\hat\beta_k

Let’s check this with a simulated AR(1) series. In an AR(1), Y_{t-2} affects Y_t only through Y_{t-1}, so the PACF at lag 2 should be about zero while the ACF at lag 2 is not:

rng_seed = 42
ar1 = ArmaProcess(ar=[1, -0.7], ma=[1]).generate_sample(nsample=500, distrvs=np.random.default_rng(rng_seed).standard_normal)
ar1 = pd.Series(ar1)

lags_df = pd.DataFrame({'y': ar1, 'lag1': ar1.shift(1), 'lag2': ar1.shift(2)}).dropna()
beta2 = sm.OLS(lags_df['y'], sm.add_constant(lags_df[['lag1', 'lag2']])).fit().params['lag2']

print(f"ACF at lag 2:                       {ar1.autocorr(2):.3f}   (theory: 0.7^2 = 0.49)")

ACF at lag 2:                       0.459   (theory: 0.7^2 = 0.49)

print(f"PACF at lag 2 (statsmodels):        {pacf(ar1, nlags=2, method='ols')[2]:.3f}")

PACF at lag 2 (statsmodels):        -0.112

print(f"PACF at lag 2 (regression coefficient): {beta2:.3f}")

PACF at lag 2 (regression coefficient): -0.112

Note: ArmaProcess writes the AR polynomial as 1-\phi_1L, so \phi_1=0.7 is entered as ar=[1, -0.7]. The MA polynomial is 1+\theta_1L, so \theta_1 is entered with its own sign: ma=[1, θ1].

7.5.3 Patterns (signatures) of AR, MA and ARMA processes

Let’s simulate four processes where we know the true model, and look at their ACF and PACF:

processes = {
    'AR(1): phi1=0.7':           ArmaProcess(ar=[1, -0.7], ma=[1]),
    'AR(2): phi1=0.5, phi2=0.3': ArmaProcess(ar=[1, -0.5, -0.3], ma=[1]),
    'MA(1): theta1=0.7':         ArmaProcess(ar=[1], ma=[1, 0.7]),
    'ARMA(1,1): phi1=0.7, theta1=0.4': ArmaProcess(ar=[1, -0.7], ma=[1, 0.4]),
}
sim_series = {name: proc.generate_sample(nsample=500, distrvs=np.random.default_rng(1).standard_normal)
              for name, proc in processes.items()}

fig, axes = plt.subplots(4, 2, figsize=(12, 13))
for i, (name, y) in enumerate(sim_series.items()):
    plot_acf(y, lags=20, zero=False, ax=axes[i, 0], title=f'ACF  - {name}')
    plot_pacf(y, lags=20, zero=False, ax=axes[i, 1], title=f'PACF - {name}', method='ywm')

<Figure size 1200x1300 with 8 Axes>
<Figure size 1200x1300 with 8 Axes>
<Figure size 1200x1300 with 8 Axes>
<Figure size 1200x1300 with 8 Axes>
<Figure size 1200x1300 with 8 Axes>
<Figure size 1200x1300 with 8 Axes>
<Figure size 1200x1300 with 8 Axes>
<Figure size 1200x1300 with 8 Axes>

plt.tight_layout()
plt.show()

The shaded area is the 95% confidence band: spikes outside it are significantly different from zero. The patterns are:

True model	ACF	PACF
AR(p)	decays gradually	cuts off after lag p
MA(q)	cuts off after lag q	decays gradually (often alternating sign)
ARMA(p,q)	decays gradually	decays gradually

So the PACF suggests p and the ACF suggests q. With real data the patterns are noisier, so these rules give us a starting point, not the final answer.

7.5.4 In-class exercise: identify the mystery series

Run the next chunk without reading the parameters (they are hidden in a list). Look at the ACF and PACF of each mystery series and propose a model:

_secret = [([1, -0.8], [1]), ([1], [1, -0.6]), ([1, -0.6, 0.3], [1])]
mystery = {f'Mystery {k+1}': ArmaProcess(*_secret[k]).generate_sample(400, distrvs=np.random.default_rng(100 + k).standard_normal)
           for k in range(3)}

fig, axes = plt.subplots(3, 2, figsize=(12, 9))
for i, (name, y) in enumerate(mystery.items()):
    plot_acf(y, lags=20, zero=False, ax=axes[i, 0], title=f'ACF  - {name}')
    plot_pacf(y, lags=20, zero=False, ax=axes[i, 1], title=f'PACF - {name}', method='ywm')

<Figure size 1200x900 with 6 Axes>
<Figure size 1200x900 with 6 Axes>
<Figure size 1200x900 with 6 Axes>
<Figure size 1200x900 with 6 Axes>
<Figure size 1200x900 with 6 Axes>
<Figure size 1200x900 with 6 Axes>

plt.tight_layout()
plt.show()

Now estimate the model you proposed for each series and check whether the estimated coefficients are significant. For example, if you think Mystery 1 is an AR(1):

res_m1 = SARIMAX(mystery['Mystery 1'], order=(1, 0, 0), trend='c').fit(disp=False)
print(res_m1.summary().tables[1])

==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept     -0.0066      0.049     -0.133      0.894      -0.103       0.090
ar.L1          0.7802      0.030     26.397      0.000       0.722       0.838
sigma2         0.9722      0.075     12.976      0.000       0.825       1.119
==============================================================================

The coefficients are estimated by Maximum Likelihood (ML). Unlike OLS, ML has no closed formula for ARMA models: the software searches iteratively for the coefficients that make the observed data most likely. MA terms make OLS impossible because the past shocks \varepsilon_{t-k} are not observed; they must be inferred together with the coefficients.

7.6 Real data: the IGAE of Mexico

The IGAE (Indicador Global de la Actividad Económica) is published monthly by INEGI. It tracks the economic activity of Mexico, similar to the GDP, but it is published every month and sooner than the quarterly GDP. We use the series without seasonal adjustment, because we want to model the seasonality ourselves.

7.6.1 Downloading the IGAE with the INEGI API

INEGI offers a free API to download its indicators. We use the Python package INEGIpy, which simplifies the calls to this API. Install it first (only once):

pip install INEGIpy==2.0.0

The API requires a token. For now you can use my token, which is in the code below. You can request your own token for free at https://www.inegi.org.mx/servicios/api_indicadores.html. Never publish a token, for example in a public GitHub repository.

Note: in December 2025 INEGI changed the ID of the IGAE. The ID that we use, 6207136901, is the IGAE index (base 2018=100), original figures (not seasonally adjusted), total of the economy.

# pip install INEGIpy==2.0.0
from INEGIpy import Indicadores

# You can ask for your own token at the INEGI web page. Here is my token for now:
TOKEN = 'c4559227-83d3-aeb4-1a8d-bee73f627140'
inegi = Indicadores(TOKEN)

# IGAE mensual: Índice base 2018=100, cifras originales, Total
# (INEGI changed this ID in Dec 2025)
IGAE_ID = "6207136901"
igaedf = inegi.obtener_df(indicadores=IGAE_ID, nombres="IGAE")

obtener_df returns a data frame with the dates as index and one column called IGAE. We make sure that the index is a date type, that each date is the first day of its month, and that the values are numeric:

igaedf.index = pd.to_datetime(igaedf.index)
igaedf.index = igaedf.index.to_period('M').to_timestamp()   # first day of each month
igaedf = igaedf[['IGAE']].apply(pd.to_numeric, errors='coerce')
igaedf = igaedf.sort_index().asfreq('MS')
igaedf['logIGAE'] = np.log(igaedf['IGAE'])
igaedf.tail()

                  IGAE   logIGAE
fechas                          
2026-04-01  106.613699  4.669212
2026-05-01  108.549177  4.687203
2026-06-01  107.390831  4.676475
2026-07-01  109.262414  4.693752
2026-08-01  108.317094  4.685063

fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)
axes[0].plot(igaedf.index, igaedf['IGAE'])

[<matplotlib.lines.Line2D object at 0x000001D6ACD163C0>]

axes[0].set_title('IGAE (original units)')

Text(0.5, 1.0, 'IGAE (original units)')

axes[1].plot(igaedf.index, igaedf['logIGAE'])

[<matplotlib.lines.Line2D object at 0x000001D6ACD16510>]

axes[1].set_title('log IGAE')

Text(0.5, 1.0, 'log IGAE')

for ax in axes:
    for year in [1994, 2000, 2006, 2012, 2018, 2024]:
        ax.axvline(pd.Timestamp(f'{year}-12-01'), color='green', linestyle='--', linewidth=0.8)
    ax.grid(True)

<matplotlib.lines.Line2D object at 0x000001D6ACD16120>
<matplotlib.lines.Line2D object at 0x000001D6ACD16270>
<matplotlib.lines.Line2D object at 0x000001D6ACD16660>
<matplotlib.lines.Line2D object at 0x000001D6ACD167B0>
<matplotlib.lines.Line2D object at 0x000001D6ACD16900>
<matplotlib.lines.Line2D object at 0x000001D6ACD16A50>
<matplotlib.lines.Line2D object at 0x000001D6ACD16BA0>
<matplotlib.lines.Line2D object at 0x000001D6ACD16CF0>
<matplotlib.lines.Line2D object at 0x000001D6ACD16E40>
<matplotlib.lines.Line2D object at 0x000001D6ACD16F90>
<matplotlib.lines.Line2D object at 0x000001D6ACD170E0>
<matplotlib.lines.Line2D object at 0x000001D6ACD17230>

plt.tight_layout()
plt.show()

The green lines mark the start of each 6-year federal administration. We see a growing trend, a strong seasonal pattern, and the crises of 1995, 2009 and 2020. On a log scale, equal vertical distances mean equal % changes, so the log plot is the right one to compare the depth of crises from different decades.

How deep was the COVID crisis? Because the data are not seasonally adjusted, we compare each month with the same month of the previous year (the annual % change):

annual_log_change = igaedf['logIGAE'].diff(12)
worst_month = annual_log_change.idxmin()
worst_log = float(annual_log_change.min())
worst_simple = float(np.exp(worst_log) - 1)
worst_label = worst_month.strftime('%B %Y')
print(f"Worst annual change: {worst_label}: log difference = {worst_log:.3f}, simple % change = {100*worst_simple:.1f}%")

Worst annual change: May 2020: log difference = -0.276, simple % change = -24.1%

The worst month was May 2020, with an annual log change of -0.276, which is a simple % change of -24.1%. Remember from Part A that log changes and simple % changes are close only for small changes.

7.6.2 Step 1 – Identification: finding a stationary version

For monthly economic series we start with the seasonal difference of the log (the annual % change), because it removes both the trend and the seasonality. We compare it with the log level using the ADF test:

igaedf['annualgrowth'] = igaedf['logIGAE'].diff(12)

p_adf_level = adf_test(igaedf['logIGAE'], 'ct')['p-value']        # adf_test was defined in Part A
p_adf_annual = adf_test(igaedf['annualgrowth'], 'c')['p-value']
igae_stationary = "stationary" if p_adf_annual < 0.05 else "NOT stationary (consider an additional first difference, d=1)"
print(f"ADF p-value, log level (ct):     {p_adf_level:.4f}")

ADF p-value, log level (ct):     0.2475

print(f"ADF p-value, annual growth (c):  {p_adf_annual:.4f}")

ADF p-value, annual growth (c):  0.0011

plt.figure(figsize=(12, 4))

<Figure size 1200x400 with 0 Axes>

plt.plot(igaedf.index, 100 * igaedf['annualgrowth'])

[<matplotlib.lines.Line2D object at 0x000001D6AC978440>]

plt.axhline(0, color='black', linewidth=0.8)

<matplotlib.lines.Line2D object at 0x000001D6AC9786E0>

plt.title('IGAE annual % change (seasonal difference of the log)')

Text(0.5, 1.0, 'IGAE annual % change (seasonal difference of the log)')

plt.grid(True)
plt.show()

For the annual % change, the ADF p-value is 0.0011, so we treat the series as stationary.

This defines the differencing parameters. We model the log of the IGAE with one seasonal difference and no first difference:

d=0,\qquad D=1,\qquad s=12

7.6.3 Step 1 (cont.) – ACF and PACF of the stationary series

ag = igaedf['annualgrowth'].dropna()
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(ag, lags=36, zero=False, ax=axes[0], title='ACF - IGAE annual growth')

<Figure size 1200x400 with 2 Axes>

plot_pacf(ag, lags=36, zero=False, ax=axes[1], title='PACF - IGAE annual growth', method='ywm')

<Figure size 1200x400 with 2 Axes>

plt.tight_layout()
plt.show()


acf1 = float(ag.autocorr(1))
pacf_vals, pacf_ci = pacf(ag, nlags=24, alpha=0.05, method='ywm')
signif = (pacf_ci[:, 0] > 0) | (pacf_ci[:, 1] < 0)
# number of consecutive significant PACF lags starting at lag 1:
p_suggest = int(np.argmin(signif[1:])) if not signif[1:].all() else 24

The ACF decays slowly from a lag-1 autocorrelation of 0.77: the annual growth is very persistent (if the economy is growing faster than a year ago this month, it probably will be next month too). The PACF has 1 significant spike(s) at the first lags and then falls to near zero. This is an AR signature, so we start with a small p and no MA terms. As a first simple model we use p=1, and let the diagnostics tell us what is missing.

8 PART D. Building and evaluating a SARIMA model

8.1 Train / test split

To evaluate forecasts honestly, we hide the last 24 months and use them only at the end. The split must be chronological. A random split would let the model “see the future” during training.

igae_log = igaedf['logIGAE'].dropna()
h_test = 24
train = igae_log.iloc[:-h_test]
test = igae_log.iloc[-h_test:]
print(f"Training: {train.index[0]:%Y-%m} to {train.index[-1]:%Y-%m} ({len(train)} months)")

Training: 1993-01 to 2024-08 (380 months)

print(f"Test:     {test.index[0]:%Y-%m} to {test.index[-1]:%Y-%m} ({len(test)} months)")

Test:     2024-09 to 2026-08 (24 months)

Two helper functions: one to fit a SARIMA model, and one to run the diagnostics.

def fit_sarima(y, order, seasonal_order, trend='c', exog=None):
    """Fit a SARIMA(X) model by maximum likelihood."""
    return SARIMAX(y, order=order, seasonal_order=seasonal_order,
                   trend=trend, exog=exog).fit(disp=False)

def diagnose(res, name, plots=True):
    """Coefficient table, Ljung-Box test on residuals and AIC; residual ACF/PACF plots."""
    s = res.model.seasonal_periods or 1
    burn = res.model.k_diff + s * res.model.k_seasonal_diff + 1   # drop the first residuals (initialization)
    resid = res.resid.iloc[burn:]
    m = res.model
    k = m.k_ar_params + m.k_ma_params + m.k_seasonal_ar_params + m.k_seasonal_ma_params   # number of ARMA terms
    lb = acorr_ljungbox(resid, lags=[12, 24], model_df=max(k, 0))
    coefs = pd.DataFrame({'coef': res.params, 'p-value': res.pvalues}).drop('sigma2').round(4)
    print(f"===== {name} =====")
    print(coefs)
    print(f"Ljung-Box p-value (lag 12): {lb['lb_pvalue'].iloc[0]:.4f} | (lag 24): {lb['lb_pvalue'].iloc[1]:.4f}")
    print(f"AIC = {res.aic:.2f}")
    if plots:
        fig, axes = plt.subplots(1, 2, figsize=(12, 3.5))
        plot_acf(resid, lags=36, zero=False, ax=axes[0], title=f'Residual ACF - {name}')
        plot_pacf(resid, lags=36, zero=False, ax=axes[1], title=f'Residual PACF - {name}', method='ywm')
        plt.tight_layout()
        plt.show()
    return resid, lb

The Ljung-Box test (Ljung and Box 1978) checks whether a group of autocorrelations of the residuals is jointly zero. H_0: the residuals are white noise (no autocorrelation left). We want a high p-value (above 0.05). A low p-value means that the model is leaving predictable structure in the errors.

8.2 Step 2 & 3 – Estimation and diagnostics: calibrating the model

8.2.1 Model 1: SARIMA(1,0,0)(0,1,0)₁₂

res1 = fit_sarima(train, (1, 0, 0), (0, 1, 0, 12))
resid1, lb1 = diagnose(res1, 'Model 1: (1,0,0)(0,1,0,12)')

===== Model 1: (1,0,0)(0,1,0,12) =====
             coef  p-value
intercept  0.0045   0.0038
ar.L1      0.7676   0.0000
Ljung-Box p-value (lag 12): 0.0000 | (lag 24): 0.0000
AIC = -1536.81

res1_acf12 = float(resid1.autocorr(12))

The AR(1) coefficient is significant. But look at the residuals: the autocorrelation at lag 12 is -0.32. A significant negative spike at the seasonal lag after a seasonal difference is the classic signal for a seasonal MA term (Q=1). It means that the seasonal difference “overcorrects”: part of the shock of 12 months ago is undone this year.

8.2.2 Model 2: SARIMA(1,0,0)(0,1,1)₁₂

res2 = fit_sarima(train, (1, 0, 0), (0, 1, 1, 12))
resid2, lb2 = diagnose(res2, 'Model 2: (1,0,0)(0,1,1,12)')

===== Model 2: (1,0,0)(0,1,1,12) =====
             coef  p-value
intercept  0.0027      0.0
ar.L1      0.8600      0.0
ma.S.L12  -0.8221      0.0
Ljung-Box p-value (lag 12): 0.0000 | (lag 24): 0.0000
AIC = -1656.18

res2_acf1 = float(resid2.autocorr(1))

The seasonal MA term is significant and the lag-12 spike disappears. Now look at lag 1: its residual autocorrelation is -0.15. A remaining significant spike at lag 1 (with an AR(1) already in the model) suggests adding an MA(1) term (q=1).

8.2.3 Model 3: SARIMA(1,0,1)(0,1,1)₁₂

res3 = fit_sarima(train, (1, 0, 1), (0, 1, 1, 12))
resid3, lb3 = diagnose(res3, 'Model 3: (1,0,1)(0,1,1,12)')

===== Model 3: (1,0,1)(0,1,1,12) =====
             coef  p-value
intercept  0.0014   0.0012
ar.L1      0.9253   0.0000
ma.L1     -0.2607   0.0000
ma.S.L12  -0.8293   0.0000
Ljung-Box p-value (lag 12): 0.0006 | (lag 24): 0.0000
AIC = -1669.05

lb3_p24 = float(lb3['lb_pvalue'].iloc[1])

The Ljung-Box p-value at lag 24 is 0. If it is above 0.05, we cannot reject that the residuals are white noise. One or two isolated spikes slightly outside the band are expected by chance: with 36 lags and a 95% band, about 2 spikes will be “significant” even for pure white noise.

8.2.4 Confirming with an information criterion (AIC)

Calibrating by hand is a good way to learn, but we should check alternatives in a systematic way. The AIC (Akaike Information Criterion) rewards the fit of a model but penalizes the number of parameters: lower is better. We fit all combinations of small orders and rank them:

from itertools import product

def sarima_grid(y, d, D, s, max_p=2, max_q=2, max_P=1, max_Q=1, trend='c'):
    """Fit all SARIMA(p,d,q)(P,D,Q,s) combinations up to the given orders; rank by AIC."""
    rows = []
    for p, q, P, Q in product(range(max_p + 1), range(max_q + 1), range(max_P + 1), range(max_Q + 1)):
        try:
            res = fit_sarima(y, (p, d, q), (P, D, Q, s), trend=trend)
            burn = d + s * D + 1
            k = p + q + P + Q
            lb_p = acorr_ljungbox(res.resid.iloc[burn:], lags=[2 * s], model_df=k)['lb_pvalue'].iloc[0]
            rows.append({'order': (p, d, q), 'seasonal': (P, D, Q, s), 'AIC': res.aic,
                         'BIC': res.bic, f'Ljung-Box p (lag {2*s})': lb_p, 'n params': k})
        except Exception:
            pass
    return pd.DataFrame(rows).sort_values('AIC').reset_index(drop=True)

grid = sarima_grid(train, d=0, D=1, s=12)
grid.head(8).round(3)

       order       seasonal       AIC       BIC  Ljung-Box p (lag 24)  n params
0  (1, 0, 1)  (1, 1, 1, 12) -1670.167 -1646.718                   0.0         4
1  (2, 0, 2)  (1, 1, 1, 12) -1669.651 -1638.386                   0.0         6
2  (2, 0, 0)  (1, 1, 1, 12) -1669.252 -1645.803                   0.0         4
3  (1, 0, 1)  (0, 1, 1, 12) -1669.046 -1649.506                   0.0         3
4  (2, 0, 2)  (0, 1, 1, 12) -1668.887 -1641.530                   0.0         5
5  (1, 0, 2)  (1, 1, 1, 12) -1668.598 -1641.242                   0.0         5
6  (2, 0, 1)  (0, 1, 1, 12) -1666.479 -1643.030                   0.0         4
7  (2, 0, 1)  (1, 1, 1, 12) -1665.462 -1638.105                   0.0         5

best_order = grid.loc[0, 'order']
best_seasonal = grid.loc[0, 'seasonal']
best_label = f"SARIMA{best_order}{best_seasonal[:3]}_12"
aic_model3 = float(res3.aic)
aic_best = float(grid.loc[0, 'AIC'])

The model with the lowest AIC is SARIMA(1, 0, 1)(1, 1, 1)_12 (AIC = -1670.2), while our hand-calibrated Model 3 has AIC = -1669. When two models have similar AIC (a difference of less than about 2–4 points), we prefer the simpler one, and we only consider models whose residuals pass the Ljung-Box test. Below, we keep Model 3 for interpretation and compare both models in the test period.

8.3 Interpreting the model

Let w_t=\ln IGAE_t-\ln IGAE_{t-12} be the annual growth. In statsmodels notation, Model 3 is

(1-\phi_1L)\,w_t = c + (1+\theta_1L)(1+\Theta_1L^{12})\,\varepsilon_t

where L is the lag operator (L^kY_t=Y_{t-k}). Expanding the product of the MA polynomials:

w_t=c+\phi_1w_{t-1}+\varepsilon_t+\theta_1\varepsilon_{t-1}+\Theta_1\varepsilon_{t-12}+\theta_1\Theta_1\varepsilon_{t-13}

The last term appears because the regular and seasonal MA terms multiply. That is why it is called a multiplicative seasonal model.

c_hat = float(res3.params['intercept'])
phi1_hat = float(res3.params['ar.L1'])
theta1_hat = float(res3.params['ma.L1'])
Theta1_hat = float(res3.params['ma.S.L12'])
mean_growth = c_hat / (1 - phi1_hat)                 # long-run mean of the annual growth
half_life = np.log(0.5) / np.log(abs(phi1_hat))      # months
print(f"c = {c_hat:.5f}, phi1 = {phi1_hat:.4f}, theta1 = {theta1_hat:.4f}, Theta1 = {Theta1_hat:.4f}")

c = 0.00142, phi1 = 0.9253, theta1 = -0.2607, Theta1 = -0.8293

print(f"Long-run mean annual growth = {100*mean_growth:.2f}%   |   half-life of a shock to growth = {half_life:.1f} months")

Long-run mean annual growth = 1.90%   |   half-life of a shock to growth = 8.9 months

\hat\phi_1= 0.925: the annual growth is very persistent. Half of a shock to the growth rate disappears after about 8.9 months. If the economy grows above its average this month, it will most likely keep growing above average for the next months, and the same holds for bad months. For a stationary model we need |\phi_1|<1.
\hat c= 0.0014 is not the mean growth. As in the AR(1) of Part A, the long-run mean is c/(1-\phi_1)= 1.9% per year, the average annual growth of the Mexican economy implied by the model.
\hat\theta_1= -0.261: the effect of last month’s shock. A negative value means that part of last month’s surprise is reversed this month.
\hat\Theta_1= -0.829: the seasonal MA. A value close to -1 means that the seasonal pattern is very stable from one year to the next, so the seasonal difference overcorrects and the model undoes most of the shock of the same month last year.

8.4 Step 4 – Forecasting and out-of-sample evaluation

8.4.1 Benchmark: the seasonal naive forecast

Any model must beat a simple benchmark. For seasonal data the natural benchmark is the seasonal naive forecast: each month will be equal to the same month of the last observed year. It is the seasonal version of the random walk of Workshop 3.

def metrics(actual_log, pred_log):
    """RMSE and MAPE computed on the original scale of the IGAE."""
    a, p = np.exp(actual_log), np.exp(pred_log)
    return {'RMSE': float(np.sqrt(np.mean((a - p) ** 2))),
            'MAPE (%)': float(100 * np.mean(np.abs((a - p) / a)))}

# Seasonal naive: repeat the last 12 months of the training data
last_year = train.iloc[-12:].values
snaive = pd.Series(np.tile(last_year, int(np.ceil(h_test / 12)))[:h_test], index=test.index)

fc3 = res3.get_forecast(h_test)
res_best = fit_sarima(train, best_order, best_seasonal)
fc_best = res_best.get_forecast(h_test)

eval_table = pd.DataFrame({
    'Seasonal naive': metrics(test, snaive),
    'Model 3 (1,0,1)(0,1,1)': metrics(test, fc3.predicted_mean),
    f'Best AIC {best_label}': metrics(test, fc_best.predicted_mean),
}).T.round(3)
eval_table

                                       RMSE  MAPE (%)
Seasonal naive                        1.813     1.316
Model 3 (1,0,1)(0,1,1)                3.198     2.718
Best AIC SARIMA(1, 0, 1)(1, 1, 1)_12  3.672     3.248

ci3 = fc3.conf_int(alpha=0.05)
plt.figure(figsize=(12, 5))

<Figure size 1200x500 with 0 Axes>

plt.plot(np.exp(igae_log['2019-01-01':]), color='black', label='Actual IGAE')

[<matplotlib.lines.Line2D object at 0x000001D6ABC3B4D0>]

plt.plot(np.exp(fc3.predicted_mean), color='tab:red', label='Model 3 forecast')

[<matplotlib.lines.Line2D object at 0x000001D6ABC3B620>]

plt.fill_between(test.index, np.exp(ci3.iloc[:, 0]), np.exp(ci3.iloc[:, 1]), color='tab:red', alpha=0.15, label='95% interval')

<matplotlib.collections.FillBetweenPolyCollection object at 0x000001D6ABB40F50>

plt.plot(np.exp(snaive), color='tab:grey', linestyle='--', label='Seasonal naive')

[<matplotlib.lines.Line2D object at 0x000001D6ABC3B770>]

plt.axvline(test.index[0], color='black', linestyle=':')

<matplotlib.lines.Line2D object at 0x000001D6ABC3B230>

plt.title('Out-of-sample forecasts for the last 24 months')

Text(0.5, 1.0, 'Out-of-sample forecasts for the last 24 months')

plt.legend(loc='upper left')

<matplotlib.legend.Legend object at 0x000001D6ABC3B380>

plt.grid(True)
plt.show()

MAPE (mean absolute percentage error) is the average forecast error as a % of the actual value; RMSE (root mean squared error) is in IGAE points and penalizes large errors more. Discuss: does the SARIMA beat the seasonal naive benchmark? By how much?

8.4.2 One-step-ahead comparison with a machine-learning model

The previous evaluation forecasts up to 24 months ahead from one origin. A common alternative is one-step-ahead evaluation: each month we forecast only the next month, using all the information available up to the previous month. We compare:

Seasonal naive: \hat Y_t=Y_{t-12}.
SARIMA (Model 3), with the coefficients estimated only with the training data.
Gradient boosting (a tree-based ML model) that predicts the annual growth using the last 12 values of the annual growth as features.

Notice that we give the ML model the stationary version (the growth rate). Tree models cannot extrapolate: if we gave it the IGAE level, it could never predict a value above the maximum it saw in training (Part B).

from sklearn.ensemble import GradientBoostingRegressor

# SARIMA one-step-ahead: same coefficients (estimated on train), applied to the full series
res3_full = res3.apply(igae_log)
sarima_1step = res3_full.get_prediction(start=test.index[0], end=test.index[-1], dynamic=False).predicted_mean

# Seasonal naive one-step-ahead
snaive_1step = igae_log.shift(12).loc[test.index]

# Gradient boosting on lags of the annual growth
g = igae_log.diff(12)
ml_data = pd.concat({'g': g, **{f'lag{k}': g.shift(k) for k in range(1, 13)}}, axis=1).dropna()
X_train = ml_data.loc[:train.index[-1]].drop(columns='g')
y_train = ml_data.loc[:train.index[-1], 'g']
X_test = ml_data.loc[test.index].drop(columns='g')
gbm = GradientBoostingRegressor(n_estimators=300, max_depth=2, learning_rate=0.05, random_state=0)
gbm.fit(X_train, y_train)

GradientBoostingRegressor(learning_rate=0.05, max_depth=2, n_estimators=300,
                          random_state=0)

gbm_1step = pd.Series(gbm.predict(X_test), index=test.index) + igae_log.shift(12).loc[test.index]

onestep_table = pd.DataFrame({
    'Seasonal naive': metrics(test, snaive_1step),
    'SARIMA Model 3': metrics(test, sarima_1step),
    'Gradient boosting (lags of growth)': metrics(test, gbm_1step),
}).T.round(3)
onestep_table

                                     RMSE  MAPE (%)
Seasonal naive                      1.759     1.262
SARIMA Model 3                      1.488     1.130
Gradient boosting (lags of growth)  2.207     1.396

Discuss: Which model wins one step ahead? A more complex model is not automatically better. With a few hundred monthly observations, a well-specified SARIMA is often hard to beat. ML models usually shine when there are many related series (for example, thousands of products) or many useful external features.

8.4.3 Final forecast

For the final forecast we re-estimate the model with all the data and forecast 24 months ahead:

res_final = fit_sarima(igae_log, (1, 0, 1), (0, 1, 1, 12))
fc_final = res_final.get_forecast(24)
ci_final = fc_final.conf_int(alpha=0.05)

plt.figure(figsize=(12, 5))

<Figure size 1200x500 with 0 Axes>

plt.plot(np.exp(igae_log['2015-01-01':]), color='black', label='IGAE')

[<matplotlib.lines.Line2D object at 0x000001D6A6C61550>]

plt.plot(np.exp(fc_final.predicted_mean), color='tab:red', label='Forecast')

[<matplotlib.lines.Line2D object at 0x000001D6A6C61BE0>]

plt.fill_between(ci_final.index, np.exp(ci_final.iloc[:, 0]), np.exp(ci_final.iloc[:, 1]), color='tab:red', alpha=0.15, label='95% interval')

<matplotlib.collections.FillBetweenPolyCollection object at 0x000001D6AC34D810>

plt.title('IGAE: forecast for the next 24 months')

Text(0.5, 1.0, 'IGAE: forecast for the next 24 months')

plt.legend(loc='upper left')

<matplotlib.legend.Legend object at 0x000001D6A6C61940>

plt.grid(True)
plt.show()

We model the log of the IGAE, so we transform the forecast back with e^x. The exponential of the forecast of the log is the median forecast of the IGAE (the mean would be slightly higher). The interval widens with the horizon, just like the random walk fan chart of Workshop 3.

8.5 Adding an external variable: ARIMAX

Mexico’s economy is strongly linked to the US economy. The US stock market is a leading indicator: prices move with the expectations of investors about future activity. Can the S&P 500 of last month help to explain the IGAE? We download the S&P 500 from Yahoo Finance with the get_yahoo function of Part B and compute monthly averages:

sp500_daily = get_yahoo('^GSPC', start='1990-01-01', end='2026-10-01')

How SARIMAX treats the external variable

SARIMAX estimates a regression with SARIMA errors: \ln IGAE_t=\beta X_t+u_t, where the error u_t follows the SARIMA model. The differences (D=1) are applied to the whole equation. So if Y is the log IGAE and X is the log of the S&P 500, the model relates the annual growth of the IGAE with the annual % change of the S&P 500. Pass X in the same form as Y (log level), not as a growth rate.

log_sp = np.log(sp500_daily.resample('MS').mean()).asfreq('MS')   # log of the monthly average
x_lag1 = log_sp.shift(1).rename('log_sp500_lag1')                   # last month's S&P 500

x_in = x_lag1.reindex(igae_log.index)
sample = x_in.notna()
y_x = igae_log[sample]
x_x = x_in[sample]

res_x = fit_sarima(y_x, (1, 0, 1), (0, 1, 1, 12), exog=x_x)
res_nox = fit_sarima(y_x, (1, 0, 1), (0, 1, 1, 12))      # same sample, no exogenous variable
_ = diagnose(res_x, 'ARIMAX: Model 3 + lagged log S&P 500', plots=False)

===== ARIMAX: Model 3 + lagged log S&P 500 =====
                  coef  p-value
intercept       0.0009   0.0191
log_sp500_lag1  0.1036   0.0000
ar.L1           0.9176   0.0000
ma.L1          -0.3090   0.0000
ma.S.L12       -0.7824   0.0000
Ljung-Box p-value (lag 12): 0.0002 | (lag 24): 0.0000
AIC = -1815.19


beta_x = float(res_x.params['log_sp500_lag1'])
beta_x_p = float(res_x.pvalues['log_sp500_lag1'])
aic_x, aic_nox = float(res_x.aic), float(res_nox.aic)

The coefficient of the external variable is \hat\beta= 0.104 (p-value = 0). It is an elasticity between growth rates: if the annual % change of the S&P 500 last month is 1 percentage point higher, the annual growth of the IGAE this month changes by about 0.1 percentage points. Stock prices are much more volatile than economic activity, so we expect a small coefficient. The AIC is -1815.2 with the external variable and -1793.3 without it, on the same sample.

The catch for forecasting. To forecast the IGAE 24 months ahead with this model, we need the future values of the S&P 500, which we do not know. Using last month’s value gives us one month of known data; for longer horizons we must forecast the external variable too. Here we use the random walk with drift of Workshop 3:

future_idx = pd.date_range(igae_log.index[-1] + pd.offsets.MonthBegin(1), periods=24, freq='MS')
# months of the S&P 500 we must forecast (shift(1) needs the series to extend to the last forecast month)
months_needed = len(pd.date_range(log_sp.index[-1] + pd.offsets.MonthBegin(1), future_idx[-1], freq='MS'))

drift_sp = float(log_sp.loc['2000-01-01':].diff().mean())          # monthly drift of the random walk
sp_future_idx = pd.date_range(log_sp.index[-1] + pd.offsets.MonthBegin(1), periods=max(months_needed, 0), freq='MS')
sp_future = pd.Series(float(log_sp.iloc[-1]) + drift_sp * np.arange(1, len(sp_future_idx) + 1), index=sp_future_idx)
log_sp_ext = pd.concat([log_sp, sp_future])
x_future = log_sp_ext.shift(1).reindex(future_idx).to_frame('log_sp500_lag1')

fc_x = res_x.get_forecast(24, exog=x_future)

plt.figure(figsize=(12, 5))

<Figure size 1200x500 with 0 Axes>

plt.plot(np.exp(igae_log['2015-01-01':]), color='black', label='IGAE')

[<matplotlib.lines.Line2D object at 0x000001D6A6BB2CF0>]

plt.plot(np.exp(fc_final.predicted_mean), color='tab:red', label='SARIMA forecast')

[<matplotlib.lines.Line2D object at 0x000001D6A6BB27B0>]

plt.plot(np.exp(fc_x.predicted_mean), color='tab:blue', label='ARIMAX forecast (S&P 500 random-walk forecast as input)')

[<matplotlib.lines.Line2D object at 0x000001D6A6BB2660>]

plt.title('IGAE forecasts with and without the S&P 500')

Text(0.5, 1.0, 'IGAE forecasts with and without the S&P 500')

plt.legend(loc='upper left')

<matplotlib.legend.Legend object at 0x000001D6A6BB2F90>

plt.grid(True)
plt.show()

The ARIMAX forecast inherits the errors of the forecast of the external variable. An external variable helps most when it is known in advance: calendar effects, planned promotions, or leading indicators whose lag is at least as long as the forecast horizon.

9 CHALLENGE 3 (Parts C and D)

Work in your Notebook. Use the helper functions of this workshop (adf_test, fit_sarima, diagnose, sarima_grid, metrics, get_yahoo).

9.1 Challenge 3.1: Forecasting the sales of a brand

Use the monthly volume sales of brand 6 (qbrand6tons) of the salesbrands3 dataset of Part A. With monthly data, s=12.

Identification. Plot the log sales. Test the log level and the seasonal difference of the log (annual % change) with the ADF test. Decide d and D and justify. Look at the ACF/PACF of the stationary series and propose a starting model.
Calibration. Hold out the last 12 months as the test set. Calibrate the model on the training data step by step: for each model report the coefficients, the residual ACF/PACF and the Ljung-Box p-value, and explain why you add or remove each term. Then run sarima_grid and compare your model with the best AIC models. (Note: if you end up with d=1 and D=1, use trend='n'.)
Evaluation. Compare your model with the seasonal naive benchmark on the test set (RMSE and MAPE). Does it win?
Forecast and interpretation. Re-estimate with all the data and forecast the next 12 months after the last month of the dataset, with 95% intervals. Interpret every coefficient, the implied long-run mean annual growth of sales, and the half-life of a shock to growth.

9.2 Challenge 3.2: Does the exchange rate help to forecast the IGAE?

Download the USD/MXN exchange rate from Yahoo Finance (MXN=X, pesos per dollar) with get_yahoo, from 2000-01-01 to 2026-09-30, and compute monthly averages.

Add the log of the USD/MXN monthly average lagged one month as an external variable to Model 3 of the IGAE, SARIMA(1,0,1)(0,1,1)₁₂ (remember: pass it in the same form as Y, the log level). Interpret \hat\beta: what does a depreciation of the peso imply for the growth of economic activity, according to your model?
Compare the AIC with and without the external variable on the same sample. Then compare the RMSE/MAPE on the same 24-month test set used in class, both 24 months ahead and one month ahead. For the 24-month forecast you will need forecasts of the exchange rate: use a random walk with drift estimated only with data of the training period.
Conclude: is the ARIMAX worth it 1 month ahead? And 24 months ahead? Which external variable looks more useful for the IGAE, the exchange rate or the S&P 500 of the class example? Explain why the answer may differ by horizon.

9.3 Challenge 3.3: Decision log (required)

Include a table with one row for each modeling decision you made in Challenge 3.1 (transformation, d, D, each term added or removed, final model choice), with the evidence that supported it (test p-value, ACF/PACF spike, AIC difference).

What I will look for: correct code and results; decisions justified with evidence; interpretations that use your numbers; honest out-of-sample evaluation; and clear explanations in your own words.

10 References

Box, George E. P., and Gwilym M. Jenkins. 1970. Time Series Analysis: Forecasting and Control. Holden-Day.

Dickey, David A., and Wayne A. Fuller. 1979. “Distribution of the Estimators for Autoregressive Time Series with a Unit Root.” Journal of the American Statistical Association 74 (366): 427–31.

Engle, Robert F., and Clive W. J. Granger. 1987. “Co-Integration and Error Correction: Representation, Estimation, and Testing.” Econometrica 55 (2): 251–76.

Granger, Clive W. J., and Paul Newbold. 1974. “Spurious Regressions in Econometrics.” Journal of Econometrics 2 (2): 111–20.

Kwiatkowski, Denis, Peter C. B. Phillips, Peter Schmidt, and Yongcheol Shin. 1992. “Testing the Null Hypothesis of Stationarity Against the Alternative of a Unit Root.” Journal of Econometrics 54 (1–3): 159–78.

Ljung, Greta M., and George E. P. Box. 1978. “On a Measure of Lack of Fit in Time Series Models.” Biometrika 65 (2): 297–303.

Murray, Michael P. 1994. “A Drunk and Her Dog: An Illustration of Cointegration and Error Correction.” The American Statistician 48 (1): 37–39.