Coverage for src/jsharpe/sharpe/quadrature.py: 100%

1"""Gaussian-expectation quadrature and order-statistic moments. 

2 

3This module groups the Gauss-Hermite expectation helper and the moments 

4of the maximum of k standard normals, which together form the 

5expectation core relied upon by the Sharpe-ratio variance routines. 

6""" 

7# ruff: noqa: N802, N806 

8 

9from collections.abc import Callable 

10 

11import numpy as np 

12import scipy 

13 

14 

15def make_expectation_gh( 

16 n_nodes: int = 200, 

17) -> Callable[[Callable[[np.ndarray], np.ndarray]], float]: 

18 """Create an expectation function using Gauss-Hermite quadrature. 

19 

20 Returns a function that computes E[g(Z)] where Z ~ N(0,1) using 

21 Gauss-Hermite quadrature with the specified number of nodes. 

22 

23 The approximation is: E[g(Z)] ≈ (1/√π) Σ w_i g(√2 t_i) 

24 where (t_i, w_i) are Gauss-Hermite nodes and weights. 

25 

26 Args: 

27 n_nodes: Number of quadrature nodes. Higher values increase accuracy. 

28 Default 200. 

29 

30 Returns: 

31 Function E(g) that computes E[g(Z)] for Z ~ N(0,1). 

32 

33 Example: 

34 >>> E = make_expectation_gh(n_nodes=100) 

35 >>> # E[Z^2] should be 1 for standard normal 

36 >>> bool(abs(E(lambda z: z**2) - 1.0) < 1e-6) 

37 True 

38 >>> # E[Z] should be 0 for standard normal 

39 >>> bool(abs(E(lambda z: z)) < 1e-10) 

40 True 

41 """ 

42 nodes, weights = np.polynomial.hermite.hermgauss(n_nodes) 

43 scale = np.sqrt(2.0) 

44 norm = 1.0 / np.sqrt(np.pi) 

45 x = scale * nodes 

46 

47 def E(g: Callable[[np.ndarray], np.ndarray]) -> float: 

48 """Compute E[g(Z)] for Z ~ N(0,1) via Gauss-Hermite quadrature. 

49 

50 Args: 

51 g: Function to integrate against the standard normal density. 

52 

53 Returns: 

54 Approximation of E[g(Z)]. 

55 """ 

56 vals = g(x) 

57 return float(norm * np.dot(weights, vals)) 

58 

59 return E 

60 

61 

62E_under_normal = make_expectation_gh(n_nodes=200) 

63 

64 

65def moments_Mk(k: int, *, rho: float = 0) -> tuple[float, float, float]: 

66 """Compute moments of the maximum of k standard normal random variables. 

67 

68 Computes E[M_k], E[M_k^2], and Var[M_k] where M_k = max(Z_1, ..., Z_k) 

69 and Z_i are standard normal. Supports both independent (rho=0) and 

70 equi-correlated (rho > 0) cases. 

71 

72 For the correlated case with equi-correlation rho: 

73 Z_i = sqrt(rho) * X + sqrt(1-rho) * Y_i 

74 M = sqrt(rho) * X + sqrt(1-rho) * max(Y_i) 

75 

76 Args: 

77 k: Number of standard normal random variables. 

78 rho: Equi-correlation coefficient in [0, 1). Default is 0 (independent). 

79 

80 Returns: 

81 Tuple of (Ez, Ez2, var): 

82 - Ez: Expected value E[M_k] 

83 - Ez2: Second moment E[M_k^2] 

84 - var: Variance Var[M_k] 

85 

86 Example: 

87 >>> Ez, Ez2, var = moments_Mk(1) 

88 >>> bool(abs(Ez) < 1e-10) # E[max of one N(0,1)] = 0 

89 True 

90 >>> bool(abs(var - 1.0) < 1e-10) # Var[max of one N(0,1)] = 1 

91 True 

92 """ 

93 Phi = scipy.stats.norm.cdf 

94 Ez = E_under_normal(lambda z: k * z * (Phi(z) ** (k - 1))) 

95 Ez2 = E_under_normal(lambda z: k * z**2 * (Phi(z) ** (k - 1))) 

96 var = Ez2 - Ez**2 

97 

98 Ez = (1 - rho) * Ez 

99 var = rho + (1 - rho) * var 

100 Ez2 = var + Ez**2 

101 

102 return Ez, Ez2, var