import numpy as np import pandas as pd import pytest from quantcore.models.covariance import CovarianceModel, cholesky_factor, symmetric_inv_sqrt, symmetric_sqrt from quantcore.models.simulation import EmpiricalInnovationSimulator, GBMSimulator @pytest.fixture def r(): rng = np.random.default_rng(7) L = np.array([[1.0, 0, 0], [0.6, 0.8, 0], [0.3, 0.2, 0.9]]) x = rng.standard_t(6, size=(1500, 3)) @ L.T * 0.01 return pd.DataFrame(x, index=pd.bdate_range("2019-01-01", periods=1500), columns=list("ABC")) def test_matrix_roots(r): cov = r.cov().values s = symmetric_sqrt(cov) np.testing.assert_allclose(s @ s, cov, atol=1e-12) np.testing.assert_allclose(symmetric_inv_sqrt(cov) @ s, np.eye(3), atol=1e-8) L = cholesky_factor(cov) np.testing.assert_allclose(L @ L.T, cov, atol=1e-12) def test_gbm_scenarios_have_target_moments(r): cm = CovarianceModel("rolling", window=500) mu = pd.Series([0.001, 0.0, -0.001], index=list("ABC")) scen = GBMSimulator(cm).simulate(r, n_scenarios=100_000, horizon_periods=5, mu=mu, seed=1) target = cm.estimate(r).values * 5 # Sobol' sets are rounded up to a power of two; moment matching makes the covariance exact assert scen.model == "gbm" and scen.n_scenarios == 131_072 and scen.horizon_periods == 5 assert scen.extra["n_requested"] == 100_000 and scen.extra["sampling"] == "sobol" np.testing.assert_allclose(np.cov(scen.log_returns, rowvar=False, ddof=0), target, rtol=1e-8, atol=1e-14) # arithmetic drift convention: the expected simple return over the horizon is mu * h (Ito correction) np.testing.assert_allclose(scen.expected_simple, mu.values * 5) np.testing.assert_allclose(scen.mean_simple_returns(), mu.values * 5, atol=3e-5) # log convention: the log-drift is mu * h exactly scen_log = GBMSimulator(cm).simulate(r, n_scenarios=4096, horizon_periods=5, mu=mu, seed=1, drift_convention="log") np.testing.assert_allclose(scen_log.log_returns.mean(axis=0), mu.values * 5, atol=1e-12) assert (scen_log.expected_simple > np.expm1(mu.values * 5)).all() # sigma^2/2 bonus under the log convention # plain Monte Carlo without moment matching: the old behaviour, moments only approximately right mc = GBMSimulator(cm).simulate(r, n_scenarios=100_000, horizon_periods=5, mu=mu, seed=1, sampling="mc", moment_match=False) assert mc.n_scenarios == 100_000 np.testing.assert_allclose(np.cov(mc.log_returns, rowvar=False), target, rtol=0.05, atol=1e-6) def test_innovations_are_standardised(r): sim = EmpiricalInnovationSimulator(CovarianceModel("ewm", lam=0.97)) eps = sim.innovations(r) assert eps.shape[1] == 3 and len(eps) > 1000 cov = np.cov(eps.values, rowvar=False) np.testing.assert_allclose(cov, np.eye(3), atol=0.15) def test_empirical_scenarios_match_latest_covariance(r): cm = CovarianceModel("ewm", lam=0.97) sim = EmpiricalInnovationSimulator(cm) scen = sim.simulate(r, n_scenarios=60_000, horizon_periods=3, seed=2, sampling="mc", moment_match=False) target = cm.estimate(r).values * 3 np.testing.assert_allclose(np.cov(scen.log_returns, rowvar=False), target, rtol=0.2, atol=2e-6) # zero expected return: the log-drift is the Ito correction, the mean simple return ~ 0 np.testing.assert_allclose(scen.mean_simple_returns(), 0.0, atol=5e-4) exact = sim.simulate(r, n_scenarios=60_000, horizon_periods=3, seed=2) # Sobol' ranks + moment matching assert exact.n_scenarios == 65_536 np.testing.assert_allclose(np.cov(exact.log_returns, rowvar=False, ddof=0), target, rtol=1e-8, atol=1e-14) # fat tails survive the filter: kurtosis of the innovations above the Gaussian 0 from scipy.stats import kurtosis assert kurtosis(sim.innovations(r).values).mean() > 0.3 def test_cholesky_root_variant(r): sim = EmpiricalInnovationSimulator(CovarianceModel("rolling", window=250), root="cholesky") eps = sim.innovations(r) np.testing.assert_allclose(np.cov(eps.values, rowvar=False), np.eye(3), atol=0.15) scen = sim.simulate(r, n_scenarios=5000, horizon_periods=1, seed=3, sampling="mc") assert scen.extra["root"] == "cholesky" and scen.log_returns.shape == (5000, 3) def test_restandardisation_reproduces_model_covariance(r): cm = CovarianceModel("ewm", lam=0.97, shrinkage=0.0) for mode in ("diagonal", "full"): sim = EmpiricalInnovationSimulator(cm, standardize=mode) scen = sim.simulate(r, n_scenarios=80_000, horizon_periods=1, seed=5, sampling="mc", moment_match=False) model = cm.estimate(r).values sc = np.cov(scen.log_returns, rowvar=False) np.testing.assert_allclose(np.sqrt(np.diag(sc)), np.sqrt(np.diag(model)), rtol=0.03) assert scen.extra["standardize"] == mode and scen.extra["innovation_variance_median"] > 0 with pytest.raises(ValueError): EmpiricalInnovationSimulator(cm, standardize="bogus") def test_auto_shrinkage_intensity_properties(): from quantcore.measurement.estimators import shrinkage_intensity rng = np.random.default_rng(11) assert shrinkage_intensity(rng.standard_normal((80, 40))) == 1.0 # independent -> full shrinkage L = np.linalg.cholesky(0.4 * np.eye(10) + 0.6) strong = shrinkage_intensity(rng.standard_normal((2000, 10)) @ L.T) assert 0.0 <= strong < 0.05 # strongly correlated, long sample -> little shrinkage cm = CovarianceModel("rolling", window=100, shrinkage="auto") x = pd.DataFrame(rng.standard_normal((300, 20)) @ np.linalg.cholesky(0.5 * np.eye(20) + 0.5).T) assert 0 < cm.intensity(x) < 1 and cm.describe()["shrinkage"] == "auto" def test_importance_sampling_tilt_both_models(r): from quantcore.models.simulation import TiltSpec cm = CovarianceModel("ewm", lam=0.97) exposure = np.array([30_000.0, 50_000.0, 20_000.0]) beta = 0.99 for sim in (GBMSimulator(cm), EmpiricalInnovationSimulator(cm)): plain = sim.simulate(r, n_scenarios=4096, horizon_periods=5, seed=9) tilted = sim.simulate(r, n_scenarios=4096, horizon_periods=5, seed=9, tilt=TiltSpec.at_confidence(exposure, beta, mixture=0.5)) assert plain.probabilities is None and tilted.probabilities is not None assert abs(tilted.probabilities.sum() - 1) < 1e-12 and tilted.describe()["effective_sample_size"] > 2000 # exact first two moments in both cases (weighted for the tilted set) w = tilted.weights np.testing.assert_allclose(w @ tilted.log_returns, plain.log_returns.mean(axis=0), atol=1e-12) cov_t = (tilted.log_returns - w @ tilted.log_returns).T @ ((tilted.log_returns - w @ tilted.log_returns) * w[:, None]) np.testing.assert_allclose(cov_t, np.cov(plain.log_returns, rowvar=False, ddof=0), rtol=1e-8, atol=1e-14) # many more scenarios beyond the reference portfolio's VaR, with the right total probability loss_p = -(plain.simple_returns() @ exposure) loss_t = -(tilted.simple_returns() @ exposure) var = np.quantile(loss_p, beta) assert (loss_t > var).sum() > 5 * (loss_p > var).sum() assert abs(w[loss_t > var].sum() - (1 - beta)) < 0.01 assert "importance_sampling" in tilted.extra and tilted.hash != plain.hash # the default tilt mixes z_0.95 and z_0.99 blocks, so one set serves both levels both = sim.simulate(r, n_scenarios=4096, horizon_periods=5, seed=9, tilt=TiltSpec(exposure=exposure)) w2 = both.weights loss_b = -(both.simple_returns() @ exposure) for level in (0.95, 0.99): var_l = np.quantile(loss_p, level) assert (loss_b > var_l).sum() > 3 * (loss_p > var_l).sum() assert abs(w2[loss_b > var_l].sum() - (1 - level)) < 0.01 # determinism: same seed -> same hash, other seed -> different a = GBMSimulator(cm).simulate(r, n_scenarios=512, seed=1) b = GBMSimulator(cm).simulate(r, n_scenarios=512, seed=1) c = GBMSimulator(cm).simulate(r, n_scenarios=512, seed=2) assert a.hash == b.hash != c.hash def test_scenario_distribution_reports_the_model_quantities(r): sim = GBMSimulator(CovarianceModel("ewm", lam=0.97, shrinkage=0.0)) scen = sim.simulate(r, n_scenarios=8192, horizon_periods=5, mu=pd.Series(0.0004, index=r.columns), seed=5) v = np.array([400_000.0, 300_000.0, 300_000.0]) d = scen.distribution(v, bins=40, confidence=(0.95, 0.99)) cov_h = CovarianceModel("ewm", lam=0.97, shrinkage=0.0).estimate(r).values * 5 sigma = np.sqrt(v @ cov_h @ v) assert d["volatility"] == pytest.approx(sigma, rel=0.02) # the reported mean is the mean of the *displayed* sample; moment matching is exact in log # space, so it equals the model's expectation only up to the convexity of exp assert d["expected_pnl"] == pytest.approx(float(scen.expected_simple @ v), rel=1e-4) assert sum(d["histogram"]["mass"]) == pytest.approx(1.0) assert len(d["histogram"]["edges"]) == 41 # Gaussian reference: VaR = z sigma - mean, CVaR = phi(z)/(1-beta) sigma - mean from scipy.stats import norm for beta in (0.95, 0.99): z = norm.ppf(beta) assert d["risk"][f"{beta:g}"]["var"] == pytest.approx(z * sigma - d["expected_pnl"], rel=0.03) assert d["risk"][f"{beta:g}"]["cvar"] == pytest.approx( norm.pdf(z) / (1 - beta) * sigma - d["expected_pnl"], rel=0.03) assert abs(d["skewness"]) < 0.25 and d["kurtosis"] == pytest.approx(3.0, abs=0.4) def test_a_tilted_set_reports_the_same_distribution_as_an_untilted_one(r): from quantcore.models.simulation import TiltSpec sim = GBMSimulator(CovarianceModel("ewm", lam=0.97, shrinkage=0.0)) v = np.array([400_000.0, 300_000.0, 300_000.0]) kw = dict(n_scenarios=8192, horizon_periods=5, mu=pd.Series(0.0004, index=r.columns), seed=11) plain = sim.simulate(r, **kw).distribution(v) tilted = sim.simulate(r, **kw, tilt=TiltSpec(exposure=v)).distribution(v) assert tilted["weighted"] and tilted["effective_sample_size"] < 8192 assert tilted["expected_pnl"] == pytest.approx(plain["expected_pnl"], rel=1e-6) assert tilted["volatility"] == pytest.approx(plain["volatility"], rel=0.02) assert tilted["risk"]["0.99"]["cvar"] == pytest.approx(plain["risk"]["0.99"]["cvar"], rel=0.03) for q in ("0.05", "0.5"): assert tilted["quantiles"][q] == pytest.approx(plain["quantiles"][q], rel=0.05) def test_summary_reports_price_returns_that_cannot_fall_below_minus_one(r): """A long horizon makes the two spaces look very different: a log-return of -2 is -86 % of the price. The table must carry both, and the price column must be bounded below by -100 %.""" sim = GBMSimulator(CovarianceModel("ewm", lam=0.97, shrinkage=0.0)) scen = sim.simulate(r, n_scenarios=4096, horizon_periods=63, mu=pd.Series(0.0004, index=r.columns), seed=3) t = scen.summary() for col in ("q01", "q05", "q50", "q95", "min", "max"): np.testing.assert_allclose(t[f"{col}_simple"], np.expm1(t[col]), rtol=1e-12) # a monotone map assert (t["min_simple"] > -1).all() and (t["min"] < t["min_simple"]).all() assert (t["max_simple"] > t["max"]).all() # the scenario mean in price space is above the log mean (Jensen) and matches the model's assert (t["mean_simple"] > t["mean"]).all() assert np.allclose(t["mean_simple"], t["expected_simple"], rtol=0.05)