ARIMA Time Series Forecasting for Stock Prices

A Step-by-Step Guide to Predicting AAPL Stock Prices

What is ARIMA?

ARIMA (AutoRegressive Integrated Moving Average) is a popular statistical model for time series forecasting. It combines three components:

  • AR (AutoRegressive): Uses past values to predict future values
  • I (Integrated): Differencing to make non-stationary data stationary
  • MA (Moving Average): Uses past forecast errors in the prediction

ARIMA Parameters (p, d, q)

Parameter Description How to Choose
p Number of lag observations (AR terms) Use PACF plot
d Degree of differencing Use ADF test
q Size of moving average window Use ACF plot

Why Use ARIMA for Stock Prices?

  1. Stock prices exhibit trends (non-stationary behavior)
  2. ARIMA can handle trends through differencing
  3. Walk-forward validation simulates real trading scenarios
  4. Provides interpretable forecasts

Author: Vedanth Ramanathan Repository: vr-quantfolio-intro

1. Import Dependencies 

# Core data manipulation and numerical computing
import pandas as pd
import numpy as np

# Visualization libraries
import matplotlib.pyplot as plt
import plotly.graph_objects as go
from plotly.subplots import make_subplots

# Financial data
import yfinance as yf

# Statistical tests and models
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima.model import ARIMA

# Machine learning metrics
from sklearn.metrics import mean_squared_error

# Suppress warnings for cleaner output
import warnings
warnings.filterwarnings('ignore')

print("All dependencies imported successfully.")

2. Load Stock Data

We’ll download AAPL (Apple Inc.) stock data from Yahoo Finance spanning from 2010 to 2025. This gives us approximately 15 years of daily trading data. 

# Define stock parameters
stock = "AAPL"
start_date = "2010-12-20"
end_date = "2025-12-30"

# Download stock data from Yahoo Finance
# auto_adjust=False keeps original OHLC prices (not adjusted for splits/dividends)
# multi_level_index=False gives us simple column names
df = yf.download(
    stock, 
    start=start_date, 
    end=end_date, 
    auto_adjust=False, 
    multi_level_index=False
)

# Check for and remove any missing values
print(f"Missing values: {df.isnull().values.any()}")
df = df.dropna()

# Display basic info
print(f"\n{stock} Stock Data Summary:")
print(f"   Date Range: {df.index.min().date()} to {df.index.max().date()}")
print(f"   Trading Days: {len(df):,}")
print(f"   Columns: {list(df.columns)}")

# Show first few rows
df.head()

3. Exploratory Data Analysis (EDA)

Let’s visualize the stock price to understand its behavior over time. We’ll focus on the Close price, which is what we’ll be forecasting. 

# Create interactive price chart with Plotly
fig = go.Figure()

fig.add_trace(go.Scatter(
    x=df.index, 
    y=df['Close'],
    mode='lines',
    name='Close Price',
    line=dict(color='#2196F3', width=1.5)
))

fig.update_layout(
    title=f'{stock} Stock Price History',
    xaxis_title='Date',
    yaxis_title='Price (USD)',
    template='plotly_white',
    height=500,
    hovermode='x unified'
)

fig.show()

4. Train/Test Split

For time series data, we cannot use random splits because the order of observations matters. We split chronologically: - Training Set (80%): Historical data the model learns from - Test Set (20%): Recent data to evaluate predictions 

# Split data: 80% training, 20% testing
train_size = int(len(df) * 0.8)

train = df['Close'][:train_size]
test = df['Close'][train_size:]

print(f"Training Set: {len(train):,} observations")
print(f"   From: {train.index.min().date()} to {train.index.max().date()}")
print(f"\nTest Set: {len(test):,} observations")
print(f"   From: {test.index.min().date()} to {test.index.max().date()}")

# Visualize the split
fig = go.Figure()

fig.add_trace(go.Scatter(
    x=train.index, y=train.values,
    mode='lines', name='Training Data',
    line=dict(color='#4CAF50', width=1.5)
))

fig.add_trace(go.Scatter(
    x=test.index, y=test.values,
    mode='lines', name='Test Data',
    line=dict(color='#FF5722', width=1.5)
))

fig.update_layout(
    title='Train/Test Split Visualization',
    xaxis_title='Date',
    yaxis_title='Close Price (USD)',
    template='plotly_white',
    height=450
)

fig.show()

5. Testing for Stationarity

ARIMA requires stationary data, meaning: - Constant mean over time - Constant variance over time - No seasonal patterns

We use the Augmented Dickey-Fuller (ADF) Test: - H₀ (Null): The series has a unit root (non-stationary) - H₁ (Alternative): The series is stationary

If p-value < 0.05, we reject H₀ and conclude the series is stationary. 

# Function to perform ADF test and display results
def adf_test(series, title=''):
    """
    Perform Augmented Dickey-Fuller test for stationarity.
    
    Returns:
        bool: True if stationary (p < 0.05), False otherwise
    """
    result = adfuller(series.dropna())
    
    print(f"ADF Test Results {title}")
    print("-" * 40)
    print(f"   Test Statistic: {result[0]:.4f}")
    print(f"   P-Value: {result[1]:.6f}")
    print(f"   Lags Used: {result[2]}")
    print(f"   Observations: {result[3]}")
    print("\n   Critical Values:")
    for key, value in result[4].items():
        print(f"      {key}: {value:.4f}")
    
    is_stationary = result[1] < 0.05
    print(f"\n   Conclusion: {'STATIONARY' if is_stationary else 'NON-STATIONARY'}")
    
    return is_stationary

# Test original price series
print("Testing ORIGINAL Close Prices:\n")
is_original_stationary = adf_test(train, title='(Original Prices)')

6. Differencing: Making Data Stationary

Since stock prices are non-stationary (they trend upward/downward), we apply differencing:

\[X'_t = X_t - X_{t-1}\]

This transforms prices into daily changes, which fluctuate around zero and are typically stationary.

Why this works: - Original: [100, 102, 105, 103, 108] → Trending up - Differenced: [+2, +3, -2, +5] → Fluctuates around zero 

# Apply differencing to training data
# diff() computes: current_value - previous_value
train_diff = train.diff().dropna()

# Visualize before and after differencing
fig, axes = plt.subplots(2, 1, figsize=(14, 8))

# Original prices
axes[0].plot(train.index, train.values, color='#2196F3', linewidth=1)
axes[0].set_title('Original Prices (Non-Stationary)', fontsize=12, fontweight='bold')
axes[0].set_xlabel('Date')
axes[0].set_ylabel('Price (USD)')
axes[0].grid(True, alpha=0.3)

# Differenced prices
axes[1].plot(train_diff.index, train_diff.values, color='#4CAF50', linewidth=0.8)
axes[1].axhline(y=0, color='red', linestyle='--', alpha=0.7, label='Zero Line')
axes[1].set_title('Differenced Prices (Stationary)', fontsize=12, fontweight='bold')
axes[1].set_xlabel('Date')
axes[1].set_ylabel('Daily Change (USD)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Test if differenced series is stationary
print("\n" + "="*50)
print("Testing DIFFERENCED Close Prices:\n")
is_diff_stationary = adf_test(train_diff, title='(Differenced Prices)')

7. ARIMA Model Training

Now we’ll train an ARIMA model with: - p = 1: One autoregressive term (uses 1 past value) - d = 2: Two levels of differencing (ARIMA applies its own differencing) - q = 1: One moving average term

Walk-Forward Validation

Instead of training once and predicting all test data, we use walk-forward validation: 1. Train on all available historical data 2. Predict the next day 3. Add that actual value to training set 4. Retrain and repeat

This mimics real-world trading where you retrain models as new data arrives. 

# Prepare data for walk-forward validation
# history: growing list of training observations (starts with full train set)
history = train.tolist()

# predictions: will store our model's forecasts for each test observation
predictions = []

# ARIMA parameters
order = (1, 2, 1)  # (p, d, q)

print(f"Starting Walk-Forward ARIMA Training")
print(f"   Order: {order}")
print(f"   Test observations to predict: {len(test)}")
print("-" * 50)

# Walk-forward validation loop
# For each day in test set: train model, predict, add actual, repeat
for t in range(len(test)):
    # Train ARIMA model on all historical data available so far
    model = ARIMA(history, order=order)
    model_fit = model.fit()
    
    # Forecast the next day's price
    # forecast() returns an array, we take the first (only) element
    yhat = model_fit.forecast()[0]
    predictions.append(yhat)
    
    # Add the ACTUAL observed value to history for next iteration
    # This is key: we're NOT adding our prediction, but the real value
    actual = test.iloc[t]
    history.append(actual)
    
    # Progress update every 100 iterations
    if (t + 1) % 100 == 0 or t == 0:
        print(f"   Completed {t + 1}/{len(test)} predictions")

print(f"\nTraining complete. Generated {len(predictions)} predictions.")

8. Understanding the Cumulative Sum Reversal

⚠️ Critical Concept: ARIMA with d=2 works on twice-differenced data. Its predictions are in “change-of-changes” space, NOT actual prices!

The Problem

  • Raw ARIMA predictions: Changes in changes (meaningless for trading)
  • What we need: Actual price predictions

The Solution: Cumulative Sum (cumsum)

Since differencing = subtraction, the inverse is cumulative sum (adding up).

Example:

Original prices:     [100, 102, 105, 103, 108]
1st diff (changes):  [+2, +3, -2, +5]
2nd diff (Δchanges): [+1, -5, +7]

To reverse:
cumsum([+1, -5, +7]) → [+1, -4, +3] ≈ 1st diff
cumsum([+2, +1, -4, +3]) → [2, 3, -1, 2] ... needs anchor!

We anchor the cumsum at the last known training price to get actual predictions. 

# ============================================================
# STEP 1: Apply differencing to test data (for comparison)
# ============================================================
# We difference test data to get "changes" which ARIMA predicts
test_diff = test.diff().dropna()

# ============================================================
# STEP 2: Convert predictions to pandas Series with proper index
# ============================================================
# predictions list -> Series aligned with test dates
predictions_series = pd.Series(predictions, index=test.index)

# Difference the predictions too (since ARIMA predicts d=2 space)
predictions_diff = predictions_series.diff().dropna()

# ============================================================
# STEP 3: Cumulative sum reversal - THE KEY STEP
# ============================================================
# cumsum() reverses differencing by accumulating the changes
# But we need an anchor point - the last known price before test period

# Get the anchor: last training price
anchor_price = train.iloc[-1]
print(f"Anchor price (last training value): ${anchor_price:.2f}")

# Reverse the differencing for ACTUAL test values
# cumsum() accumulates: [a, b, c] -> [a, a+b, a+b+c]
# Adding anchor shifts the whole series to correct price level
reverse_test_diff = test_diff.cumsum() + anchor_price

# Reverse the differencing for PREDICTED values
reverse_predictions = predictions_diff.cumsum() + anchor_price

print(f"\nReversed {len(reverse_predictions)} predictions back to price space.")

9. Error Calculation

We evaluate model performance using two metrics:

Mean Squared Error (MSE)

\[MSE = \frac{1}{n}\sum_{i=1}^{n}(y_i - \hat{y}_i)^2\] - Penalizes large errors heavily (squared) - Lower is better

Symmetric Mean Absolute Percentage Error (SMAPE)

\[SMAPE = \frac{100\%}{n}\sum_{i=1}^{n}\frac{|y_i - \hat{y}_i|}{(|y_i| + |\hat{y}_i|)/2}\] - Scale-independent (percentage) - Ranges from 0% (perfect) to 200% (terrible) - Symmetric: treats over/under-prediction equally 

# ============================================================
# FINAL STEP: Calculate prediction error metrics
# ============================================================
from sklearn.metrics import mean_absolute_error, mean_squared_error

# Align the series (they should match but let's be safe)
# Use intersection of indices
common_index = reverse_test_diff.index.intersection(reverse_predictions.index)
actual = reverse_test_diff.loc[common_index]
predicted = reverse_predictions.loc[common_index]

# Calculate metrics
mae = mean_absolute_error(actual, predicted)
rmse = np.sqrt(mean_squared_error(actual, predicted))
mape = np.mean(np.abs((actual - predicted) / actual)) * 100

print("Prediction Error Metrics")
print("=" * 40)
print(f"Mean Absolute Error (MAE):     ${mae:.2f}")
print(f"Root Mean Squared Error (RMSE): ${rmse:.2f}")
print(f"Mean Absolute % Error (MAPE):  {mape:.2f}%")

print("\nInterpretation")
print("-" * 40)
print(f"On average, predictions were off by ${mae:.2f}")
print(f"This represents a {mape:.1f}% error rate")

if mape < 5:
    print("Result: Excellent accuracy for stock prediction.")
elif mape < 10:
    print("Result: Good accuracy, usable for trend analysis.")
else:
    print("Result: Moderate accuracy - use with caution.")

10. Final Visualization: Predictions vs Actual

Let’s visualize how well our ARIMA model performed:

# Create comprehensive visualization
fig = make_subplots(
    rows=2, cols=2,
    subplot_titles=(
        'Full Price History with Predictions',
        'Test Period: Actual vs Predicted',
        'Prediction Errors Over Time',
        'Error Distribution'
    ),
    vertical_spacing=0.12,
    horizontal_spacing=0.1
)

# ---- Plot 1: Full history ----
fig.add_trace(
    go.Scatter(x=train.index, y=train.values, name='Training Data',
               line=dict(color='#2196F3', width=1)),
    row=1, col=1
)
fig.add_trace(
    go.Scatter(x=actual_prices.index, y=actual_prices.values, name='Actual (Test)',
               line=dict(color='#4CAF50', width=1.5)),
    row=1, col=1
)
fig.add_trace(
    go.Scatter(x=predicted_prices.index, y=predicted_prices.values, name='Predicted',
               line=dict(color='#FF5722', width=1.5, dash='dash')),
    row=1, col=1
)

# ---- Plot 2: Zoomed test period ----
fig.add_trace(
    go.Scatter(x=actual_prices.index, y=actual_prices.values, name='Actual',
               line=dict(color='#4CAF50', width=2), showlegend=False),
    row=1, col=2
)
fig.add_trace(
    go.Scatter(x=predicted_prices.index, y=predicted_prices.values, name='Predicted',
               line=dict(color='#FF5722', width=2, dash='dash'), showlegend=False),
    row=1, col=2
)

# ---- Plot 3: Errors over time ----
errors = actual_prices - predicted_prices
fig.add_trace(
    go.Scatter(x=errors.index, y=errors.values, name='Error',
               line=dict(color='#9C27B0', width=1), showlegend=False),
    row=2, col=1
)
fig.add_hline(y=0, line_dash="dash", line_color="red", row=2, col=1)

# ---- Plot 4: Error distribution ----
fig.add_trace(
    go.Histogram(x=errors.values, nbinsx=50, name='Error Distribution',
                 marker_color='#9C27B0', showlegend=False),
    row=2, col=2
)

# Update layout
fig.update_layout(
    height=800,
    title_text=f'ARIMA({order[0]},{order[1]},{order[2]}) Model Performance on {stock}',
    template='plotly_white',
    showlegend=True,
    legend=dict(x=0.02, y=0.98)
)

fig.update_xaxes(title_text='Date', row=1, col=1)
fig.update_xaxes(title_text='Date', row=1, col=2)
fig.update_xaxes(title_text='Date', row=2, col=1)
fig.update_xaxes(title_text='Error ($)', row=2, col=2)

fig.update_yaxes(title_text='Price ($)', row=1, col=1)
fig.update_yaxes(title_text='Price ($)', row=1, col=2)
fig.update_yaxes(title_text='Error ($)', row=2, col=1)
fig.update_yaxes(title_text='Count', row=2, col=2)

fig.show()

11. Key Takeaways

What We Learned:

  1. Stationarity is Essential: ARIMA requires stationary data. Stock prices are non-stationary, so we use differencing.

  2. The ARIMA Process:

    • p (AR): How many past values influence the current value
    • d (I): How many times to difference the data
    • q (MA): How many past forecast errors influence the current forecast

  3. Walk-Forward Validation: More realistic than train-once-predict-all because it mimics real trading conditions.

  4. Cumulative Sum Reversal: Essential to convert ARIMA’s difference-space predictions back to actual prices.

  5. Model Limitations:

    • ARIMA assumes linear relationships
    • Can’t capture complex patterns or external factors
    • Stock prices are notoriously hard to predict!

Next Steps:

  • Try different (p, d, q) parameters
  • Use auto_arima from pmdarima for automatic parameter selection
  • Explore SARIMA for seasonal data
  • Compare with machine learning approaches (LSTM, Prophet)
Back to top