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

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 block bootstrap, returning an ``(n, period)`` array. 

38 

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

40 indexing — no Python-level loop is required. 

41 

42 Args: 

43 values: 1-D array of historical returns. 

44 n: Number of simulation paths. 

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

46 block_size: Block length for the block bootstrap resampling. 

47 

48 Returns: 

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

50 

51 """ 

52 n_obs = values.size 

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

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

55 

56 n_blocks = math.ceil(period / block_size) 

57 max_start = max(1, n_obs - block_size + 1) 

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

59 starts = np.random.randint(0, max_start, size=(n, n_blocks)) 

60 # Build full index array: (n, n_blocks, block_size) 

61 idx = starts[:, :, np.newaxis] + np.arange(block_size)[np.newaxis, np.newaxis, :] 

62 idx = np.clip(idx, 0, n_obs - 1) 

63 # Flatten and trim to the requested period length 

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

65 

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

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

68 

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

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

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

72 

73 Args: 

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

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

76 

77 Returns: 

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

79 

80 """ 

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

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

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

84 

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

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

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

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

89 if values.size == 0: 

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

91 else: 

92 block = min(block_size, values.size) 

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

94 return paths 

95 

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

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

98 

99 Args: 

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

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

102 

103 Returns: 

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

105 cumulative return per path and asset. 

106 

107 

108 Returns NaN when: 

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

110 observations. 

111 """ 

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

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

114 return pl.DataFrame(result) 

115 

116 def montecarlo_sharpe( 

117 self, 

118 n: int = 1000, 

119 period: int = 252, 

120 periods_per_year: int | float | None = None, 

121 ) -> pl.DataFrame: 

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

123 

124 Args: 

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

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

127 periods_per_year: Annualisation factor. Defaults to the value 

128 inferred from the data. 

129 

130 Returns: 

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

132 Sharpe ratio per path and asset. 

133 

134 

135 Returns NaN when: 

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

137 has no usable observations. 

138 """ 

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

140 if ppy <= 0: 

141 raise NonPositivePeriodsPerYearError 

142 scale = math.sqrt(ppy) 

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

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

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

146 means = arr.mean(axis=1) 

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

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

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

150 return pl.DataFrame(result) 

151 

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

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

154 

155 Args: 

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

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

158 

159 Returns: 

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

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

162 

163 

164 Returns NaN when: 

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

166 observations. 

167 """ 

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

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

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

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

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

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

174 return pl.DataFrame(result) 

175 

176 def montecarlo_cagr( 

177 self, 

178 n: int = 1000, 

179 period: int = 252, 

180 periods_per_year: int | float | None = None, 

181 ) -> pl.DataFrame: 

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

183 

184 Args: 

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

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

187 periods_per_year: Annualisation factor. Defaults to the value 

188 inferred from the data. 

189 

190 Returns: 

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

192 CAGR per path and asset. 

193 

194 

195 Returns NaN when: 

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

197 or the asset has no usable observations. 

198 """ 

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

200 if ppy <= 0: 

201 raise NonPositivePeriodsPerYearError 

202 years = period / ppy 

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

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

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

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

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

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

209 return pl.DataFrame(result)