# -*- coding: utf-8 -*- """ Variance Ratio Test + Chow-Denning Joint Test for Japanese Equities ===================================================================== Implements two related tests for the random walk hypothesis on log-prices: 1. Lo-MacKinlay (1988) individual Variance Ratio (VR) test -------------------------------------------------------- For each aggregation interval q: VR(q) = Var(r_q) / (q * Var(r_1)) [should equal 1.0 under H0] Z1(q) = sqrt(T) * (VR(q) - 1) / sqrt(2*(q-1)) [homoskedastic] Z*(q) = sqrt(T) * (VR(q) - 1) / sqrt(theta(q)) [heteroskedastic] Under H0 (random walk): Z(q) --d--> N(0,1) as T → ∞ 2. Chow-Denning (1993) joint / multiple comparison test ------------------------------------------------------ Tests H0: VR(q) = 1 for ALL q simultaneously, controlling the familywise error rate via the Studentised Maximum Modulus (SMM) distribution. CD = max_q |Z*(q)| Reject H0 at level α if CD > SMM critical value (k, ∞, α). Using the heteroskedastic Z* is standard practice because it is robust to conditional heteroskedasticity (ARCH effects), which are common in equity returns. References ---------- Lo, A. W., & MacKinlay, A. C. (1988). "Stock Market Prices Do Not Follow Random Walks." Review of Financial Studies, 1(1), 41-66. Chow, K. V., & Denning, K. C. (1993). "A Simple Multiple Variance Ratio Test." Journal of Econometrics, 58(3), 385-401. Author : cting Created: Wed May 27 16:12:46 2026 """ import pandas as pd import numpy as np import matplotlib.pyplot as plt from scipy.stats import norm # --------------------------------------------------------------------------- # Helper functions # --------------------------------------------------------------------------- def read_jp(ticker: str, datadir: str) -> pd.DataFrame: """ Load a Japanese stock CSV and return a split-adjusted price series. The CSV must contain columns: Date (YYYYMMDD), Close, Split. Split values are backward-adjustment factors: when a split of factor f occurs on date d, all prices before d are divided by f so the entire series is expressed on the post-split share basis. Parameters ---------- ticker : stock ticker string (used as CSV filename stem) datadir : directory path (trailing backslash required on Windows) Returns ------- df : DataFrame with DatetimeIndex and a 'Price' column of adjusted closes """ filename = datadir + ticker + '.csv' df = pd.read_csv(filename) df.index = pd.to_datetime(df["Date"], format="%Y%m%d") n = len(df) # Cumulative split-adjustment vector (all ones until a split is found). # When split factor f is recorded at row i, multiply all earlier rows # by f so that every price reflects the current share count. spfac = np.ones(n) for i in range(n): if df.Split.iloc[i] != 0.0: print(f"Stock Split detected: {df.Date.iloc[i]} factor={df.Split.iloc[i]}") spfac[:i] *= df.Split.iloc[i] # vectorised slice — no inner loop df['Price'] = df.Close / spfac return df def autocorrelation_lag1(x: np.ndarray) -> float: """ Lag-1 sample autocorrelation using the biased (1/T) denominator, consistent with Lo-MacKinlay (1988) notation. rho_1 = [sum_{t=2}^T (x_t - xbar)(x_{t-1} - xbar)] / [sum_{t=1}^T (x_t - xbar)^2] References ---------- https://en.wikipedia.org/wiki/Autocorrelation#Estimation Parameters ---------- x : 1-D array of returns Returns ------- Lag-1 autocorrelation as a scalar float """ x = np.asarray(x, dtype=float) n = len(x) xc = x - x.mean() # np.correlate 'full' gives lags -(n-1)…0…(n-1); last n = lags 0,1,…,n-1 r = np.correlate(xc, xc, mode='full')[-n:] result = r / (x.var() * np.arange(n, 0, -1)) return float(result[1]) def heteroskedastic_variance(r1: np.ndarray, q: int) -> float: """ Compute theta(q), the heteroskedasticity-consistent asymptotic variance of sqrt(T)*(VR(q)-1), following Lo & MacKinlay (1988) eq. (8). theta(q) = 4 * sum_{k=1}^{q-1} [(q-k)/q]^2 * delta(k) where delta(k) is the heteroskedasticity-robust weight at lag k: delta(k) = [sum_{t=k+1}^{T} (r_t - mu)^2 * (r_{t-k} - mu)^2] / [sum_{t=1}^{T} (r_t - mu)^2]^2 Parameters ---------- r1 : array of single-period log returns, length T q : aggregation interval Returns ------- theta : scalar float (the asymptotic variance) """ T = len(r1) mu = r1.mean() rc = r1 - mu # demeaned returns denom = (rc ** 2).sum() ** 2 # [sum(rc^2)]^2 theta = 0.0 for k in range(1, q): # Numerator: sum of products of squared demeaned returns at lag k numerator = (rc[k:] ** 2 * rc[:T - k] ** 2).sum() delta_k = T * numerator / denom # LM88 scale by T weight = ((q - k) / q) ** 2 theta += weight * delta_k return 4.0 * theta def smm_critical_value(k: int, alpha: float = 0.05) -> float: """ Approximate the two-tailed Studentised Maximum Modulus (SMM) critical value for k tests at significance level alpha, using the Bonferroni bound as a conservative upper bound. The exact SMM quantile requires special tables (Chow & Denning 1993, Table 1). The Bonferroni approximation z_{alpha/(2k)} is standard and slightly conservative, which is acceptable for most empirical work. For exact critical values see: Stoline, M. R., & Ury, H. K. (1979). "Tables of the Studentized Maximum Modulus Distribution." Technometrics, 21(1), 87-93. Parameters ---------- k : number of variance ratios being tested jointly alpha : familywise significance level (default 0.05) Returns ------- cv : Bonferroni-adjusted critical value (scalar float) """ return float(norm.ppf(1 - alpha / (2 * k))) # --------------------------------------------------------------------------- # Configuration # --------------------------------------------------------------------------- DATADIR = 'C:\\YahooJap\\stocks\\' TICKER = '1329' # iShares MSCI Japan ETF (TSE-listed) MAX_Q = 10 # maximum aggregation interval q (inclusive) ALPHA = 0.05 # familywise significance level for Chow-Denning # --------------------------------------------------------------------------- # Data loading and return construction # --------------------------------------------------------------------------- df = read_jp(TICKER, DATADIR) log_price = np.log(df['Price']) # Single-period log returns: r_t = log(P_t) - log(P_{t-1}) r1 = log_price.diff().dropna().values # shape (T,) T = len(r1) mu_1 = r1.mean() var_1 = r1.var(ddof=0) # biased variance (consistent with LM88) # --------------------------------------------------------------------------- # Variance Ratio test loop (q = 2, 3, …, MAX_Q) # --------------------------------------------------------------------------- q_values = range(2, MAX_Q + 1) VR = {} # VR(q) — variance ratio Z_homo = {} # Z1(q) — homoskedastic test statistic (LM88 eq. 5) Z_hetero = {} # Z*(q) — heteroskedastic test statistic (LM88 eq. 8) lag1_auto = {} # lag-1 autocorrelation of q-period returns n_obs = {} # number of non-overlapping q-period observations for q in q_values: # ---- build non-overlapping q-period returns ---- idx = np.arange(0, T + 1, q) # indices into log_price (T+1 obs) p_q = log_price.iloc[idx].values r_q = np.diff(p_q) # q-period log returns M = len(r_q) lag1_auto[q] = autocorrelation_lag1(r_q) # q-period return variance (biased estimator, mean = q*mu_1 under H0) mu_q = q * mu_1 var_q = np.mean((r_q - mu_q) ** 2) # Variance Ratio VR[q] = var_q / (q * var_1) # Homoskedastic Z statistic (LM88, eq. 5) # Assumes iid returns; sensitive to ARCH/GARCH effects Z_homo[q] = np.sqrt(T) * (VR[q] - 1) / np.sqrt(2 * (q - 1)) # Heteroskedastic Z* statistic (LM88, eq. 8) # Robust to conditional heteroskedasticity; used in Chow-Denning (1993) theta = heteroskedastic_variance(r1, q) Z_hetero[q] = np.sqrt(T) * (VR[q] - 1) / np.sqrt(theta) if theta > 0 else 0.0 n_obs[q] = M print(f"q={q:2d} M={M:4d} theta={theta:.6f}") # --------------------------------------------------------------------------- # Chow-Denning (1993) joint test # --------------------------------------------------------------------------- # CD statistic = max over all q of |Z*(q)| # Reject H0: {VR(q)=1 for all q} if CD > SMM critical value k = len(q_values) # number of intervals tested CD = max(abs(Z_hetero[q]) for q in q_values) CD_q_star = max(q_values, key=lambda q: abs(Z_hetero[q])) # q that achieves max cv_5pct = smm_critical_value(k, alpha=0.05) cv_1pct = smm_critical_value(k, alpha=0.01) reject_5 = CD > cv_5pct reject_1 = CD > cv_1pct # --------------------------------------------------------------------------- # Results table # --------------------------------------------------------------------------- print("\n" + "=" * 76) print(f"{'q':>4} {'n_obs':>6} {'AC(1)':>8} {'VR(q)':>8} " f"{'Z_homo':>8} {'Z*_hetero':>10}") print("-" * 76) for q in q_values: star = " <--" if q == CD_q_star else "" print(f"{q:4d} {n_obs[q]:6d} {lag1_auto[q]:8.4f} {VR[q]:8.4f} " f"{Z_homo[q]:8.2f} {Z_hetero[q]:10.2f}{star}") print("=" * 76) print("\nIndividual test critical values (two-tailed):") print(" |Z| > 1.96 → reject H0 at 5% (homoskedastic)") print(" |Z| > 2.58 → reject H0 at 1% (homoskedastic)") print(f"\nChow-Denning Joint Test (k={k} intervals, Bonferroni approximation):") print(f" CD = max|Z*(q)| = {CD:.4f} (achieved at q={CD_q_star})") print(f" Critical values: 5% = {cv_5pct:.4f}, 1% = {cv_1pct:.4f}") print(f" Decision at 5%: {'REJECT H0 — random walk rejected' if reject_5 else 'Fail to reject H0 — consistent with random walk'}") print(f" Decision at 1%: {'REJECT H0 — random walk rejected' if reject_1 else 'Fail to reject H0 — consistent with random walk'}") # --------------------------------------------------------------------------- # Plot # --------------------------------------------------------------------------- qs = list(q_values) vrs = [VR[q] for q in qs] zs = [Z_hetero[q] for q in qs] fig, axes = plt.subplots(1, 2, figsize=(13, 5)) fig.suptitle(f'Variance Ratio & Chow-Denning Test — {TICKER}', fontsize=13) # -- Panel 1: Variance Ratio profile -- ax1 = axes[0] ax1.axhline(1.0, color='grey', linestyle='--', linewidth=1, label='H₀: VR = 1') ax1.plot(qs, vrs, '-^r', linewidth=1.5, markersize=6, label='VR(q)') ax1.set_xlabel('Aggregation interval q (trading days)') ax1.set_ylabel('Variance Ratio VR(q)') ax1.set_title('Lo-MacKinlay Variance Ratios') ax1.legend() ax1.grid(True, alpha=0.3) # -- Panel 2: Heteroskedastic Z* statistics with CD critical value bands -- ax2 = axes[1] ax2.axhline( cv_5pct, color='orange', linestyle='--', linewidth=1, label=f'CD 5% c.v. = ±{cv_5pct:.2f}') ax2.axhline(-cv_5pct, color='orange', linestyle='--', linewidth=1) ax2.axhline( cv_1pct, color='red', linestyle=':', linewidth=1, label=f'CD 1% c.v. = ±{cv_1pct:.2f}') ax2.axhline(-cv_1pct, color='red', linestyle=':', linewidth=1) ax2.axhline(0.0, color='grey', linestyle='-', linewidth=0.8) ax2.bar(qs, zs, color=['tomato' if abs(z) > cv_5pct else 'steelblue' for z in zs], alpha=0.75, label='Z*(q) heteroskedastic') ax2.scatter([CD_q_star], [Z_hetero[CD_q_star]], color='black', zorder=5, label=f'CD statistic = {CD:.2f}') ax2.set_xlabel('Aggregation interval q (trading days)') ax2.set_ylabel('Z* statistic') ax2.set_title('Chow-Denning Test Statistics') ax2.legend(fontsize=8) ax2.grid(True, alpha=0.3) plt.tight_layout() plt.show()