Coverage for src/jquantstats/_stats/_montecarlo.py: 100%

71 statements  

« prev     ^ index     » next       coverage.py v7.17.1, created at 2026-11-03 15:52 +0000

1"""Monte Carlo simulation statistics.""" 

2 

3from __future__ import annotations 

4 

5import math 

6from typing import TYPE_CHECKING 

7 

8import numpy as np 

9import polars as pl 

10 

11from ..exceptions import NonPositivePeriodsPerYearError, NonPositiveWindowError 

12 

13if TYPE_CHECKING: 

14 from ..data import Data 

15 

16 

17class _MonteCarloStatsMixin: 

18 """Mixin providing Monte Carlo simulation statistics for return series.""" 

19 

20 _data: Data 

21 all: pl.DataFrame 

22 

23 if TYPE_CHECKING: 

24 from .._protocol import DataLike 

25 

26 data: DataLike 

27 

28 @staticmethod 

29 def _validate_positive_integer(name: str, value: int) -> int: 

30 """Validate that *value* is a positive integer.""" 

31 if isinstance(value, bool) or not isinstance(value, int) or value <= 0: 

32 raise NonPositiveWindowError(name) 

33 return value 

34 

35 @staticmethod 

36 def _block_bootstrap_paths(values: np.ndarray, n: int, period: int, block_size: int) -> np.ndarray: 

37 """Generate *n* return paths via circular block bootstrap, returning an ``(n, period)`` array. 

38 

39 Block starts are drawn uniformly over all observations and a block that 

40 runs past the last observation wraps around to the first (Politis and 

41 Romano, 1992). Every observation is therefore drawn with equal 

42 probability. A non-circular moving-block bootstrap would draw the first 

43 and last observations at ``1 / block_size`` of the rate of interior 

44 ones, biasing any statistic of a series whose behaviour changes through 

45 the sample. 

46 

47 All *n* paths are sampled simultaneously using fully vectorised numpy 

48 indexing — no Python-level loop is required. 

49 

50 Args: 

51 values: 1-D array of historical returns. 

52 n: Number of simulation paths. 

53 period: Length of each simulated path (number of return observations). 

54 block_size: Block length for the block bootstrap resampling. 

55 

56 Returns: 

57 np.ndarray: Float64 array of shape ``(n, period)``. 

58 

59 """ 

60 n_obs = values.size 

61 if n_obs == 0: # pragma: no cover 

62 return np.full((n, period), np.nan, dtype=np.float64) 

63 

64 n_blocks = math.ceil(period / block_size) 

65 # Draw all block starting indices at once: shape (n, n_blocks) 

66 starts = np.random.randint(0, n_obs, size=(n, n_blocks)) 

67 # Build full index array: (n, n_blocks, block_size), wrapping past the end 

68 idx = (starts[:, :, np.newaxis] + np.arange(block_size)[np.newaxis, np.newaxis, :]) % n_obs 

69 # Flatten and trim to the requested period length 

70 return np.asarray(values[idx].reshape(n, -1)[:, :period]) 

71 

72 def _simulate_distribution(self, n: int, period: int) -> dict[str, np.ndarray]: 

73 """Prepare validated inputs and sample *n* block-bootstrap paths per asset. 

74 

75 Returns a dict mapping each asset column name to an ``(n, period)`` 

76 float64 array of simulated return paths. Assets with no usable data 

77 map to an ``(n, period)`` array of NaN. 

78 

79 Args: 

80 n: Number of simulation paths (must be a positive integer). 

81 period: Path length in return observations (must be a positive integer). 

82 

83 Returns: 

84 dict[str, np.ndarray]: Asset name → ``(n, period)`` paths array. 

85 

86 """ 

87 n = self._validate_positive_integer("n", n) 

88 period = self._validate_positive_integer("period", period) 

89 block_size = max(1, round(period**0.5)) 

90 

91 paths: dict[str, np.ndarray] = {} 

92 for col, series in self._data.items(): 

93 clean = series.cast(pl.Float64).drop_nulls().drop_nans() 

94 values = np.asarray(clean.to_numpy(), dtype=np.float64) 

95 if values.size == 0: 

96 paths[col] = np.full((n, period), np.nan, dtype=np.float64) 

97 else: 

98 block = min(block_size, values.size) 

99 paths[col] = self._block_bootstrap_paths(values, n, period, block) 

100 return paths 

101 

102 def montecarlo(self, n: int = 1000, period: int = 252) -> pl.DataFrame: 

103 """Simulate cumulative returns across *n* block-bootstrap paths. 

104 

105 Args: 

106 n: Number of Monte Carlo paths. Defaults to 1000. 

107 period: Simulation horizon in return observations. Defaults to 252. 

108 

109 Returns: 

110 pl.DataFrame: Shape ``(n, n_assets)`` — one simulated terminal 

111 cumulative return per path and asset. 

112 

113 

114 Returns NaN when: 

115 Entries are NaN for assets with no usable (non-null, non-NaN) 

116 observations. 

117 """ 

118 paths = self._simulate_distribution(n=n, period=period) 

119 result = {col: np.prod(1.0 + arr, axis=1) - 1.0 for col, arr in paths.items()} 

120 return pl.DataFrame(result) 

121 

122 def montecarlo_sharpe( 

123 self, 

124 n: int = 1000, 

125 period: int = 252, 

126 periods_per_year: int | float | None = None, 

127 ) -> pl.DataFrame: 

128 """Simulate the Sharpe-ratio distribution across block-bootstrap paths. 

129 

130 Args: 

131 n: Number of Monte Carlo paths. Defaults to 1000. 

132 period: Simulation horizon in return observations. Defaults to 252. 

133 periods_per_year: Annualisation factor. Defaults to the value 

134 inferred from the data. 

135 

136 Returns: 

137 pl.DataFrame: Shape ``(n, n_assets)`` — one simulated annualised 

138 Sharpe ratio per path and asset. 

139 

140 

141 Returns NaN when: 

142 Entries are NaN when a path's standard deviation is zero or the asset 

143 has no usable observations. 

144 """ 

145 ppy = self._data._periods_per_year if periods_per_year is None else periods_per_year 

146 if ppy <= 0: 

147 raise NonPositivePeriodsPerYearError 

148 scale = math.sqrt(ppy) 

149 paths = self._simulate_distribution(n=n, period=period) 

150 result: dict[str, np.ndarray] = {} 

151 for col, arr in paths.items(): 

152 means = arr.mean(axis=1) 

153 stds = arr.std(axis=1, ddof=1) 

154 with np.errstate(invalid="ignore", divide="ignore"): 

155 result[col] = np.where(stds == 0.0, np.nan, means / stds * scale) 

156 return pl.DataFrame(result) 

157 

158 def montecarlo_drawdown(self, n: int = 1000, period: int = 252) -> pl.DataFrame: 

159 """Simulate the maximum-drawdown distribution across block-bootstrap paths. 

160 

161 Args: 

162 n: Number of Monte Carlo paths. Defaults to 1000. 

163 period: Simulation horizon in return observations. Defaults to 252. 

164 

165 Returns: 

166 pl.DataFrame: Shape ``(n, n_assets)`` — one simulated maximum 

167 drawdown per path and asset (values in ``[-1, 0]``). 

168 

169 

170 Returns NaN when: 

171 Entries are NaN for assets with no usable (non-null, non-NaN) 

172 observations. 

173 """ 

174 paths = self._simulate_distribution(n=n, period=period) 

175 result: dict[str, np.ndarray] = {} 

176 for col, arr in paths.items(): 

177 nav = np.cumprod(1.0 + arr, axis=1) 

178 hwm = np.maximum.accumulate(nav, axis=1) 

179 result[col] = np.min(nav / hwm - 1.0, axis=1) 

180 return pl.DataFrame(result) 

181 

182 def montecarlo_cagr( 

183 self, 

184 n: int = 1000, 

185 period: int = 252, 

186 periods_per_year: int | float | None = None, 

187 ) -> pl.DataFrame: 

188 """Simulate the CAGR distribution across block-bootstrap paths. 

189 

190 Args: 

191 n: Number of Monte Carlo paths. Defaults to 1000. 

192 period: Simulation horizon in return observations. Defaults to 252. 

193 periods_per_year: Annualisation factor. Defaults to the value 

194 inferred from the data. 

195 

196 Returns: 

197 pl.DataFrame: Shape ``(n, n_assets)`` — one simulated annualised 

198 CAGR per path and asset. 

199 

200 

201 Returns NaN when: 

202 Entries are NaN when a path's terminal compound return is non-positive 

203 or the asset has no usable observations. 

204 """ 

205 ppy = self._data._periods_per_year if periods_per_year is None else periods_per_year 

206 if ppy <= 0: 

207 raise NonPositivePeriodsPerYearError 

208 years = period / ppy 

209 paths = self._simulate_distribution(n=n, period=period) 

210 result: dict[str, np.ndarray] = {} 

211 for col, arr in paths.items(): 

212 totals = np.prod(1.0 + arr, axis=1) 

213 with np.errstate(invalid="ignore"): 

214 result[col] = np.where(totals > 0, totals ** (1.0 / years) - 1.0, np.nan) 

215 return pl.DataFrame(result)