library(ROI); library(ROI.plugin.glpk)
## Warning: package 'ROI' was built under R version 4.4.3
## ROI: R Optimization Infrastructure
## Registered solver plugins: nlminb, glpk.
## Default solver: auto.
## Warning: package 'ROI.plugin.glpk' was built under R version 4.4.3
cvar_shortfall <- function(R, W0=1, T=1.0, alpha=0.95, rho=1.0){
  S <- nrow(R); D <- ncol(R)
  # Vars: x(1..D)>=0, zeta (VaR proxy), u(1..S)>=0
  nvar <- D + 1 + S
  # Objective: minimize rho * [ zeta + (1/((1-alpha)S)) * sum u_s ]
  obj <- c(rep(0, D), rho, rep(rho/((1-alpha)*S), S))
  
  # (1) sum x = 1
  A1 <- matrix(0, 1, nvar); A1[1, 1:D] <- 1
  dir <- "=="; rhs <- 1
  
  # (2) u_s >= T - W_s - zeta
  # => -u_s <= W_s + zeta - T
  # W_s = W0*(1 + R_s x)
  A2 <- matrix(0, S, nvar)
  A2[, 1:D] <-  W0 * R         # x part (moved to LHS)
  A2[, D+1] <-  1              # zeta
  A2[, (D+2):(D+1+S)] <- -diag(S)  # -u_s
  b2 <- rep(T - W0, S)
  
  L <- rbind(A1, A2)
  dir <- c(dir, rep(">=", S))  # (A2 * vars) >= b2  is u >= T - W - zeta moved correctly
  rhs <- c(rhs, b2)
  
  lb <- c(rep(0, D), -Inf, rep(0, S))
  ub <- rep(Inf, nvar)
  
  model <- OP(objective = obj,
              constraints = L_constraint(L, dir, rhs),
              bounds = V_bound(li=1:nvar, lb=lb, ui=1:nvar, ub=ub),
              maximum = FALSE)
  sol <- ROI_solve(model, solver="glpk")
  list(x = sol$solution[1:D],
       zeta = sol$solution[D+1],
       u = sol$solution[(D+2):(D+1+S)],
       objective = sol$objval)
}

# Example:
set.seed(2)
S <- 5000; D <- 3
R <- cbind(rnorm(S,0.08,0.15), rnorm(S,0.03,0.005), rnorm(S,0.02,0.02))
res <- cvar_shortfall(R, W0=1, T=1.02, alpha=0.975, rho=1.0)
print(res$x); print(res$objective)
## [1] 0.01628024 0.98371976 0.00000000
## [1] 0.008879157

1) What is the “3 × #” matrix?

\[ [ R ;=; \begin{bmatrix} r_{E,1} & r_{C,1} & r_{B,1} \ r_{E,2} & r_{C,2} & r_{B,2} \ \vdots & \vdots & \vdots \ r_{E,S} & r_{C,S} & r_{B,S} \end{bmatrix} ] \]

Each row is one scenario; each column is one asset’s return in that scenario. It is not “3×3” unless you literally have S = 3 scenarios.

The “3×3” you sometimes see is the covariance matrix \((\Sigma)\) of the three assets—that’s a different object (assets×assets), used if you generate scenarios parametrically.

2) How do I “ascertain” the entries of (R)? (Three practical ways)

A) Historical bootstrap (fast + keeps real dependence)

  1. Collect aligned historical returns for Shares, Cash, Bonds at your horizon (e.g., monthly).
  2. Sample rows with replacement to create (S) scenarios.
  3. That sampled table is your (R).

Pros: zero modelling, realistic co-moves. Cons: only repeats the past.

B) Parametric simulation (needs \((\mu)\) and \((\Sigma)\))

  1. Estimate mean vector \((\mu \in \mathbb{R}^3)\) and covariance \((\Sigma \in \mathbb{R}^{3\times3})\) from history (or judgment).

  2. Draw (S) scenarios from a joint distribution:

    • Multivariate Normal: \((R_s \sim \mathcal{N}(\mu,\Sigma))\), or
    • Multivariate Student-t (fatter tails): \((R_s \sim t_\nu(\mu,\Sigma))\).

  3. Stack the draws: that’s your (R).

Pros: easy to stress (change \((\mu,\Sigma)\), add fat tails). Cons: normal can understate crashes unless you use t or add jumps.

This is where the “3×3” shows up: \((\Sigma)\) is 3×3 because you have 3 assets.

C) Factor model (economically interpretable stress)

  1. Choose factors (F) (e.g., Market (MKT), Term (TERM), Credit (CRD)).
  2. Estimate betas (B) (assets×factors) and idiosyncratic stdevs \((\sigma_\varepsilon)\). \[ [ R ;=; \alpha \mathbf{1}^\top + B F^\top + \varepsilon ] \]
  3. Create factor scenarios (e.g., “+300 bps rates”, “credit spread blow-out”), sample \((\varepsilon)\), then compute asset returns.

Pros: lets you craft named stress events; bonds respond to TERM/CRD directly. Cons: you must pick factors and betas.

3) How to set each asset’s returns sensibly

Shares (E)

  • Historical arithmetic return at your horizon (e.g., monthly) or a forward view (equity risk premium + risk-free).
  • Stress scenarios: deep drawdowns (−20%, −35%), volatility spikes.

Cash (C)

  • Near-risk-free for the period. Monthly cash ≈ current overnight/cash rate / 12 (or the prevailing “cash management account” yield).
  • Stress: usually stable; you can nudge it for policy shocks.

Bonds / Term Deposit (B)

  • Term deposit: near-fixed return over the horizon (quoted TD rate for the term).

  • Bond fund (market-to-market): use duration/convexity with a yield shock:

    \[ [ r_B \approx\underbrace{\text{carry}}_{\text{coupon/price yield}} ;-; D ,\Delta y ;+; \tfrac{1}{2}C,(\Delta y)^2 ;-; \Delta \text{spread} ] \]

    where (D) is duration, (C) convexity, (y) the rate move. Stress scenarios: sharp rate up (neg return), spread widening (credit).

If minimum lots make bonds infeasible for small investors, set the bond column to 0 or omit it—your solver will cope.

4) Where does “average wealth” enter?

\[ [ W_s = W_0\left(1 + r_{E,s}x_E + r_{C,s}x_C + r_{B,s}x_B\right),\quad \bar W = \frac{1}{S}\sum_s W_s. ] \]

5) Sanity-check workflow (simple and reliable)

  1. Pick horizon (monthly is common).

  2. Build (R) (S×3) with one of A/B/C:

    • A: bootstrap past monthly triplets.
    • B: draw from \((\mathcal{N}(\mu,\Sigma)\)) or \((t_\nu)\).
    • C: script factor shocks; map to assets via betas.

  3. Inspect columns (means, stdevs) and correlations—do they look reasonable?

  4. Add stress rows explicitly (e.g., “GFC-like”, “rate +300 bps”), even if they’re extreme outliers.

  5. Run your chosen stress solver \((maximin/minimax/CVaR)\).

6) One-screen pseudo-code (Python) to make (R) two ways

python
import numpy as np

S = 10000  # scenarios
# Historical monthly means (illustrative):
mu = np.array([0.008, 0.0003, 0.003])           # shares, cash, bond
sd = np.array([0.045, 0.0008, 0.015])
rho = np.array([[1.00, 0.10, -0.20],            # assumed corr
                [0.10, 1.00,  0.10],
                [-0.20,0.10,  1.00]])
Sigma = np.outer(sd, sd) * rho

# (B) Parametric scenarios
R_param = np.random.multivariate_normal(mu, Sigma, size=S)

# (A) Bootstrap scenarios (if you have real columns E, C, B as arrays)
# idx = np.random.randint(0, len(E), size=S)
# R_boot = np.column_stack([E[idx], C[idx], B[idx]])

# Add explicit stress rows
stress = np.array([
    [-0.25,  0.0003, -0.05],  # equity crash + bond selloff
    [-0.35,  0.0003, +0.02],  # worse equity, bond rally
])
R = np.vstack([R_param, stress])

That’s it: R is S×3. Feed it to any of your stress solvers.

TL;DR

  • R is scenarios×assets (S×3), not 3×3; the 3×3 you’re thinking of is the covariance used when generating scenarios.
  • Average wealth is the mean of scenario wealths; not needed for maximin/minimax/CVaR, which look at worst or tail.
  • Populate (R) by bootstrap, parametric \(((\mu,\Sigma))\), or factor mapping—and feel free to append hand-crafted stress rows.

Great — here’s what that result means:

Weights \((x)\):

Objective value (z) for minimax shortfall:

### How to read it (Minimax Shortfall vs floor (T))

Recall the model: minimize the worst shortfall \((z)\) subject to \((z \ge T - W_s,\ \forall s)\), with \((W_s = 1 + r_s^\top x)\) (assuming \((W_0=1)\)).

Why the weights look like this

  • The optimizer stays almost entirely in cash, with a tiny equity tilt (~1.6%) to help reduce the largest gap to the floor without worsening the worst case too much.
  • Bonds got 0% likely because, under your stress set, their worst tails (or rate shocks) deteriorate the worst-case outcome more than cash does.

What would change the outcome

  • Raise (T) (a tougher floor): (z) will generally increase; weights may push even more to cash (or 0% equity) to protect the worst case.
  • Lower (T) (easier floor): you may see a bit more equity/bond risk, because it’s easier to keep the worst shortfall small.
  • Switch to CVaR shortfall: you’ll typically get more diversified weights (not all about the single worst scenario) because you’re optimizing tail-average rather than the single worst point.

If you tell me your actual (T) (floor) and confirm the three asset return generators/columns, I can translate (z) into the exact worst-case wealth and suggest a small parameter tweak (e.g., CVaR level ()) to get a more balanced allocation if that’s your goal.