Forward Daily, Weekly and Monthly Return Density Estimation Using Fat-Tailed Models with Skew and Kurtosis Adjustments
Forecasting future price levels in financial markets isn’t about trying to predict a single outcome. That carries too much error.
Instead, we want to understand the distribution of possible future outcomes. This article outlines a method to do that in a parametric way.
Here, we estimate future return distributions using rolling statistical moments, i.e. mean, volatility, skewness, and kurtosis .
We then feed those into a Normal Inverse Gaussian distribution. This allows us to generate realistic, fat-tailed forward-looking price bands.
The complete Python notebook for the analysis is provided below.

1. Modeling NIG Forward Return Distributions
Financial returns are not always symmetric or move in straight lines. They spike, cluster, and skew.
The model presented here is built on such structure: the shape of recent return distributions, captured through rolling statistics.
We then project those statistical characteristics forward using a fat-tailed distribution.
1.1. Log Returns at the Base
We begin with log returns:

This ensures additive behavior over time horizons and simplifies the treatment of compounded returns.
1.2. Rolling Moment Estimation
From the rolling window, we estimate four key statistics:
- Mean μ
- Standard deviation σ
- Skewness γ
- Kurtosis κ
Formally:
These higher-order moments capture asymmetry and fat tails, which are much needed traits in real market returns.
1.3. Parameterizing the NIG Distribution
To approximate the return distribution, we use the ‘Normal Inverse Gaussian’ distribution.
Formally, the NIG PDF for a random variable x is:
- α > 0 → tail heaviness (larger = thinner tails)
- β ∈ (−α,α) → skewness (positive = right-skewed)
- δ > 0 → scale (volatility-like term)
- μ ∈ R → location (mean or median-like center)
- K1​→ modified Bessel function of the second kind (used in heavy-tailed models)
The above may seem like a complex formula, but essentially the purpose is to capture heavy tails, skewness and partially captures volatility clustering (via a rolling window) without simulation or complex fitting.
1.4. Estimating Parameters from Rolling Moments
We then use rolling log-return windows (e.g. 100 days) to estimate the key moments, i.e. mean, standard deviation, skewness, and kurtosis.
These map heuristically into the NIG parameters as follows:
These mappings are empirical approximations that link the intuitive moments of the data to the corresponding dimensions of the NIG.
The scale δ is set directly to the observed standard deviation, while the shape and skew parameters scale the observed kurtosis and skewness into NIG values.
1.5. Scaling by Forecast Horizon
For each forecast step h, we scale the parameters:
This assumes independent returns across time steps, consistent with standard stochastic modeling.
1.6. Converting to Price Probabilities
To translate the return distribution into price space, we:
1. Define a range around the current price using ATR:

2. Discretize this range into bins
3. Compute the midpoint return for each bin and convert to log-space
4. Evaluate the NIG PDF at each bin to get probabilities
2. Python Implementation
The process below is structured to match the methodology just outlined.
2.1. Parameters
We begin by specifying key inputs for the analysis.
These parameters control the modeling horizon, granularity, smoothing windows, and visualization scale.
Adjust as needed. Each parameter is exaplained as a comment in the snippet.
import numpy as np
import pandas as pd
import yfinance as yf
import matplotlib.pyplot as plt
from matplotlib.dates import AutoDateLocator, AutoDateFormatter
### PARAMETERS
TICKER = "TSLA"
START_DATE = "2020-01-01"
END_DATE = "2025-07-13"
INTERVAL = "1d" # Data granularity ("1d" = daily candles); used by yfinance
FREQUENCY = "monthly" # Time horizon unit for forward projections: "daily", "weekly", or "monthly"
NUM_HORIZONS = 5 # Number of forward steps; increasing adds longer-term projections
WINDOW_SIZE = 100 # Rolling window length; increasing smooths moment estimates, decreasing makes it more reactive
ATR_LENGTH = 14 # ATR window size; larger values reduce noise in range estimation
ATR_MULT = 3.0 # Range multiplier; higher values widen forecast bands, lower values narrow them
BIN_SIZE = 19 # Number of return bins; higher = finer granularity, lower = faster but coarser
BOX_WIDTH = 1.0 # Visual width of histogram bars; increase for wider boxes
MIN_ALPHA = 0.30 # Minimum opacity for low-probability bins; lower = more transparent tails
MAX_ALPHA = 0.90 # Maximum opacity for high-probability bins; higher = more emphasis on mode
PRICE_WINDOW_DAYS = 180 # Number of recent days shown in price chart; higher = longer. This is only for visual context
unit_map = {"daily": (1, "d"), "weekly": (5, "w"), "monthly": (21, "m")}
unit_days, unit_label = unit_map[FREQUENCY]
2.2 Data Retrieval and Log Returns
We download historical OHLC data using the yfinance library and compute daily log returns.
### LOG RETURNS AT THE BASE
df = yf.download(
TICKER, start=START_DATE, end=END_DATE,
interval=INTERVAL, auto_adjust=True, progress=False
)
if isinstance(df.columns, pd.MultiIndex):
df.columns = df.columns.get_level_values(0)
df["ret"] = np.log(df["Close"] / df["Close"].shift(1))
rets = df["ret"].dropna()
2.3 Rolling Moments Estimation
We estimate the mean, standard deviation, skewness, and kurtosis from a rolling window of historical returns.
### ROLLING MOMENT ESTIMATION
mu_hat = rets.rolling(WINDOW_SIZE).mean().iloc[-1]
sigma_hat = rets.rolling(WINDOW_SIZE).std(ddof=0).iloc[-1]
m3 = (rets - mu_hat)**3
m4 = (rets - mu_hat)**4
skewness = m3.rolling(WINDOW_SIZE).mean().iloc[-1] / sigma_hat**3
kurtosis = m4.rolling(WINDOW_SIZE).mean().iloc[-1] / sigma_hat**4
2.4 NIG PDF Distribution and Parameters
We define the NIG PDF using an analytical approximation of the modified Bessel function.
Parameters are literally mapped from empirical moments:

Newsletter