b5da424dfe
The portfolio data, risk and mean-CVaR optimisation service, prepared for students who run it on their own DGX Spark without root, Docker or Tailscale. - deploy/apptainer/qc.sh: install, start/stop/status/logs, cli/api/python, data import/export and Yahoo setup/update. cuOpt and quantcore run as two Apptainer instances of the user and talk over 127.0.0.1; the API and the test bench listen on 127.0.0.1 only. Builds quantcore.sif without root, or falls back to a plain Python image plus a virtualenv when the Spark allows no unprivileged builds. Detects taken ports and suggests free ones. - deploy/apptainer/quantcore.def + constraints.txt: the image, with the library versions of the verified production container pinned. - cuOpt under Apptainer: NVIDIA's image sets HOME=/opt/cuopt, read-only in a SIF, so the solver died on its first job; qc.sh gives it a writable home. - solvers: a cuOpt that disappears during a solve now falls back to HiGHS instead of answering 500 (requests' ConnectionError escaped the fallback). - client: idempotent GETs retry a dropped keep-alive connection, which made `quantcore-api update --wait` fail when a poll met uvicorn's 5 s timeout. - config: cuOpt defaults to http://127.0.0.1:5000 and `quantcore serve` to 127.0.0.1; the Docker images keep passing their own values. - docs/SPARK_APPTAINER.md: the students' guide; README, USAGE, ARCHITECTURE and SPARK_DEPLOYMENT point to it. Tested on a DGX Spark (GB10, Apptainer 1.5.4, unprivileged user) in both modes: selftest on every endpoint, optimisation on cuOpt, Yahoo update, port conflicts; 139 tests pass in the container against the live cuOpt. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
198 lines
11 KiB
Python
198 lines
11 KiB
Python
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)
|