Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Time-Series Models for Financial Data › Implementation: Analyzing S&P 500 Returns

Michael BrenndoerferDecember 9, 202512 min read

Quantitative FinanceAll chapters

Part 8 of 15Econometrics and Statistical Modeling

Chapter 41Time-Series Models for Financial Data › Implementation: Analyzing S&P 500 Returns

Reading section 9 of 13 · Implementation: Analyzing S&P 500 Returns · 12 minSection 9 of 13 · 12 min

Section 9 of 13 · 12 min read

Time-Series Models for Financial Data › Implementation: Analyzing S&P 500 ReturnsOne section · Sections in this chapter9 / 13 ↓

Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Let's apply these concepts to real financial data, analyzing daily S&P 500 returns.

In[8]:

Code

import numpy as np

np.random.seed(42)

Simulating Financial DataLink Copied

We'll work with a simulated dataset that mimics S&P 500 characteristics, then show how to identify and fit ARIMA models.

In[9]:

Code

import numpy as np
import pandas as pd

# Simulate a price series with realistic characteristics
n_days = 1000
mu_daily = 0.0003  # ~7.5% annual return
sigma_daily = 0.012  # ~19% annual volatility

# Generate returns with a modest AR(1) structure for an identifiable example
returns = np.zeros(n_days)
phi = 0.25
for t in range(1, n_days):
    returns[t] = (
        mu_daily
        + phi * (returns[t - 1] - mu_daily)
        + sigma_daily * np.random.randn()
    )

# Generate price series
prices = 100 * np.exp(np.cumsum(returns))

# Create date index
dates = pd.date_range(start="2020-01-01", periods=n_days, freq="B")
price_series = pd.Series(prices, index=dates, name="Price")
return_series = pd.Series(returns, index=dates, name="Returns")

Visualizing Prices vs. ReturnsLink Copied

Out[10]:

Visualization

Line chart of simulated asset price rising from roughly 100 with an upward drift over 1000 days. — Simulated asset price series over 1000 trading days. The upward drift and non-constant level illustrate why prices are non-stationary and unsuitable for direct ARMA modeling.

Line chart of simulated daily returns oscillating around zero with the mean indicated by a dashed reference line. — Simulated daily return series over 1000 trading days. Returns fluctuate around a near-zero constant mean with no discernible trend, satisfying the stationarity requirement for ARMA modeling.

The price series shows a clear upward trend, violating the constant mean requirement for stationarity. Returns, however, fluctuate around a stable mean near zero, making them suitable for ARMA modeling.

Computing ACF and PACFLink Copied

In[11]:

Code

def compute_acf(series, max_lag=20):
    """Compute sample autocorrelation function."""
    n = len(series)
    mean = series.mean()
    var = series.var()
    acf = np.zeros(max_lag + 1)

    for k in range(max_lag + 1):
        cov = (
            np.sum((series[k:] - mean) * (series[:-k] - mean)) / n
            if k > 0
            else var
        )
        acf[k] = cov / var

    return acf


def compute_pacf(series, max_lag=20):
    """Compute sample partial autocorrelation using Durbin-Levinson."""
    acf = compute_acf(series, max_lag)
    pacf = np.zeros(max_lag + 1)
    pacf[0] = 1.0
    pacf[1] = acf[1]

    phi = np.zeros((max_lag + 1, max_lag + 1))
    phi[1, 1] = acf[1]

    for k in range(2, max_lag + 1):
        num = acf[k] - sum(phi[k - 1, j] * acf[k - j] for j in range(1, k))
        den = 1 - sum(phi[k - 1, j] * acf[j] for j in range(1, k))
        phi[k, k] = num / den
        for j in range(1, k):
            phi[k, j] = phi[k - 1, j] - phi[k, k] * phi[k - 1, k - j]
        pacf[k] = phi[k, k]

    return pacf

In[12]:

Code

max_lag = 20
acf_returns = compute_acf(returns, max_lag)
pacf_returns = compute_pacf(returns, max_lag)

# Confidence bounds (approximate 95% under null of white noise)
conf_bound = 1.96 / np.sqrt(n_days)

Out[13]:

Visualization

Bar chart of sample ACF for simulated returns with dashed confidence bounds, showing modest but significant autocorrelation at lag 1. — Sample ACF of simulated returns up to lag 20. The gradual decay from a significant lag-1 value is consistent with an autoregressive process rather than a moving average cutoff.

Bar chart of sample PACF for simulated returns with dashed confidence bounds, showing a prominent spike at lag 1 and insignificant values at higher lags. — Sample PACF of simulated returns up to lag 20. The sharp spike at lag 1 followed by values within the dashed 95% confidence bounds identifies the process as AR(1).

The ACF shows the lag-1 autocorrelation is significant, with subsequent lags decaying. The PACF shows a significant spike at lag 1 followed by values within the confidence bands. This pattern, with PACF cutting off after lag 1 and ACF decaying, strongly suggests an AR(1) model.

Estimating an AR(1) ModelLink Copied

We'll estimate the AR(1) parameters using ordinary least squares, which is consistent for stationary AR models.

In[14]:

Code

# Prepare data for AR(1) estimation
Y = returns[1:]  # Y_t
Y_lag = returns[:-1]  # Y_{t-1}

# OLS estimation: Y_t = c + phi * Y_{t-1} + epsilon_t
X = np.column_stack([np.ones(len(Y_lag)), Y_lag])
beta = np.linalg.lstsq(X, Y, rcond=None)[0]
c_hat, phi_hat = beta

# Compute residuals
Y_hat = X @ beta
residuals = Y - Y_hat
sigma_hat = np.std(residuals)

# Compute standard errors
XtX_inv = np.linalg.inv(X.T @ X)
var_beta = sigma_hat**2 * XtX_inv
se_c = np.sqrt(var_beta[0, 0])
se_phi = np.sqrt(var_beta[1, 1])

Out[15]:

Console

AR(1) Model Estimation Results
========================================
Constant (c):     0.000454  (SE: 0.000372)
AR(1) coef (φ):   0.2433  (SE: 0.0307)
Residual std:     0.011748

Implied mean:     0.000600
Stationarity:     Yes (|φ| = 0.2433)

The estimated AR(1) coefficient is close to our simulated value of 0.25. The model is stationary since $|\hat{\phi}| < 1$ .

Diagnostic CheckingLink Copied

In[16]:

Code

from scipy import stats


# Ljung-Box test for residual autocorrelation
def ljung_box_test(residuals, max_lag=10):
    """Compute Ljung-Box Q statistic."""
    n = len(residuals)
    acf_resid = compute_acf(residuals, max_lag)

    Q = (
        n
        * (n + 2)
        * sum(acf_resid[k] ** 2 / (n - k) for k in range(1, max_lag + 1))
    )
    p_value = 1 - stats.chi2.cdf(Q, df=max_lag)

    return Q, p_value


Q_stat, p_val = ljung_box_test(residuals)
acf_residuals = compute_acf(residuals, max_lag)

Out[17]:

Console

Residual Diagnostics
========================================
Ljung-Box Q(10):  7.45
p-value:          0.6828

Residual mean:    0.000000
Residual std:     0.011748
Skewness:         0.1188
Kurtosis:         0.0656

The Ljung-Box test checks whether residual autocorrelations are collectively significant. A high p-value suggests the residuals behave like white noise, indicating an adequate model.

Out[18]:

Visualization

Time series plot of AR(1) model residuals oscillating around zero with no obvious trend or pattern. — AR(1) residuals over the 1000-day sample period. The absence of systematic patterns or clustering around zero supports the white noise assumption for the model innovations.

Bar chart of residual ACF with all bars within the dashed 95% confidence bounds, indicating white noise residuals. — Sample ACF of AR(1) residuals up to lag 20. All bars fall within the dashed 95% confidence bounds, confirming that no significant serial correlation remains after fitting the AR(1) model.

ForecastingLink Copied

With the estimated model, we can generate forecasts with confidence intervals.

In[19]:

Code

# Generate forecasts
forecast_horizon = 20
Y_T = returns[-1]  # Last observed return
mu_unconditional = c_hat / (1 - phi_hat)

forecasts = np.zeros(forecast_horizon)
forecast_var = np.zeros(forecast_horizon)

for h in range(forecast_horizon):
    # Point forecast
    forecasts[h] = mu_unconditional + (phi_hat ** (h + 1)) * (
        Y_T - mu_unconditional
    )

    # Forecast variance
    forecast_var[h] = (
        sigma_hat**2 * (1 - phi_hat ** (2 * (h + 1))) / (1 - phi_hat**2)
    )

# 95% confidence intervals
lower_95 = forecasts - 1.96 * np.sqrt(forecast_var)
upper_95 = forecasts + 1.96 * np.sqrt(forecast_var)

# 80% confidence intervals
lower_80 = forecasts - 1.28 * np.sqrt(forecast_var)
upper_80 = forecasts + 1.28 * np.sqrt(forecast_var)

Out[20]:

Visualization

Line plot showing point forecasts converging to mean with widening confidence bands. — Twenty-day return forecasts from the AR(1) model with 80% and 95% confidence intervals. The point forecast rapidly converges to the dashed unconditional mean, while the nested confidence intervals widen to the unconditional variance limit.

The forecasts converge quickly to the unconditional mean because the AR(1) coefficient is small. The confidence intervals widen initially but stabilize as they approach the unconditional variance.

Using statsmodels for ARIMALink Copied

In practice, we use established libraries for ARIMA estimation. The statsmodels package provides the ARIMA tools used below.

In[21]:

Code

from statsmodels.tsa.stattools import adfuller

# Augmented Dickey-Fuller test for stationarity
adf_result = adfuller(returns)

Out[22]:

Console

Augmented Dickey-Fuller Test for Returns
=============================================
ADF Statistic:    -24.6312
p-value:          0.000000
Critical values:
  1%:         -3.4369
  5%:         -2.8644
  10%:         -2.5683

Conclusion: Stationary (reject unit root)

The ADF test strongly rejects the null hypothesis of a unit root, confirming that returns are stationary and suitable for ARMA modeling without differencing.

In[23]:

Code

# Fit ARIMA models with different specifications
from statsmodels.tsa.arima.model import ARIMA

models = {}
for p in range(3):
    for q in range(3):
        if p == 0 and q == 0:
            continue  # Skip trivial model
        try:
            model = ARIMA(returns, order=(p, 0, q))
            fitted = model.fit()
            models[(p, 0, q)] = {
                "aic": fitted.aic,
                "bic": fitted.bic,
                "model": fitted,
            }
        except:
            pass  # Skip failed fits

# Sort models by BIC for comparison
sorted_models = sorted(models.items(), key=lambda x: x[1]["bic"])

Out[24]:

Console

Model Comparison (ARIMA)
==================================================
Model           AIC             BIC            
--------------------------------------------------
ARIMA(1,0,0)      -6045.13        -6030.41       
ARIMA(0,0,1)      -6041.88        -6027.15       
ARIMA(1,0,1)      -6043.13        -6023.50       
ARIMA(2,0,0)      -6043.13        -6023.50       
ARIMA(0,0,2)      -6042.45        -6022.82       
ARIMA(2,0,1)      -6041.13        -6016.59       
ARIMA(1,0,2)      -6039.96        -6015.42       
ARIMA(2,0,2)      -6039.59        -6010.14

Information criteria help select the best model balancing fit and parsimony. The BIC, with its heavier penalty for complexity, often selects simpler models appropriate for financial data.

In[25]:

Code

# Fit best model (AR(1) based on BIC)
if sorted_models:
    best_order = sorted_models[0][0]
    best_model = models[best_order]["model"]
else:
    # Fallback to AR(1) if no models were successfully fit
    from statsmodels.tsa.arima.model import ARIMA

    model = ARIMA(returns, order=(1, 0, 0))
    best_model = model.fit()
    best_order = (1, 0, 0)

Out[26]:

Console


Best Model: ARIMA(1, 0, 0)
==================================================
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const          0.0006      0.000      1.217      0.223      -0.000       0.002
ar.L1          0.2433      0.032      7.723      0.000       0.182       0.305
sigma2         0.0001   6.08e-06     22.679      0.000       0.000       0.000
==============================================================================

The model summary confirms that the AR(1) coefficient is statistically significant and close to the true value of 0.25. The low p-values indicate strong evidence against zero coefficients.

ARIMA for Price SeriesLink Copied

When working with prices directly, we need ARIMA with $d=1$ . This models returns (first differences) using ARMA.

In[27]:

Code

# Fit ARIMA(1,1,0) to prices, which is equivalent to AR(1) on returns
price_model = ARIMA(prices, order=(1, 1, 0))
price_fitted = price_model.fit()

# Generate price forecasts
price_forecast = price_fitted.get_forecast(steps=20)
price_pred = price_forecast.predicted_mean
price_ci = price_forecast.conf_int(alpha=0.05)

Out[28]:

Visualization

Price forecast plot with confidence bands widening over forecast horizon. — Price forecasts from an ARIMA(1,1,0) model over a 20-day horizon. The shaded 95% confidence interval widens significantly faster than for returns, illustrating the accumulation of uncertainty in integrated I(1) processes.

The confidence interval for price forecasts widens dramatically compared to return forecasts. This reflects the accumulation of return forecast errors: uncertainty compounds over time when predicting integrated series.