This Python 3 function uses the standard library and implements the drift-allowed, debiased, overlapping robust test. It accepts log-price levels and differences them once. Feeding it returns would instead test changes in returns. For raw prices, first verify P > 0 and calculate log(P). The function rejects missing values and constant or numerically unresolved returns; it never deletes observations or fills gaps.
Complete synthetic inputs and estimator. All printed p-values are asymptotic arithmetic illustrations.python
from itertools import accumulate
from math import erfc, fsum, isfinite, sqrt
def variance_ratio(log_prices, q):
x = [float(value) for value in log_prices]
n = len(x) - 1 # n returns, n + 1 log prices
if isinstance(q, bool) or not isinstance(q, int) or not 2 <= q < n:
raise ValueError('Require integer 2 <= q < number of returns')
if not all(isfinite(value) for value in x):
raise ValueError('Missing or nonfinite log price')
r = [x[t] - x[t - 1] for t in range(1, n + 1)]
if not all(isfinite(value) for value in r):
raise ValueError('Return overflow')
mu = fsum(r) / n
e = [value - mu for value in r]
scale = max(abs(value) for value in e)
if scale <= 1e-12 * max(1.0, max(abs(value) for value in r)):
raise ValueError('Constant or numerically unresolved returns')
# Common scaling cancels in VR and delta; avoids fourth-power overflow.
u = [value / scale for value in e]
squares = [value * value for value in u]
s = fsum(squares)
v1 = s / (n - 1)
blocks = [fsum(u[start:start + q]) for start in range(n - q + 1)]
m = q * (n - q + 1) * (1 - q / n)
vq = fsum(value * value for value in blocks) / m
vr = vq / v1
# delta has NO n factor here; theta estimates Var(VR - 1).
delta = [fsum(squares[t] * squares[t - j] for t in range(j, n))
/ s**2 for j in range(1, q)]
theta = fsum((2 * (q - j) / q)**2 * delta[j - 1]
for j in range(1, q))
if not isfinite(theta) or theta <= 0:
raise ValueError('Undefined test variance')
se = sqrt(theta)
z = (vr - 1) / se
p = erfc(abs(z) / sqrt(2)) # two-sided asymptotic normal p-value
se_homo = sqrt(2 * (2 * q - 1) * (q - 1) / (3 * q * n))
return dict(n=n, blocks=len(blocks), m=m, vr=vr, theta=theta,
se=se, z=z, p=p, z_homo=(vr - 1) / se_homo)
# Complete synthetic log-return sample: each integer is 0.001 log units.
units = [2, 1, -1, 3, 2, -2, 1, 0, 4, 3, -3, -1, 2, 1, -2, 0,
8, 5, -6, -4, 7, 3, -5, -2, 6, 4, -7, -3, 5, 2, -4, -1]
log_prices = [4.0] + [4.0 + value / 1000 for value in accumulate(units)]
if __name__ == '__main__':
print('log_prices:', log_prices)
for horizon in (2, 4, 8):
print(horizon, variance_ratio(log_prices, horizon))
