Multi-Step Forecasting With Fourier Transform in Python
How to Use Python to Develop a Fourier Time Series Model
If you’ve ever worked with time series data, you’ve probably heard some version of this line:
“Just take the Fourier transform.”
It’s often said with the same casual confidence as “just normalize the data” or “just smooth it.” And yet, for many people, the Fourier transform remains strangely abstract. They know it turns time into frequency. They know it’s useful. But they don’t really know what it does.
Stop Guessing. Start Trading with The Signal Beyond report 🚀
COT Data • Technical charts • Equity Arbitrage • Seasonality
A Layman’ Terms Introduction to Fourier
A time series is simple at heart. It’s just a list of values indexed by time.
The temperature every hour
The voltage measured from a sensor
The closing price of a stock each day
The air pressure inside a machine as it runs
You can draw it as a line: time on the x-axis, value on the y-axis.
In this view, the signal feels concrete. You can point to spikes, trends, noise, and sudden changes. You can say “something happened here” and mean a specific moment. But this view also hides structure.
Some patterns are easy to see in time. Others are almost invisible.
When you look at a time series, you are implicitly asking:
“What is the value at each moment?”
The Fourier transform asks a different question:
“What oscillations make this signal up?”
Instead of asking when things happen, it asks how often things repeat. That’s the key mental shift.
You are not transforming the signal into something else. You are re-describing the same signal in a different language.
The Fourier transform assumes something radical but surprisingly effective:
Any time series can be written as a sum of pure waves.
Not approximately. Exactly, under broad conditions. Each of those waves has:
A frequency (how fast it oscillates)
An amplitude (how strong it is)
A phase (where it starts)
The transform figures out which waves you need, and how much of each.
Low-frequency waves capture slow trends.
High-frequency waves capture rapid changes.
Everything in between fills in the detail.
After the Fourier transform, you no longer have a value for every time step.
Instead, you have a spectrum.
For each frequency, the spectrum tells you how much that frequency contributes to the signal. This answers questions like:
Is there a dominant cycle?
Are there hidden periodic components?
Is the signal mostly smooth or mostly jagged?
Where does the noise live?
Importantly, this information can be hard or impossible to see in the time domain.
A weak but consistent oscillation can be buried in noise over time, yet appear as a clean spike in the frequency domain.
When you apply the standard Fourier transform, you lose direct information about when things happen.
The spectrum doesn’t say:
When a frequency starts
When it stops
Whether it was present only briefly
It only says:
“Across the whole signal, this frequency exists this much.”
This is why Fourier analysis works best when the signal is roughly stationary, meaning its statistical properties don’t change much over time.
If the signal evolves, the transform smears everything together.
This tradeoff is not a flaw. It’s a consequence of the question you chose to ask.
One of the most practical things the Fourier transform lets you do is filtering.
The idea is simple:
Transform the signal to frequency space
Reduce or remove certain frequencies
Transform back to time
Low-pass filters remove rapid fluctuations. High-pass filters remove slow trends. Band-pass filters isolate specific ranges.
In the time domain, these operations look mysterious. In the frequency domain, they’re often just multiplication.
This is one of the deepest insights Fourier gives you:
Convolution in time becomes multiplication in frequency.
That single fact powers enormous chunks of signal processing.
Noise and Structure
Many kinds of noise are high-frequency by nature. They change rapidly, randomly, and without structure.
Many meaningful signals change more slowly. The Fourier transform often separates these cleanly. Not always. Not magically. But often enough to matter.
This is why denoising, smoothing, and compression so often start in frequency space.
Another useful mental model is projection. Imagine your time series as a vector in a very high-dimensional space.
The sine and cosine waves form a basis for that space. The Fourier transform computes how much of your signal points in each of those basis directions.
You’re not destroying the signal. You’re expressing it in a different coordinate system.
Inverse Transform
The Fourier transform alone is interesting. The inverse Fourier transform is what makes it powerful. It guarantees that, unless you deliberately throw information away, you can come back.
Time → frequency → time.
Same signal. Different view. This reversibility is what lets you edit signals safely in frequency space, then return them to the real world.
At its core, the Fourier transform does not make your data smarter.
It makes you ask a smarter question. Instead of “what happened when,” you ask:
“What patterns exist?”
“What rhythms are present?”
“What structure repeats?”
Sometimes that’s exactly the question you needed all along.
And sometimes, seeing your data this way changes how you understand the system behind it.
That, more than any formula, is why the Fourier transform endures.
Forecasting Multiple Steps Into the Future Using Fourier
The code seen in this section builds a time-series forecasting model that separates a price series into two parts:
A simple, extrapolatable trend
A repeating (oscillatory) residual pattern
It then forecasts each part forward in time and recombines them.
The key idea is:
Trends should be modeled with simple models that extrapolate cleanly, and periodic structure should be modeled with Fourier components, but only on what’s left after the trend is removed.
This avoids many common mistakes people make when applying Fourier transforms directly to raw price data.
The model first converts prices to log prices.
Why?
Log prices behave more like additive processes
Linear trends in log space correspond to exponential growth in price space
Residuals tend to be more stable and symmetric
From this point on, everything happens in log space until the final plot.
Step 1: Estimate the trend using rolling linear regression
Instead of fitting a global trend over all history, the code:
Takes the most recent window of data
Fits a straight line to it
Extracts:
The current level (where the trend is now)
The slope (how fast it’s moving)
This is important because:
Linear trends extrapolate cleanly
No parameters are re-optimized during forecasting
The trend is deliberately kept simple and boring
Think of this as the model’s answer to:
“If nothing fancy happens, where is this series already headed?”
Step 2: Remove the trend to get residuals
Next, the code subtracts the estimated trend from recent log prices.
What remains is the residual:
Cycles
Mean reversion
Oscillations
Short-term structure
This step is crucial. Fourier methods assume stationarity. Raw prices are not stationary. Residuals after detrending are much closer.
So instead of asking:
“What frequencies exist in the price?”
It asks:
“What frequencies exist in the deviations from trend?”
That’s a much safer question.
Step 3: Fourier analysis on the residuals only
Now the Fourier transform comes in. What the code does conceptually:
Takes a fixed window of recent residuals
Applies a Hann window to reduce edge artifacts
Computes the Fourier transform
Keeps only the top k strongest frequencies
Stores each frequency’s:
Frequency
Amplitude
Phase
Important methodological choices here:
No refitting during forecasting
Fixed number of components
No attempt to “optimize” frequencies
No use of the DC (mean) component
This avoids overfitting and forces the model to capture only the strongest, most persistent cycles.
Step 4: Project those cycles forward
Each retained Fourier component is just a cosine wave.
The model:
Extends those cosine waves forward in time
Keeps their frequency, amplitude, and phase fixed
Sums them together to form a residual forecast
This is equivalent to saying:
“If these cycles have existed recently, assume they continue unchanged.”
That assumption may be wrong, but it’s explicit and testable.
Step 5: Extrapolate the trend forward
Separately, the trend is extended forward using the same slope estimated earlier.
No re-estimation. No adjustment. Just a straight continuation.
This gives a trend forecast.
Step 6: Combine trend + residual forecast
The final forecast is simply:
forecast = trend_forecast + residual_forecastStill in log space.
Only at plotting time does the code convert back to prices using exp.
Step 7: Compare against dumb baselines
This is a critical part of the design. The code always benchmarks against:
Persistence: tomorrow equals today
Slope continuation: last price + slope × time
It then computes RMSE ratios like:
Fourier vs persistence
Fourier vs slope
This answers the only question that matters:
“Is this more complex model actually better than trivial alternatives?”
If it doesn’t beat them, it fails honestly.
The following is the resulting graph from the forecast on a specific asset.
Use the following code to implement the experiment.
import numpy as np
from dataclasses import dataclass
@dataclass
class ForecastResult:
forecast: np.ndarray # Final forecast (log price)
trend_forecast: np.ndarray # Trend component
residual_forecast: np.ndarray # Fourier residual component
persistence: np.ndarray # Baseline: last value repeated
slope_continuation: np.ndarray # Baseline: linear extrapolation
def rolling_linear_regression(y: np.ndarray, window: int) -> tuple[float, float]:
"""
Fit linear regression on the last `window` points.
Returns (intercept, slope) where intercept is at t=window-1.
"""
n = min(len(y), window)
y_win = y[-n:]
t = np.arange(n)
t_mean = t.mean()
y_mean = y_win.mean()
slope = np.sum((t - t_mean) * (y_win - y_mean)) / np.sum((t - t_mean) ** 2)
intercept = y_mean - slope * t_mean
# Return level at end of window and slope
level = intercept + slope * (n - 1)
return level, slope
def extract_top_k_fourier(residual: np.ndarray, k: int) -> list[tuple[float, float, float]]:
"""
Apply Hann window and extract top k Fourier components.
Returns list of (frequency_idx, amplitude, phase) tuples.
"""
n = len(residual)
# Apply Hann window to reduce spectral leakage
window = np.hanning(n)
windowed = residual * window
# FFT
fft = np.fft.rfft(windowed)
freqs = np.fft.rfftfreq(n)
# Get magnitudes (skip DC component at index 0)
magnitudes = np.abs(fft[1:])
# Find top k components by magnitude
top_indices = np.argsort(magnitudes)[-k:][::-1]
components = []
for idx in top_indices:
actual_idx = idx + 1 # Offset for skipped DC
amp = 2 * np.abs(fft[actual_idx]) / (n * np.mean(window)) # Correct for window
phase = np.angle(fft[actual_idx])
freq = freqs[actual_idx]
components.append((freq, amp, phase))
return components
def project_fourier_components(
components: list[tuple[float, float, float]],
n_history: int,
n_forecast: int
) -> np.ndarray:
"""
Extend sinusoids forward from end of history window.
"""
t_forecast = np.arange(n_history, n_history + n_forecast)
projection = np.zeros(n_forecast)
for freq, amp, phase in components:
projection += amp * np.cos(2 * np.pi * freq * t_forecast + phase)
return projection
def forecast_fourier(
prices: np.ndarray,
n_steps: int,
trend_window: int = 60,
fourier_window: int = 120,
k_components: int = 5
) -> ForecastResult:
"""
Forecast n_steps ahead using trend + Fourier residual decomposition.
Parameters:
-----------
prices : array
Raw price series
n_steps : int
Number of steps to forecast
trend_window : int
Window for rolling linear regression (default 60)
fourier_window : int
Window for Fourier analysis (default 120)
k_components : int
Number of Fourier components to keep (default 5, typical range 3-10)
Returns:
--------
ForecastResult with forecast and baselines
"""
# Work in log space
log_prices = np.log(prices)
# Step 1: Fit trend via rolling linear regression
level, slope = rolling_linear_regression(log_prices, trend_window)
# Step 2: Compute residual = log_price - trend
n = len(log_prices)
trend_window_actual = min(n, trend_window)
t_trend = np.arange(trend_window_actual)
trend_in_window = (level - slope * (trend_window_actual - 1)) + slope * t_trend
# Use fourier_window for residual analysis
fourier_n = min(n, fourier_window)
log_recent = log_prices[-fourier_n:]
# Reconstruct trend over fourier window
trend_full = level + slope * (np.arange(fourier_n) - (fourier_n - 1))
residual = log_recent - trend_full
# Step 3: Fourier decomposition on residual (with Hann window)
components = extract_top_k_fourier(residual, k_components)
# Step 4: Project components forward
residual_forecast = project_fourier_components(components, fourier_n, n_steps)
# Step 5: Extrapolate trend
t_forecast = np.arange(1, n_steps + 1)
trend_forecast = level + slope * t_forecast
# Step 6: Combine
forecast = trend_forecast + residual_forecast
# Baselines
persistence = np.full(n_steps, log_prices[-1])
slope_continuation = log_prices[-1] + slope * t_forecast
return ForecastResult(
forecast=forecast,
trend_forecast=trend_forecast,
residual_forecast=residual_forecast,
persistence=persistence,
slope_continuation=slope_continuation
)
def evaluate_forecast(actual: np.ndarray, result: ForecastResult) -> dict:
"""
Compare forecast against baselines using RMSE.
"""
log_actual = np.log(actual)
def rmse(pred, true):
return np.sqrt(np.mean((pred - true) ** 2))
return {
'fourier_rmse': rmse(result.forecast, log_actual),
'persistence_rmse': rmse(result.persistence, log_actual),
'slope_rmse': rmse(result.slope_continuation, log_actual),
'fourier_vs_persistence': rmse(result.forecast, log_actual) / rmse(result.persistence, log_actual),
'fourier_vs_slope': rmse(result.forecast, log_actual) / rmse(result.slope_continuation, log_actual),
}
import yfinance as yf
import pandas as pd
def get_timeseries(ticker: str = None, ema_span= None) -> pd.DataFrame:
df = yf.download(ticker, start="2000-01-01", interval="1d")
df = df[["Close", "Low", "High", "Open"]].rename(columns={"Close": "close", "Low": "low", "High": "high", "Open": "open"})
df.columns = df.columns.get_level_values(0)
if ema_span is not None:
df['close_raw'] = df['close']
df['close'] = df['close'].ewm(span=ema_span, adjust=False).mean()
return df
symbol = "TSLA"
df = get_timeseries(ticker=symbol)
df = df['close']
prices = np.array(df)
# Demo
if __name__ == "__main__":
import matplotlib.pyplot as plt
# Split: use first 280 for fitting, last 20 for validation
n_forecast = 50
len_df = len(prices)
train_prices = prices[:len_df-n_forecast]
test_prices = prices[len_df-n_forecast:]
# Forecast
result = forecast_fourier(
train_prices,
n_steps=n_forecast,
trend_window=20,
fourier_window=50,
k_components=50
)
# Evaluate
metrics = evaluate_forecast(test_prices, result)
print("Forecast Evaluation (RMSE in log space):")
print(f" Fourier model: {metrics['fourier_rmse']:.6f}")
print(f" Persistence: {metrics['persistence_rmse']:.6f}")
print(f" Slope continuation:{metrics['slope_rmse']:.6f}")
print()
print("Relative performance (< 1.0 means Fourier wins):")
print(f" vs Persistence: {metrics['fourier_vs_persistence']:.3f}")
print(f" vs Slope: {metrics['fourier_vs_slope']:.3f}")
# Plot
fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=False)
# Top plot: Full context + forecast
ax1 = axes[0]
n_context = 60 # Show last 60 points of training data for context
t_context = np.arange(-n_context, 0)
t_forecast = np.arange(0, n_forecast)
# Training context
ax1.plot(t_context, train_prices[-n_context:], 'k-', lw=1.5, label='History')
# Actual test data
ax1.plot(t_forecast, test_prices, 'ko-', lw=2, markersize=5, label='Actual')
# Forecasts (convert from log space)
ax1.plot(t_forecast, np.exp(result.forecast), 'b-', lw=2, label='Fourier')
ax1.plot(t_forecast, np.exp(result.persistence), 'r--', lw=1.5, alpha=0.7, label='Persistence')
ax1.plot(t_forecast, np.exp(result.slope_continuation), 'g--', lw=1.5, alpha=0.7, label='Slope')
ax1.axvline(0, color='gray', linestyle=':', alpha=0.5)
ax1.set_xlabel('Steps from forecast origin')
ax1.set_ylabel('Price')
ax1.set_title('Forecast Comparison')
ax1.legend(loc='upper left')
ax1.grid(True, alpha=0.3)
# Bottom plot: Decomposition of Fourier forecast
ax2 = axes[1]
ax2.plot(t_forecast, np.exp(result.trend_forecast), 'c-', lw=2, label='Trend extrapolation')
ax2.plot(t_forecast, result.residual_forecast, 'm-', lw=2, label='Residual projection')
ax2.axhline(0, color='gray', linestyle=':', alpha=0.5)
ax2.set_xlabel('Steps from forecast origin')
ax2.set_ylabel('Value')
ax2.set_title('Forecast Decomposition (Trend in price space, Residual in log space)')
ax2.legend(loc='upper left')
ax2.grid(True, alpha=0.3)
# Add metrics annotation
metrics_text = (
f"RMSE (log): Fourier={metrics['fourier_rmse']:.4f}, "
f"Persist={metrics['persistence_rmse']:.4f}, Slope={metrics['slope_rmse']:.4f}\n"
f"Fourier vs Persist: {metrics['fourier_vs_persistence']:.3f}x, "
f"Fourier vs Slope: {metrics['fourier_vs_slope']:.3f}x"
)
fig.text(0.5, 0.02, metrics_text, ha='center', fontsize=10,
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
plt.tight_layout()
plt.subplots_adjust(bottom=0.12)You can also find the code on my GitHub. The following Figure shows another forecast.
Make sure to understand the limitations. You won’t be able to always find forecasts that look like the actual outcomes, especially with a simple model like this.
This model assumes:
Long-term movement is dominated by a smooth trend
Short-term structure is mostly oscillatory
Cycles persist longer than noise
Trend and cycles should be modeled separately
Simplicity beats over-optimization
It does not assume:
Mean reversion
Regime switching
Predictable shocks
That markets are “efficient” or “inefficient”
It just asks:
“If the recent past decomposes cleanly into trend + cycles, what happens if we extend both?”
Master Deep Learning techniques tailored for time series and market analysis🔥
My book breaks it all down from basic machine learning to complex multi-period LSTM forecasting while going through concepts such as fractional differentiation and forecasting thresholds. Get your copy here 📖!