edgekit.timeseries

Time-series analysis and causal feature engineering — the diagnostics you run before you build a strategy (is this series trending, mean-reverting, or a random walk?) and the machinery that turns raw OHLC into a leakage-free feature matrix a rule or a model can read.

What's inside. Return and volatility construction (returns, log_returns, rolling_volatility, ewma_volatility, realized_vol), the autocorrelation / persistence battery (acf, pacf, autocorr, ljung_box, variance_ratio, hurst_exponent), stationarity testing (adf_stat, stationarity_report), transforms (difference, seasonal_decompose), and the causal feature builders (make_features, add_lags).

→numpy + pandas only — no statsmodels, no scipy

Every estimator here is implemented from scratch on top of numpy and pandas: the ACF/PACF (Durbin-Levinson), the Lo-MacKinlay variance ratio, the Hurst exponent, and the ADF t-statistic. That keeps edgekit dependency-light and the arithmetic auditable, at the cost of the extras a full stats package would give you (p-value tables, lag-order selection). Compare the returned statistics to the standard critical values noted below.

Returns & volatility#

returns#

Period-over-period returns from a price series — simple by default, log when log=True. The first element is nan (no prior price to difference against). Simple returns aggregate across assets; log returns aggregate across time — pick the one that matches what you are summing.

returns(prices, log=False) -> np.ndarray   # rets[0] = nan

Param	Type	Default	Meaning
`prices`	array / Series	—	The price level series (e.g. close).
`log`	bool	False	False = simple pct-change; True = log returns.

Returns: a numpy array aligned to prices, with [0] = nan.

import edgekit as ek
rets = ek.timeseries.returns(bars.close)            # simple
lrets = ek.timeseries.returns(bars.close, log=True) # log

log_returns#

Log returns as a named shortcut for returns(prices, log=True). Additive across time and roughly symmetric, which is why most statistical machinery (vol, ACF, stationarity tests) is run on log returns rather than raw prices.

log_returns(prices) -> np.ndarray          # [0] = nan

prices — the price level series.

Returns: a numpy array of log returns, [0] = nan.

rolling_volatility#

Rolling standard deviation of returns over window bars — the trailing estimate of how big a normal move is. Pass periods_per_year to annualise (multiply by sqrt(periods_per_year)), so a daily series becomes a comparable annual vol.

rolling_volatility(rets, window=20, periods_per_year=None) -> np.ndarray

Param	Type	Default	Meaning
`rets`	array	—	The return series (e.g. from returns()).
`window`	int	20	Rolling lookback in bars.
`periods_per_year`	int \| None	None	If given, annualise by sqrt(periods_per_year).

Returns: a numpy array of rolling vol (leading positions NaN until the window fills).

vol = ek.timeseries.rolling_volatility(rets, 20, periods_per_year=252)  # annualised daily vol

ewma_volatility#

RiskMetrics exponentially-weighted volatility: variance updated as sigma2[t] = lam*sigma2[t-1] + (1-lam)*r[t]^2. Reacts to a vol shock faster than a flat rolling window because recent squared returns carry more weight. lam=0.94 is the RiskMetrics daily default.

ewma_volatility(rets, lam=0.94) -> np.ndarray

rets — the return series.
lam — decay factor in (0,1); larger = longer memory (0.94 daily, ~0.97 monthly).

Returns: a numpy array of EWMA volatility (per-bar sigma).

realized_vol#

Realized volatility over window bars — the root of summed squared returns, the backward-looking volatility actually delivered over the window. The empirical counterpart to a volatility forecast.

realized_vol(rets, window=20) -> np.ndarray

rets — the return series.
window — the lookback (default 20).

Returns: a numpy array of realized vol.

Autocorrelation & persistence#

acf#

Autocorrelation function — the correlation of the series with itself at lags 0..nlags.acf[0] is always 1. A slow decay signals persistence (trend); a sharp cut-off after a few lags signals short-memory (MA) structure. Compare each bar to the significance band ±1.96/sqrt(n) — spikes outside it are significant.

acf(x, nlags=40) -> np.ndarray             # acf[0] = 1

x — the series (usually returns, not prices).
nlags — number of lags to compute (default 40).

Returns: a numpy array of length nlags+1, [0]=1.

ac = ek.timeseries.acf(rets, nlags=40)
band = 1.96 / len(rets)**0.5   # significance band

pacf#

Partial autocorrelation function via Durbin-Levinson — the correlation at lag k after removing the effect of all shorter lags. Where the ACF gauges MA order, the PACF gauges AR order: a PACF that cuts off after lag p points to an AR(p) process.

pacf(x, nlags=40) -> np.ndarray            # Durbin-Levinson

x — the series.
nlags — number of lags (default 40).

Returns: a numpy array of partial autocorrelations, [0]=1.

autocorr#

Scalar autocorrelation at a single lag — the convenience form when you only need one number (e.g. lag-1 autocorrelation of returns, the sign of which hints trend vs reversal). It is the value stationarity_report stores under autocorr_lag1.

autocorr(x, lag=1) -> float

x — the series.
lag — the lag to evaluate (default 1).

Returns: a float autocorrelation coefficient.

ljung_box#

Ljung-Box portmanteau test for joint autocorrelation across the first lagslags — one number that answers “is there anylinear dependence here, or is this white noise?”. Run it on returns (are they predictable?) and on squared returns (is there volatility clustering?).

ljung_box(x, lags=10) -> dict

x — the series.
lags — number of lags to pool (default 10).

Returns: a dict with keys "stat" (the Q statistic) and "dof" (degrees of freedom = lags). Compare stat to a chi-squared distribution with dof degrees of freedom — a large stat rejects white noise.

lb = ek.timeseries.ljung_box(rets, lags=10)
# lb == {"stat": ..., "dof": 10}; compare stat to chi2(dof)

variance_ratio#

Lo-MacKinlay variance ratio on log prices: the variance of k-period returns divided by k× the variance of 1-period returns. Under a random walk this is 1; positive autocorrelation (trend) inflates it above 1, mean-reversion pulls it below 1.

variance_ratio(prices, k=2) -> float       # >1 trend, <1 revert, ~1 random walk

prices — the price series (the function takes logs internally).
k — the aggregation horizon (default 2).

Returns: a float variance ratio.

vr = ek.timeseries.variance_ratio(bars.close, k=5)   # >1 => trending on a 5-bar horizon

hurst_exponent#

The Hurst exponent from a rescaled-range / dispersion fit across lags. ≈0.5 is a random walk, >0.5 is persistent (trending, long memory), <0.5 is anti-persistent (mean-reverting). A model-free companion to the variance ratio.

hurst_exponent(x, max_lag=50) -> float     # ~0.5 walk, >0.5 persist, <0.5 revert

x — the series (prices or a spread).
max_lag — the largest lag in the scaling fit (default 50).

Returns: a float Hurst exponent, typically in (0,1).

Stationarity#

adf_stat#

The Augmented Dickey-Fuller t-statistic — tests the null of a unit root (non-stationary). A sufficiently negative statistic rejects the unit root and calls the series stationary. Rule of thumb: < -2.86 is significant at the 5% level, < -3.43 at 1%.

adf_stat(x, lags=1) -> float               # < ~-2.86 (5%) => stationary

x — the series to test.
lags — number of augmenting lagged-difference terms (default 1).

Returns: a float ADF t-statistic (more negative = more stationary).

!It is a statistic, not a p-value

adf_stat returns the raw t-stat only — there is no scipy on the dependency list to look up an exact p-value. Compare it against the standard MacKinnon critical values above.

stationarity_report#

A one-call diagnostic that runs the ADF test, the variance ratio, the Hurst exponent, and lag-1 autocorrelation together and folds them into a single verdict — the fastest way to characterise a fresh series before you decide whether to build a trend model, a reversion model, or nothing.

stationarity_report(x) -> dict

x — the series (price or return).

Returns: a dict with keys "adf_stat", "variance_ratio", "hurst", "autocorr_lag1", and "verdict" — where verdict is one of "trending", "mean-reverting", or "random-walk".

rep = ek.timeseries.stationarity_report(bars.close)
rep["verdict"]   # "trending" | "mean-reverting" | "random-walk"

Transforms#

difference#

The d-th order difference of a series — the standard way to stationarise a level series (a first difference of prices is close cousin to returns). Differencing removes a stochastic trend so downstream stats stop lying to you.

difference(x, d=1) -> np.ndarray

x — the series.
d — order of differencing (default 1).

Returns: a numpy array of the differenced series (shorter / NaN-padded by d).

seasonal_decompose#

Classic additive decomposition into trend, seasonal, and residual components over a fixed period — useful for spotting a day-of-week or intraday-session cycle in returns or volume.

seasonal_decompose(series, period) -> dict

series — a pandas Series with a meaningful index.
period — the seasonal period in bars (e.g. 5 for weekly on daily data).

Returns: a dict with keys "trend", "seasonal", and "resid" — each a pandas Series aligned to the input.

dec = ek.timeseries.seasonal_decompose(daily_ret, period=5)
dec["seasonal"]   # the repeating weekly component

Feature engineering#

✕make_features is causal — do not re-introduce look-ahead

make_features shifts every column by one bar before returning it, so row i holds only information known at the close of bar i-1. That is deliberate: the feature matrix is safe to hand straight to a rule or a model without leaking the current bar.

The corollary: do not add your own un-lagged columns next to it, and do not un-shift what it gives you. Mixing a causal matrix with a same-bar feature is the classic way look-ahead sneaks back into a model that otherwise looked clean. See the causality contract.

make_features#

Turn an OHLC frame into a causal feature matrix in one call: returns, rolling volatility, momentum, RSI, ATR%, distance-from-moving-average z-scores, range%, calendar dummies, and lagged returns — every column shifted one bar. Internally it reuses atr, rsi, and zscore from edgekit.indicators, so a feature means exactly what the indicator of the same name means.

make_features(df, windows=(5, 20, 60), lags=(1, 2, 3, 5), calendar=True) -> pd.DataFrame

Param	Type	Default	Meaning
`df`	DataFrame	—	OHLC frame (open/high/low/close, optionally volume).
`windows`	tuple	(5, 20, 60)	Lookbacks for rolling vol / momentum / MA-distance.
`lags`	tuple	(1, 2, 3, 5)	Return lags to include as columns.
`calendar`	bool	True	Add day-of-week / month calendar features.

Returns: a pandas DataFrame of causal features indexed like df (every column shifted one bar; leading rows contain NaN until the longest window fills — drop or impute before modelling).

X = ek.timeseries.make_features(bars, windows=(10, 50), lags=(1, 2, 3))
X = X.dropna()   # trim the warm-up rows

add_lags#

Append lagged copies of chosen columns to a frame — the low-level primitive behind the lagsargument of make_features, handy when you want lagged versions of a specific signal you built yourself. Every added column is shifted, so it stays causal.

add_lags(df, cols, lags=(1, 2, 3)) -> pd.DataFrame

Param	Type	Default	Meaning
`df`	DataFrame	—	The frame to extend.
`cols`	list[str]	—	Column names to lag.
`lags`	tuple	(1, 2, 3)	Which lags to add per column.

Returns: a new DataFrame with one extra column per (col, lag) pair.

df = ek.timeseries.add_lags(df, ["ret", "vol"], lags=(1, 2, 5))

Dynamic hedging & regime models#

The estimators above characterise a single series. This block models how two series co-move and how a series switches regime: a Kalman-filtered dynamic hedge ratio, the Ornstein-Uhlenbeck fit that a mean-reverting spread lives or dies on, a GARCH(1,1) volatility model, and a two-state Gaussian HMM. All numpy — no filterpy, arch, or hmmlearn.

→These are fitted models — check convergence, don't trust blindly

Unlike the closed-form diagnostics above, these run iterative estimation (Kalman recursion, MLE, Baum-Welch). Sanity-check the outputs: a GARCH persistence at or above 1 is a non-stationary fit, an OU half_life longer than your holding horizon means the spread reverts too slowly to trade, and HMM state labels are arbitrary (state 0 vs 1 is not fixed across runs — key off means / vars).

kalman_hedge#

A Kalman filter that estimates a time-varying hedge ratio between two series — the dynamic-beta replacement for a static OLS hedge. As the relationship drifts, beta tracks it, and spread is the residual you trade for mean reversion. The workhorse behind a pairs / stat-arb spread.

kalman_hedge(y, x, delta=1e-4, r=1.0) -> dict

Param	Type	Default	Meaning
`y`	array (T,)	—	The dependent series (the leg you hedge).
`x`	array (T,)	—	The explanatory series (the hedging leg).
`delta`	float	1e-4	State-transition noise; larger = beta adapts faster.
`r`	float	1.0	Observation-noise variance.

Returns: a dict with keys "beta" (the time-varying hedge ratio, per bar), "alpha" (the time-varying intercept), and "spread" (the residual y − (alpha + beta·x) to trade).

import edgekit as ek
kf = ek.timeseries.kalman_hedge(y, x, delta=1e-5)
z = ek.indicators.zscore(kf["spread"], 20)   # trade the spread when |z| is large

ou_params#

Fit an Ornstein-Uhlenbeck process to a spread and read off its mean-reversion speed. The half_life — how long it takes the spread to close half the gap to its mean — is the single number that tells you whether a spread is tradeable and on what horizon.

ou_params(spread) -> dict

spread — a mean-reverting series (e.g. kalman_hedge()["spread"]).

Returns: a dict with keys "kappa" (mean-reversion speed), "theta" (long-run mean), "sigma" (instantaneous volatility), and "half_life" (ln(2)/kappa, in bars).

ou = ek.timeseries.ou_params(kf["spread"])
ou["half_life"]   # bars to revert halfway; too long => not tradeable

garch11#

A GARCH(1,1) volatility model fit by maximum likelihood — captures volatility clustering (calm begets calm, shocks beget shocks) and gives a one-step-ahead vol forecast. Use it for conditional vol-targeting and risk overlays where a flat rolling window reacts too slowly.

garch11(rets, max_iter=200) -> dict

Param	Type	Default	Meaning
`rets`	array (T,)	—	The return series to model.
`max_iter`	int	200	Maximum likelihood-optimiser iterations.

Returns: a dict with keys "omega", "alpha", "beta" (the three GARCH parameters), "cond_vol" (the fitted conditional-volatility path), "forecast" (next-step vol), and "persistence" (alpha + beta; approaches 1 as shocks decay slowly).

g = ek.timeseries.garch11(rets)
g["forecast"], g["persistence"]   # persistence >= 1 is a degenerate (non-stationary) fit

hmm_two_state#

A two-state Gaussian hidden Markov model fit by Baum-Welch — the simplest regime detector, splitting a series into (typically) a calm state and a turbulent one. Use the decoded states to gate a strategy on or off by regime.

hmm_two_state(x, iters=100, seed=0) -> dict

Param	Type	Default	Meaning
`x`	array (T,)	—	The observed series (usually returns).
`iters`	int	100	Baum-Welch (EM) iterations.
`seed`	int	0	RNG seed for the initial parameter guess.

Returns: a dict with keys "states" (the most-likely state per bar, 0/1), "means" and "vars" (per-state Gaussian mean and variance), "trans" (the 2×2 transition matrix), and "prob" (the smoothed probability of each state per bar).

!State labels are not stable across runs

Which regime is called 0 vs 1 depends on the random init — always identify the turbulent regime by its higher variance in vars, never by its index.

h = ek.timeseries.hmm_two_state(rets, iters=200)
turbulent = int(np.argmax(h["vars"]))   # the high-vol state
trade_on = h["states"] != turbulent     # sit out the turbulent regime