ch5 factor model large

import%20marimo%0A%0A__generated_with%20%3D%20%220.18.3%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20!%5BMOSEK%20ApS%5D(https%3A%2F%2Fwww.mosek.com%2Fstatic%2Fimages%2Fbranding%2Fwebgraphmoseklogocolor.png%20)%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Imports%20and%20configuration%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20import%20sys%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20pandas%20as%20pd%0A%20%20%20%20import%20matplotlib%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20from%20matplotlib.colors%20import%20LinearSegmentedColormap%0A%0A%20%20%20%20from%20mosek.fusion%20import%20Model%2C%20Domain%2C%20Expr%2C%20ObjectiveSense%2C%20Matrix%2C%20Var%2C%20SolutionStatus%0A%20%20%20%20import%20mosek.fusion.pythonic%20%20%20%23%20Requires%20MOSEK%20%3E%3D%2010.2%0A%0A%20%20%20%20%23%20Options%0A%20%20%20%20np.set_printoptions(precision%3D5%2C%20linewidth%3D120%2C%20suppress%3DTrue)%0A%20%20%20%20pd.set_option(%22display.max_rows%22%2C%20None)%0A%20%20%20%20plt.rcParams%5B%22figure.figsize%22%5D%20%3D%20%5B12%2C%208%5D%0A%0A%20%20%20%20%23%20Diagnostic%0A%20%20%20%20print(f%22Python%3A%20%7Bsys.version%7D%22)%0A%20%20%20%20print(f%22marimo%3A%20%7Bmo.__version__%7D%2C%20matplotlib%3A%20%7Bmatplotlib.__version__%7D%2C%20pandas%3A%20%7Bpd.__version__%7D%2C%20numpy%3A%20%7Bnp.__version__%7D%2C%20mosek%3A%20%7BModel.getVersion()%7D%22)%0A%20%20%20%20return%20Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20SolutionStatus%2C%20np%2C%20plt%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Define%20the%20optimization%20model%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20We%20create%20a%20function%20to%20randomize%20factor%20models%2C%20i.%20e.%2C%20large%20random%20covariance%20matrices%20with%20only%20a%20few%20significant%20eigenvalues.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20random_factor_model(N%2C%20K%2C%20T)%3A%0A%20%20%20%20%20%20%20%20%23%20Generate%20K%20%2B%20N%20zero%20mean%20factors%2C%20with%20block%20covariance%3A%20%0A%20%20%20%20%20%20%20%20%23%20-%20K%20x%20K%20weighted%20diagonal%20block%20for%20the%20factors%2C%20%0A%20%20%20%20%20%20%20%20%23%20-%20N%20x%20N%20white%20noise%20(uncorrelated%20to%20the%20factors)%0A%20%20%20%20%20%20%20%20S_F%20%3D%20np.diag(np.sqrt(range(1%2C%20K%20%2B%201)))%0A%20%20%20%20%20%20%20%20Cov%20%3D%20np.block(%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%5BS_F%2C%20%20%20%20%20%20%20%20%20%20%20%20%20%20np.zeros((K%2C%20N))%5D%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Bnp.zeros((N%2C%20K))%2C%20np.eye(N)%5D%0A%20%20%20%20%20%20%20%20%5D)%0A%20%20%20%20%20%20%20%20Y%20%3D%20np.random.default_rng(seed%3D1).multivariate_normal(np.zeros(K%20%2B%20N)%2C%20Cov%2C%20T).T%0A%20%20%20%20%20%20%20%20Z_F%20%3D%20Y%5B%3AK%2C%20%3A%5D%0A%0A%20%20%20%20%20%20%20%20%23%20Generate%20random%20factor%20model%20parameters%0A%20%20%20%20%20%20%20%20B%20%3D%20np.random.default_rng(seed%3D2).normal(size%3D(N%2C%20K))%0A%20%20%20%20%20%20%20%20a%20%3D%20np.random.default_rng(seed%3D3).normal(loc%3D1%2C%20size%3D(N%2C%201))%0A%20%20%20%20%20%20%20%20e%20%3D%20Y%5BK%3A%2C%20%3A%5D%0A%0A%20%20%20%20%20%20%20%20%23%20Generate%20N%20time-series%20from%20the%20factors%0A%20%20%20%20%20%20%20%20Z%20%3D%20a%20%2B%20B%20%40%20Z_F%20%2B%20e%0A%0A%20%20%20%20%20%20%20%20%23%20Residual%20covariance%0A%20%20%20%20%20%20%20%20S_theta%20%3D%20np.cov(e)%0A%20%20%20%20%20%20%20%20diag_S_theta%20%3D%20np.diag(S_theta)%0A%0A%20%20%20%20%20%20%20%20%23%20Optimization%20parameters%0A%20%20%20%20%20%20%20%20m%20%3D%20np.mean(Z%2C%20axis%3D1)%0A%20%20%20%20%20%20%20%20S%20%3D%20np.cov(Z)%0A%20%20%20%20%20%20%20%20%23print(np.linalg.eigvalsh(np.corrcoef(Z))%5B-20%3A%5D)%0A%0A%20%20%20%20%20%20%20%20return%20m%2C%20S%2C%20B%2C%20S_F%2C%20diag_S_theta%0A%20%20%20%20return%20(random_factor_model%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Here%20we%20define%20the%20optimization%20model%20in%20MOSEK%20Fusion.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20SolutionStatus)%3A%0A%20%20%20%20%23%20Solve%20optimization%0A%20%20%20%20def%20Markowitz(N%2C%20m%2C%20G%2C%20gamma2)%3A%0A%20%20%20%20%20%20%20%20with%20Model(%22markowitz%22)%20as%20M%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Settings%0A%20%20%20%20%20%20%20%20%20%20%20%20%23M.setLogHandler(sys.stdout)%20%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Decision%20variable%20(fraction%20of%20holdings%20in%20each%20security)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20The%20variable%20x%20is%20restricted%20to%20be%20positive%2C%20which%20imposes%20the%20constraint%20of%20no%20short-selling.%20%20%20%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20M.variable(%22x%22%2C%20N%2C%20Domain.greaterThan(0.0))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Budget%20constraint%0A%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('budget'%2C%20Expr.sum(x)%20%3D%3D%201.0)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Objective%20%0A%20%20%20%20%20%20%20%20%20%20%20%20M.objective('obj'%2C%20ObjectiveSense.Maximize%2C%20x.T%20%40%20m)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Imposes%20a%20bound%20on%20the%20risk%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20isinstance(G%2C%20tuple)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20G_factor%20%3D%20G%5B0%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20g_specific%20%3D%20G%5B1%5D%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20factor_risk%20%3D%20G_factor.T%20%40%20x%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20specific_risk%20%3D%20Expr.mulElm(g_specific%2C%20x)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20total_risk%20%3D%20Expr.vstack(factor_risk%2C%20specific_risk)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('risk'%2C%20Expr.vstack(gamma2**0.5%2C%20total_risk)%2C%20Domain.inQCone())%0A%20%20%20%20%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20M.constraint('risk'%2C%20Expr.vstack(gamma2**0.5%2C%20G.T%20%40%20x)%2C%20Domain.inQCone())%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Solve%20optimization%0A%20%20%20%20%20%20%20%20%20%20%20%20M.solve()%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Check%20if%20the%20solution%20is%20an%20optimal%20point%0A%20%20%20%20%20%20%20%20%20%20%20%20solsta%20%3D%20M.getPrimalSolutionStatus()%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20(solsta%20!%3D%20SolutionStatus.Optimal)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20See%20https%3A%2F%2Fdocs.mosek.com%2Flatest%2Fpythonfusion%2Faccessing-solution.html%20about%20handling%20solution%20statuses.%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20raise%20Exception(%22Unexpected%20solution%20status!%22)%20%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20returns%20%3D%20M.primalObjValue()%0A%20%20%20%20%20%20%20%20%20%20%20%20portfolio%20%3D%20x.level()%0A%20%20%20%20%20%20%20%20%20%20%20%20time%20%3D%20M.getSolverDoubleInfo(%22optimizerTime%22)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20returns%2C%20time%0A%20%20%20%20return%20(Markowitz%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Run%20the%20optimization%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Define%20the%20parameters%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20parameters%20are%20the%20number%20of%20factors%20%24K%24%20and%20the%20risk%20limit%20%24%5Cgamma%5E2%24.%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%20%23%20Risk%20limit%0A%20%20%20%20gamma2%20%3D%200.1%0A%0A%20%20%20%20%23%20Number%20of%20factors%0A%20%20%20%20K%20%3D%2010%0A%20%20%20%20return%20K%2C%20gamma2%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Solve%20the%20optimization%20problem%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Next%20we%20generate%20random%20factor%20structured%20covariance%20matrices%20of%20different%20sizes%2C%20and%20solve%20the%20portfolio%20optimization%20both%20when%20we%20utilize%20the%20factor%20structure%20and%20when%20we%20do%20Cholesky%20factorization%20on%20it%20instead.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(K%2C%20Markowitz%2C%20gamma2%2C%20np%2C%20random_factor_model)%3A%0A%20%20%20%20%23%20Generate%20runtime%20data%20%0A%20%20%20%20%23%20NOTE%3A%20This%20can%20have%20a%20long%20runtime%2C%20depending%20on%20the%20range%20given%20for%20n%20below!%0A%20%20%20%20list_runtimes_orig%20%3D%20%5B%5D%0A%20%20%20%20list_runtimes_factor%20%3D%20%5B%5D%0A%20%20%20%20for%20n%20in%20range(5%2C%2013)%3A%0A%20%20%20%20%20%20%20%20N%20%3D%202**n%0A%20%20%20%20%20%20%20%20T%20%3D%2010000%0A%20%20%20%20%20%20%20%20m%2C%20S%2C%20B%2C%20S_F%2C%20diag_S_theta%20%3D%20random_factor_model(N%2C%20K%2C%20T)%0A%0A%20%20%20%20%20%20%20%20F%20%3D%20np.linalg.cholesky(S_F)%0A%20%20%20%20%20%20%20%20G_factor%20%3D%20B%20%40%20F%0A%20%20%20%20%20%20%20%20g_specific%20%3D%20np.sqrt(diag_S_theta)%0A%0A%20%20%20%20%20%20%20%20G_orig%20%3D%20np.linalg.cholesky(S)%0A%0A%20%20%20%20%20%20%20%20optimum_orig%2C%20runtime_orig%20%3D%20Markowitz(N%2C%20m%2C%20G_orig%2C%20gamma2)%0A%20%20%20%20%20%20%20%20optimum_factor%2C%20runtime_factor%20%3D%20Markowitz(N%2C%20m%2C%20(G_factor%2C%20g_specific)%2C%20gamma2)%0A%20%20%20%20%20%20%20%20list_runtimes_orig.append((N%2C%20runtime_orig))%0A%20%20%20%20%20%20%20%20list_runtimes_factor.append((N%2C%20runtime_factor))%0A%0A%20%20%20%20tup_N_orig%2C%20tup_time_orig%20%3D%20list(zip(*list_runtimes_orig))%0A%20%20%20%20tup_N_factor%2C%20tup_time_factor%20%3D%20list(zip(*list_runtimes_factor))%0A%20%20%20%20return%20tup_N_factor%2C%20tup_N_orig%2C%20tup_time_factor%2C%20tup_time_orig%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Plot%20results%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(plt%2C%20tup_N_factor%2C%20tup_N_orig%2C%20tup_time_factor%2C%20tup_time_orig)%3A%0A%20%20%20%20%23%20Runtime%20plot%0A%20%20%20%20plt.plot(tup_N_orig%2C%20tup_time_orig%2C%20%22-o%22)%0A%20%20%20%20plt.plot(tup_N_factor%2C%20tup_time_factor%2C%20%22-o%22)%0A%20%20%20%20plt.xlabel(%22N%22)%0A%20%20%20%20plt.ylabel(%22runtime%20(s)%22)%0A%20%20%20%20ax%20%3D%20plt.gca()%0A%20%20%20%20ax.set_xscale('log'%2C%20base%3D2)%0A%20%20%20%20ax.set_yscale('log')%0A%20%20%20%20ax.grid()%0A%20%20%20%20legend%20%3D%20%5B%22Cholesky%22%2C%20%22factor%20model%22%5D%0A%20%20%20%20plt.legend(legend)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%3Ca%20rel%3D%22license%22%20href%3D%22http%3A%2F%2Fcreativecommons.org%2Flicenses%2Fby%2F4.0%2F%22%3E%3Cimg%20alt%3D%22Creative%20Commons%20License%22%20style%3D%22border-width%3A0%22%20src%3D%22https%3A%2F%2Fi.creativecommons.org%2Fl%2Fby%2F4.0%2F80x15.png%22%20%2F%3E%3C%2Fa%3E%3Cbr%20%2F%3EThis%20work%20is%20licensed%20under%20a%20%3Ca%20rel%3D%22license%22%20href%3D%22http%3A%2F%2Fcreativecommons.org%2Flicenses%2Fby%2F4.0%2F%22%3ECreative%20Commons%20Attribution%204.0%20International%20License%3C%2Fa%3E.%20The%20**MOSEK**%20logo%20and%20name%20are%20trademarks%20of%20%3Ca%20href%3D%22http%3A%2F%2Fmosek.com%22%3EMosek%20ApS%3C%2Fa%3E.%20The%20code%20is%20provided%20as-is.%20Compatibility%20with%20future%20release%20of%20**MOSEK**%20or%20the%20%60Fusion%20API%60%20are%20not%20guaranteed.%20For%20more%20information%20contact%20our%20%5Bsupport%5D(mailto%3Asupport%40mosek.com).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20return%20(mo%2C)%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A