factor

%23%20%2F%2F%2F%20script%0A%23%20requires-python%20%3D%20%22%3E%3D3.12%22%0A%23%20dependencies%20%3D%20%5B%0A%23%20%20%20%20%20%22marimo%3D%3D0.14.13%22%2C%0A%23%20%20%20%20%20%22numpy%3D%3D2.3.0%22%2C%0A%23%20%20%20%20%20%22plotly%3D%3D6.7.0%22%2C%0A%23%20%20%20%20%20%22polars%3D%3D1.44.1%22%2C%0A%23%20%20%20%20%20%22jquantstats%3D%3D0.11.0%22%2C%0A%23%20%20%20%20%20%22cvx-linalg%3E%3D0.9.3%22%2C%0A%23%20%20%20%20%20%22cvxcla%22%0A%23%20%5D%0A%23%0A%23%20%5Btool.uv.sources%5D%0A%23%20cvxcla%20%3D%20%7B%20path%20%3D%20%22..%2F..%2F..%22%2C%20editable%3Dtrue%20%7D%0A%23%0A%23%20%2F%2F%2F%0A%22%22%22RMT-cleaned%20efficient%20frontier%20with%20the%20factor%20covariance%20backend.%22%22%22%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.25.1%22%0Aapp%20%3D%20marimo.App()%0A%0Awith%20app.setup%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20polars%20as%20pl%0A%20%20%20%20from%20jquantstats%20import%20Data%0A%0A%20%20%20%20from%20cvxcla%20import%20CLA%2C%20FactorCovariance%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20From%20Marchenko-Pastur%20to%20Woodbury%0A%0A%20%20%20%20We%20trace%20the%20**exact**%20efficient%20frontier%20of%20an%20RMT-cleaned%20covariance%0A%20%20%20%20without%20ever%20forming%20an%20%24n%20%5Ctimes%20n%24%20matrix.%0A%0A%20%20%20%201.%20Simulate%20returns%20with%20a%20latent%20factor%20structure.%0A%20%20%20%202.%20Clip%20the%20sample%20eigenvalues%20at%20the%20Marchenko-Pastur%20edge.%0A%20%20%20%203.%20The%20cleaned%20covariance%20is%20diagonal-plus-low-rank%2C%0A%20%20%20%20%20%20%20%24%5CSigma%20%3D%20%5Cbar%7Bd%7D%5C%2C%20I%20%2B%20V_k%20(%5CLambda_k%20-%20%5Cbar%7Bd%7D%20I)%20V_k%5E%5Ctop%24%2C%0A%20%20%20%20%20%20%20which%20is%20exactly%20what%20%60FactorCovariance%60%20solves%20via%20the%20Woodbury%0A%20%20%20%20%20%20%20identity%20in%20%24O(nk)%24%20memory.%0A%20%20%20%204.%20Hand%20it%20to%20%60CLA%60%20and%20plot%20the%20frontier.%0A%20%20%20%205.%20Check%20the%20operator%20against%20the%20data%3A%20the%20volatility%20the%20Woodbury%0A%20%20%20%20%20%20%20quadratic%20form%20predicts%20for%20a%20frontier%20portfolio%20has%20to%20equal%20the%0A%20%20%20%20%20%20%20volatility%20its%20realised%20return%20series%20actually%20shows.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20n_slider%20%3D%20mo.ui.slider(100%2C%201000%2C%20step%3D100%2C%20value%3D500%2C%20label%3D%22Number%20of%20assets%22)%0A%20%20%20%20n_slider%0A%20%20%20%20return%20(n_slider%2C)%0A%0A%0A%40app.function(hide_code%3DTrue)%0Adef%20simulate_returns(rng%2C%20t%2C%20n%2C%20k_true%3D10)%3A%0A%20%20%20%20%22%22%22Simulate%20t%20observations%20of%20n%20asset%20returns%20with%20k_true%20latent%20factors.%22%22%22%0A%20%20%20%20exposures%20%3D%20rng.standard_normal((n%2C%20k_true))%20*%200.3%0A%20%20%20%20factor_returns%20%3D%20rng.standard_normal((t%2C%20k_true))%0A%20%20%20%20idiosyncratic%20%3D%20rng.standard_normal((t%2C%20n))%0A%20%20%20%20return%20factor_returns%20%40%20exposures.T%20%2B%20idiosyncratic%0A%0A%0A%40app.function(hide_code%3DTrue)%0Adef%20clip_covariance(returns)%3A%0A%20%20%20%20%22%22%22Clean%20the%20sample%20covariance%20by%20Marchenko-Pastur%20eigenvalue%20clipping.%0A%0A%20%20%20%20Eigenvalues%20above%20the%20MP%20upper%20edge%20are%20kept%3B%20the%20remainder%20are%20replaced%0A%20%20%20%20by%20their%20average%2C%20preserving%20the%20trace.%20The%20result%20is%20the%0A%20%20%20%20diagonal-plus-low-rank%20model%20d%20*%20I%20%2B%20U%20%40%20diag(delta)%20%40%20U.T.%0A%20%20%20%20%22%22%22%0A%20%20%20%20t%2C%20n%20%3D%20returns.shape%0A%20%20%20%20sample%20%3D%20returns.T%20%40%20returns%20%2F%20t%0A%20%20%20%20eigenvalues%2C%20eigenvectors%20%3D%20np.linalg.eigh(sample)%0A%0A%20%20%20%20%23%20Marchenko-Pastur%20upper%20edge%20for%20variance%20sigma%5E2%20and%20aspect%20ratio%20n%2Ft%0A%20%20%20%20noise_variance%20%3D%20np.median(eigenvalues)%20%2F%20(1%20-%20np.sqrt(n%20%2F%20t))%20**%202%0A%20%20%20%20edge%20%3D%20noise_variance%20*%20(1%20%2B%20np.sqrt(n%20%2F%20t))%20**%202%0A%0A%20%20%20%20keep%20%3D%20eigenvalues%20%3E%20edge%0A%20%20%20%20d_bar%20%3D%20float(eigenvalues%5B~keep%5D.mean())%0A%0A%20%20%20%20u%20%3D%20eigenvectors%5B%3A%2C%20keep%5D%0A%20%20%20%20delta%20%3D%20eigenvalues%5Bkeep%5D%20-%20d_bar%0A%20%20%20%20return%20FactorCovariance(d%3Dnp.full(n%2C%20d_bar)%2C%20u%3Du%2C%20delta%3Ddelta)%0A%0A%0A%40app.cell%0Adef%20_(n_slider)%3A%0A%20%20%20%20rng%20%3D%20np.random.default_rng(42)%0A%20%20%20%20n%20%3D%20n_slider.value%0A%20%20%20%20returns%20%3D%20simulate_returns(rng%2C%20t%3D2%20*%20n%2C%20n%3Dn)%0A%20%20%20%20covariance%20%3D%20clip_covariance(returns)%0A%20%20%20%20mo.md(f%22Kept%20**%7Bcovariance.k%7D**%20factors%20out%20of%20%7Bn%7D%20sample%20eigenvalues.%22)%0A%20%20%20%20return%20covariance%2C%20n%2C%20returns%2C%20rng%0A%0A%0A%40app.cell%0Adef%20_(covariance%2C%20n%2C%20rng)%3A%0A%20%20%20%20frontier%20%3D%20CLA(%0A%20%20%20%20%20%20%20%20mean%3Drng.uniform(0.0%2C%200.1%2C%20n)%2C%0A%20%20%20%20%20%20%20%20covariance%3Dcovariance%2C%0A%20%20%20%20%20%20%20%20lower_bounds%3Dnp.zeros(n)%2C%0A%20%20%20%20%20%20%20%20upper_bounds%3Dnp.ones(n)%2C%0A%20%20%20%20%20%20%20%20a%3Dnp.ones((1%2C%20n))%2C%0A%20%20%20%20%20%20%20%20b%3Dnp.ones(1)%2C%0A%20%20%20%20).frontier%0A%20%20%20%20mo.md(f%22The%20exact%20frontier%20has%20**%7Blen(frontier)%7D**%20turning%20points.%22)%0A%20%20%20%20return%20(frontier%2C)%0A%0A%0A%40app.cell%0Adef%20_(frontier)%3A%0A%20%20%20%20frontier.plot(volatility%3DTrue)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Does%20the%20operator%20agree%20with%20the%20data%3F%0A%0A%20%20%20%20%60FactorCovariance%60%20never%20forms%20the%20%24n%20%5Ctimes%20n%24%20matrix%2C%20so%20the%20volatilities%0A%20%20%20%20plotted%20above%20come%20out%20of%20the%20Woodbury%20identity%20rather%20than%20a%20dense%0A%20%20%20%20quadratic%20form.%20That%20is%20worth%20checking%20against%20the%20returns%20themselves.%0A%0A%20%20%20%20Below%2C%20each%20frontier%20portfolio%20is%20applied%20to%20the%20simulated%20history%20to%20give%20a%0A%20%20%20%20realised%20return%20series%2C%20which%0A%20%20%20%20%5BjQuantStats%5D(https%3A%2F%2Fgithub.com%2Fjebel-quant%2Fjquantstats)%20measures.%20The%0A%20%20%20%20predicted%20column%20is%20%24%5Csqrt%7Bw%5E%5Ctop%20%5CSigma%20w%7D%24%20from%20the%20operator%3B%20the%0A%20%20%20%20realised%20column%20is%20the%20sample%20standard%20deviation%20of%20the%20series%2C%20and%20their%0A%20%20%20%20ratio%20sits%20within%20a%20few%20percent%20of%201%20across%20the%20whole%20slider%20range.%20Not%0A%20%20%20%20exactly%201%3A%20clipping%20deliberately%20discards%20the%20part%20of%20the%20sample%20spectrum%0A%20%20%20%20it%20calls%20noise%2C%20so%20the%20cleaned%20%24%5CSigma%24%20is%20*not*%20the%20sample%20covariance%20of%0A%20%20%20%20this%20history.%20A%20few%20percent%20is%20that%20discarded%20noise.%20An%20order-of-magnitude%0A%20%20%20%20gap%2C%20or%20one%20that%20widened%20with%20%24n%24%2C%20would%20instead%20point%20at%20the%20low-rank%0A%20%20%20%20solve.%0A%0A%20%20%20%20Risk%20is%20all%20we%20ask%20of%20this%20table.%20The%20means%20fed%20to%20the%20CLA%20are%20a%20synthetic%0A%20%20%20%20forecast%20unrelated%20to%20the%20simulation%2C%20and%20%60simulate_returns%60%20draws%20from%20a%0A%20%20%20%20standard%20normal%20rather%20than%20at%20a%20realistic%20return%20magnitude%2C%20so%20the%0A%20%20%20%20return-%2C%20Sharpe-%20and%20drawdown-based%20metrics%20jQuantStats%20also%20offers%20would%0A%20%20%20%20be%20measuring%20the%20simulation's%20conventions%20rather%20than%20the%20frontier.%0A%20%20%20%20%60experiments%2Ffrontier_stats.py%60%20reads%20those%20metrics%20off%20real%20S%26P%20500%20data%0A%20%20%20%20instead.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.function(hide_code%3DTrue)%0Adef%20score(frontier%2C%20returns)%3A%0A%20%20%20%20%22%22%22Compare%20each%20portfolio's%20predicted%20volatility%20with%20its%20realised%20one.%0A%0A%20%20%20%20Args%3A%0A%20%20%20%20%20%20%20%20frontier%20(cvxcla.types.Frontier)%3A%20The%20traced%20efficient%20frontier.%0A%20%20%20%20%20%20%20%20returns%20(numpy.ndarray)%3A%20The%20t%20x%20n%20simulated%20return%20history.%0A%0A%20%20%20%20Returns%3A%0A%20%20%20%20%20%20%20%20polars.DataFrame%3A%20One%20row%20per%20portfolio%20with%20the%20volatility%20the%0A%20%20%20%20%20%20%20%20%20%20%20%20covariance%20operator%20predicts%2C%20the%20volatility%20jQuantStats%20measures%20on%0A%20%20%20%20%20%20%20%20%20%20%20%20the%20realised%20series%2C%20and%20their%20ratio.%0A%20%20%20%20%22%22%22%0A%20%20%20%20t%2C%20n%20%3D%20returns.shape%0A%20%20%20%20%23%20Frontier%20order%20runs%20from%20maximum%20return%20towards%20minimum%20variance%2C%20but%20read%0A%20%20%20%20%23%20the%20minimiser%20off%20the%20variance%20vector%20rather%20than%20relying%20on%20that%20order.%0A%20%20%20%20min_variance%20%3D%20frontier.weights%5Bint(np.argmin(frontier.variance))%5D%0A%20%20%20%20portfolios%20%3D%20%7B%0A%20%20%20%20%20%20%20%20%22max_sharpe%22%3A%20frontier.max_sharpe%5B1%5D%2C%0A%20%20%20%20%20%20%20%20%22min_variance%22%3A%20min_variance%2C%0A%20%20%20%20%20%20%20%20%22equal_weight%22%3A%20np.full(n%2C%201.0%20%2F%20n)%2C%0A%20%20%20%20%7D%0A%0A%20%20%20%20calendar%20%3D%20pl.date_range(pl.date(2015%2C%201%2C%201)%2C%20pl.date(2045%2C%201%2C%201)%2C%20interval%3D%221d%22%2C%20eager%3DTrue)%0A%20%20%20%20dates%20%3D%20calendar.filter(calendar.dt.weekday()%20%3C%3D%205).head(t)%0A%20%20%20%20series%20%3D%20%7Bname%3A%20returns%20%40%20weights%20for%20name%2C%20weights%20in%20portfolios.items()%7D%0A%20%20%20%20data%20%3D%20Data.from_returns(pl.DataFrame(%7B%22date%22%3A%20dates%2C%20**series%7D))%0A%0A%20%20%20%20%23%20Per-observation%20volatility%2C%20so%20it%20is%20comparable%20with%20the%20frontier's%20own%3A%0A%20%20%20%20%23%20the%20simulation%20carries%20no%20annualisation%20convention.%0A%20%20%20%20realised%20%3D%20data.stats.volatility(annualize%3DFalse)%0A%20%20%20%20%23%20The%20Woodbury%20operator's%20own%20quadratic%20form%20--%20no%20dense%20matrix%20is%20formed.%0A%20%20%20%20predicted%20%3D%20%7Bname%3A%20float(np.sqrt(w%20%40%20frontier.covariance.matvec(w)))%20for%20name%2C%20w%20in%20portfolios.items()%7D%0A%20%20%20%20return%20pl.DataFrame(%0A%20%20%20%20%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%22portfolio%22%3A%20list(portfolios)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22predicted%20volatility%22%3A%20%5Bpredicted%5Bname%5D%20for%20name%20in%20portfolios%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22realised%20volatility%22%3A%20%5Brealised%5Bname%5D%20for%20name%20in%20portfolios%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22ratio%22%3A%20%5Brealised%5Bname%5D%20%2F%20predicted%5Bname%5D%20for%20name%20in%20portfolios%5D%2C%0A%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20_(frontier%2C%20returns)%3A%0A%20%20%20%20mo.ui.table(score(frontier%2C%20returns)%2C%20selection%3DNone)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A