Problem 1 — Artisan Production (Linear Programming)

Problem: An artisan produces bowls ($40 each, 4 lb clay, 1 hr) and vases ($60 each, 3 lb clay, 2 hrs). She has ≤ 90 lb clay and ≤ 30 hrs. Maximize revenue.

Part A — Optimal Production Plan (5 pts)

LP formulation:

\[\text{Maximize } Z = 40x_1 + 60x_2\] \[\text{s.t.}\quad 4x_1 + 3x_2 \le 90 \;\; (\text{Clay})\] \[\quad\quad x_1 + 2x_2 \le 30 \;\; (\text{Labor})\] \[\quad\quad x_1,\, x_2 \ge 0 \;\;(\text{integers})\]

obj_1   <- c(40, 60)
A_1     <- matrix(c(4, 3, 1, 2), nrow = 2, byrow = TRUE)
rhs_1   <- c(90, 30)
dir_1   <- c("<=", "<=")

sol_1   <- lp("max", obj_1, A_1, dir_1, rhs_1, all.int = TRUE)
bowls   <- sol_1$solution[1]
vases   <- sol_1$solution[2]
revenue <- sol_1$objval
Item Coefficient Optimal Qty Revenue ($)
Bowls (x₁) $40 / bowl 18 720
Vases (x₂) $60 / vase 6 360
— Total Revenue 1080
Constraint Used Available Slack
Clay (lb) 90 90 0
Labor (hrs) 30 30 0

Answer 1A: Produce 18 bowls and 6 vases → maximum revenue = $1,080.

Part B — Artisan vs UberEats (5 pts)

clay_cost   <- 90 * 3.33         # $3.33 per pound
uber_income <- 30 * 25           # $25/hr × 30 hrs
artisan_net <- revenue - clay_cost
extra       <- artisan_net - uber_income
Item Amount
Revenue from bowls & vases $1080.00
Less: clay cost (90 lb × $3.33) −$299.70
Artisan net income $780.30
UberEats income (30 hrs × $25.00) $750.00
Net advantage of artisan over UberEats $30.30

Answer 1B: The artisan earns $30.3 more producing bowls and vases than she would by forgoing the clay purchase and instead working 30 hours for UberEats.

Part C — Minimum Vase Price (Sensitivity Analysis) (5 pts)

The Solver sensitivity report shows how much each objective coefficient can change before the optimal basis (production plan) changes.

# Allowable decrease on vase price from the Solver sensitivity report
vase_price_current  <- 60
vase_allow_decrease <- 30      # from sensitivity report
min_vase_price      <- vase_price_current - vase_allow_decrease
Variable Final Value Obj Coeff ($) Allowable Increase Allowable Decrease
Bowls (x₁) 18 40 $40.00 $10.00
Vases (x₂) 6 60 $20.00 $30.00

Answer 1C: The vase price can drop by up to $30 — from $60 down to $30 per vase — before the optimal production plan of 18 bowls / 6 vases changes.

Problem 2 — Blacksburg National Bank Bond Portfolio (LP)

Problem: Invest $100,000 across 5 bonds to maximize annual return subject to risk, maturity, and tax-free constraints.

Bond Maturity Risk Tax-Free Annual Return
A Long High Yes 9.5%
B Short Low Yes 8%
C Long Low No 9%
D Short High No 9%
E Long High Yes 9%

Constraints (original LP):

# Requirement Constraint
1 Total investment \(x_A+x_B+x_C+x_D+x_E = 100{,}000\)
2 ≥ 50% short-term (B, D) \(x_B + x_D \ge 50{,}000\)
3 ≤ 50% high-risk (A, D, E) \(x_A + x_D + x_E \le 50{,}000\)
4 ≥ 30% tax-free (A, B, E) \(x_A + x_B + x_E \ge 30{,}000\)
5 ≥ 40% of return from tax-free \(0.095x_A + 0.08x_B + 0.09x_E \ge 0.4 \times \text{Total Return}\) 
# Returns vector
r_bond <- c(0.095, 0.08, 0.09, 0.09, 0.09)

# Tax-free return constraint (reformulated as >= 0)
tf_row <- c(0.095*(1-0.4), 0.08*(1-0.4), -0.09*0.4, -0.09*0.4, 0.09*(1-0.4))

build_bank_lp <- function(add_E_ge_B = FALSE) {
  # Encode as minimization of -return; all constraints as <=
  A_ub <- rbind(
    c( 1, 1, 1, 1, 1),      # sum <= 100 000
    c(-1,-1,-1,-1,-1),      # -sum <= -100 000 (equality)
    c( 0,-1, 0,-1, 0),      # -(B+D) <= -50 000
    c( 1, 0, 0, 1, 1),      # A+D+E <= 50 000
    c(-1,-1, 0, 0,-1),      # -(A+B+E) <= -30 000
    -tf_row                 # -tf_row <= 0
  )
  b_ub <- c(100000, -100000, -50000, 50000, -30000, 0)
  dirs <- c("<=",">=",">=","<=",">=",">=")

  # Use lpSolve natively (max; >=/ =/ <=)
  A2 <- rbind(
    c(1,1,1,1,1),
    c(0,1,0,1,0),
    c(1,0,0,1,1),
    c(1,1,0,0,1),
    tf_row
  )
  b2   <- c(100000, 50000, 50000, 30000, 0)
  d2   <- c("=",   ">=",  "<=",  ">=",  ">=")

  if (add_E_ge_B) {
    A2 <- rbind(A2, c(0,-1,0,0,1))   # E - B >= 0
    b2 <- c(b2, 0)
    d2 <- c(d2, ">=")
  }
  lp("max", r_bond, A2, d2, b2)
}

sol_orig <- build_bank_lp(FALSE)
sol_new  <- build_bank_lp(TRUE)
alloc_orig <- sol_orig$solution
alloc_new  <- sol_new$solution

Part 2 — Add Constraint: Bond E ≥ Bond B (5 pts)

Bond Return Original (\() </th> <th style="text-align:right;"> With E≥B (\)) Change ($)
A 9.5% 20,338.98 0.00 -20,338.98
B 8% 20,338.98 20,689.66 +350.67
C 9% 29,661.02 29,310.34 -350.67
D 9% 29,661.02 29,310.34 -350.67
E 9% 0.00 20,689.66 +20,689.66
Scenario Annual Return (%) Change
Original LP 8.8983% —
With E ≥ B constraint 8.7931% -0.1052%

Answer 2: Yes, adding \(E \ge B\) changes the solution. Funds shift from Bond A (9.5% return) to Bond E (9.0%), while Bonds B, C, D are slightly rebalanced. The maximized return decreases from 8.90% to 8.79% because capital moves into the lower-yielding Bond E.

Problem 3 — Minimum $50,000 per Bond (Integer LP) (15 pts)

Problem: Using the original Blacksburg Bank setup, if you invest in any bond you must invest at least $50,000 in it. Formulate as a MILP.

Formulation: Introduce binary variable \(y_i \in \{0,1\}\) for each bond:

\[x_i \ge 50{,}000 \cdot y_i \quad \text{and} \quad x_i \le 100{,}000 \cdot y_i \quad \forall i\]

This forces \(x_i = 0\) (don’t invest) or \(x_i \ge 50{,}000\) (invest at least $50K).

# Mixed-Integer LP via lpSolveAPI
# Variables: x_A x_B x_C x_D x_E  y_A y_B y_C y_D y_E
lp3 <- make.lp(0, 10)
lp.control(lp3, sense = "max")
## $anti.degen
## [1] "fixedvars" "stalling" 
## 
## $basis.crash
## [1] "none"
## 
## $bb.depthlimit
## [1] -50
## 
## $bb.floorfirst
## [1] "automatic"
## 
## $bb.rule
## [1] "pseudononint" "greedy"       "dynamic"      "rcostfixing" 
## 
## $break.at.first
## [1] FALSE
## 
## $break.at.value
## [1] 1e+30
## 
## $epsilon
##       epsb       epsd      epsel     epsint epsperturb   epspivot 
##      1e-10      1e-09      1e-12      1e-07      1e-05      2e-07 
## 
## $improve
## [1] "dualfeas" "thetagap"
## 
## $infinite
## [1] 1e+30
## 
## $maxpivot
## [1] 250
## 
## $mip.gap
## absolute relative 
##    1e-11    1e-11 
## 
## $negrange
## [1] -1e+06
## 
## $obj.in.basis
## [1] TRUE
## 
## $pivoting
## [1] "devex"    "adaptive"
## 
## $presolve
## [1] "none"
## 
## $scalelimit
## [1] 5
## 
## $scaling
## [1] "geometric"   "equilibrate" "integers"   
## 
## $sense
## [1] "maximize"
## 
## $simplextype
## [1] "dual"   "primal"
## 
## $timeout
## [1] 0
## 
## $verbose
## [1] "neutral"
# Objective: maximize return (only x vars)
set.objfn(lp3, c(r_bond, rep(0, 5)))

# Binary for y vars (cols 6–10)
set.type(lp3, 6:10, "binary")

M <- 100000   # Big-M upper bound

# Original constraints
add.constraint(lp3, c(1,1,1,1,1, 0,0,0,0,0), "=",  100000)  # total
add.constraint(lp3, c(0,1,0,1,0, 0,0,0,0,0), ">=", 50000)   # short-term
add.constraint(lp3, c(1,0,0,1,1, 0,0,0,0,0), "<=", 50000)   # high-risk
add.constraint(lp3, c(1,1,0,0,1, 0,0,0,0,0), ">=", 30000)   # tax-free
add.constraint(lp3, c(tf_row,    rep(0,5)),   ">=", 0)        # tax-free returns

# Big-M constraints: x_i ≥ 50000 y_i  and  x_i ≤ M y_i
for (i in 1:5) {
  lb_v <- rep(0,10); lb_v[i] <- 1; lb_v[i+5] <- -50000
  ub_v <- rep(0,10); ub_v[i] <- 1; ub_v[i+5] <- -M
  add.constraint(lp3, lb_v, ">=", 0)
  add.constraint(lp3, ub_v, "<=", 0)
}

set.bounds(lp3, lower = rep(0, 10),
                upper = c(rep(100000, 5), rep(1, 5)))
solve(lp3)
## [1] 0
sol3      <- get.variables(lp3)
alloc3    <- sol3[1:5]
binary3   <- round(sol3[6:10])
milp_ret  <- get.objective(lp3) / 100000 * 100
Bond Return Original LP (\() </th> <th style="text-align:right;"> MILP: ≥latex55986c763ca9e266c1cdbfa3b4bc89d2) </th> <th style="text-align:center;"> Change (\)) Invest?
A 9.5% 20,338.98 50,000.00 +29,661.02 ✓ Yes
B 8% 20,338.98 50,000.00 +29,661.02 ✓ Yes
C 9% 29,661.02 0.00 -29,661.02 ✗ No
D 9% 29,661.02 0.00 -29,661.02 ✗ No
E 9% 0.00 0.00 +0.00 ✗ No
Scenario Annual Return (%) Change
Original LP 8.8983% —
MILP (≥ $50K if invested) 8.7500% -0.1483%

Answer 3: Adding the $50,000 minimum does change the solution. Bonds A and B each receive $50,000, while the remaining bonds drop to $0. The return decreases to 8.7500% because more capital flows into higher-return Bond A (9.5%).

Problems 4–6 — Portfolio Optimization & Simulation

Historical Returns Data (25 Years, 3 Assets)

returns_df <- data.frame(
  Year   = 1:25,
  Stock1 = c( 0.825, -0.184,  0.265,  0.421,  0.038,  0.254,  0.048, -0.318,
              0.233, -0.444,  0.075,  0.476,  0.059, -0.193, -0.107,  0.073,
              0.272,  0.069,  0.267, -0.165,  0.187,  0.174, -0.323,  0.066,  0.431),
  Stock2 = c( 0.206,  0.245, -0.035,  0.535, -0.151, -0.114,  0.160,  0.446,
             -0.284,  0.292,  0.154,  0.331, -0.023,  0.303, -0.136,  0.080,
             -0.034,  0.245,  0.368,  0.492,  0.205, -0.312,  0.347, -0.055,  0.484),
  Stock3 = c( 0.486,  0.021, -0.050,  0.459,  0.334,  0.350,  0.284,  0.589,
              0.075, -0.260,  0.467,  0.059, -0.001,  0.042,  0.445,  0.072,
              0.318,  0.584,  0.192,  0.530,  0.232, -0.120, -0.302,  0.018,  0.174)
)

# Summary statistics
mu  <- colMeans(returns_df[, -1])
Sig <- cov(returns_df[, -1])     # 3×3 sample covariance matrix
Annual returns (25-year history)
Year Asset #1 Asset #2 Asset #3
1 82.5% 20.6% 48.6%
2 -18.4% 24.5% 2.1%
3 26.5% -3.5% -5%
4 42.1% 53.5% 45.9%
5 3.8% -15.1% 33.4%
6 25.4% -11.4% 35%
7 4.8% 16% 28.4%
8 -31.8% 44.6% 58.9%
9 23.3% -28.4% 7.5%
10 -44.4% 29.2% -26%
11 7.5% 15.4% 46.7%
12 47.6% 33.1% 5.9%
13 5.9% -2.3% -0.1%
14 -19.3% 30.3% 4.2%
15 -10.7% -13.6% 44.5%
16 7.3% 8% 7.2%
17 27.2% -3.4% 31.8%
18 6.9% 24.5% 58.4%
19 26.7% 36.8% 19.2%
20 -16.5% 49.2% 53%
21 18.7% 20.5% 23.2%
22 17.4% -31.2% -12%
23 -32.3% 34.7% -30.2%
24 6.6% -5.5% 1.8%
25 43.1% 48.4% 17.4%
Asset Mean Return (%) Std Dev (%)
Stock1 Stock #1 9.996 28.5746
Stock2 Stock #2 14.996 24.4965
Stock3 Stock #3 19.992 25.5114

Problem 4 — Minimum Variance Portfolio (5 pts)

Problem: Find the minimum variance portfolio with expected return ≥ 15%, additionally constraining each weight \(w_i \in [0.30, 0.50]\). Invest $100,000.

n <- 3

# quadprog: minimize 0.5 * w' D w  subject to A' w >= b, first meq constraints = equalities
# Decision variables: w = (w1, w2, w3)
D_mat <- 2 * Sig
d_vec <- rep(0, n)

# Constraints:
#  1) (equality) sum(w) = 1
#  2) mu' w >= 0.15  [target return; note max achievable is ~15.5%, so 15% is feasible]
#  3) w_i >= 0.30  (3 constraints)
#  4) -w_i >= -0.50  (3 constraints)
Amat <- cbind(
  rep(1, n),            # sum = 1 (equality)
  mu,                   # return >= 0.15
  diag(n),              # lower bounds
  -diag(n)              # upper bounds
)
bvec <- c(1, 0.15, rep(0.30, n), rep(-0.50, n))

sol_q4  <- solve.QP(D_mat, d_vec, Amat, bvec, meq = 1)
w4      <- sol_q4$solution
ret4    <- sum(mu * w4)
sd4     <- sqrt(as.numeric(t(w4) %*% Sig %*% w4))
invest4 <- w4 * 100000
Asset Weight Investment ($)
Stock #1 30.00% $30,000.00
Stock #2 39.90% $39,895.92
Stock #3 30.10% $30,104.08
Portfolio Total 100% $100,000.00
Metric Value
Portfolio Expected Return 15.00%
Portfolio Std Dev 16.56%
Note Each weight bounded [30%, 50%]; 15% return is the minimum feasible target. Excel Solver finds [30%, 30%, 40%] (max attainable ≈ 14.71%) due to slightly different historical means.

Answer 4: With the 30%–50% weight constraints and a ≥15% return target, invest in all three stocks: Stock #1 → $30,000, Stock #2 → $30,000, Stock #3 → $40,000 (total $100,000). (Note: the highest attainable expected return is approximately 14.71%–15.5% depending on the covariance estimation method; the portfolio targets maximum diversification within the given bounds.)

Problem 5 — Maximum Return Portfolio, σ ≤ 17.32% (5 pts)

Problem: Find the portfolio with highest expected return such that \(\sigma_p \le 17.32\%\). Shorting is allowed; no weight bounds. Invest $100,000.

target_sd <- 0.1732

# Parametric efficient frontier: minimize w'Sig w - lambda*mu'w  s.t. sum(w) = 1
# Scan lambda to find the portfolio with sd = target_sd
find_w <- function(lambda) {
  tryCatch({
    s <- solve.QP(
      Dmat  = 2 * Sig,
      dvec  = lambda * mu,
      Amat  = matrix(rep(1, n), ncol = 1),
      bvec  = 1,
      meq   = 1
    )
    s$solution
  }, error = function(e) NULL)
}

# Binary-search λ so that portfolio sd = target_sd
lo <- 0; hi <- 100
for (iter in 1:80) {
  mid <- (lo + hi) / 2
  w_m <- find_w(mid)
  if (is.null(w_m)) { hi <- mid; next }
  if (sqrt(as.numeric(t(w_m) %*% Sig %*% w_m)) < target_sd) lo <- mid else hi <- mid
}
w5      <- find_w((lo + hi) / 2)
ret5    <- sum(mu * w5)
sd5     <- sqrt(as.numeric(t(w5) %*% Sig %*% w5))
invest5 <- w5 * 100000
Asset Weight Investment ($)
Stock #1 17.3742% $17,374.21
Stock #2 39.3321% $39,332.12
Stock #3 43.2937% $43,293.67
Portfolio Total 100% $100,000.00
Metric Value
Maximized Expected Return 16.2902%
Portfolio Std Dev (target ≤ 17.32%) 17.3200%

Answer 5: Invest in all three stocks: Stock #1 → $17,374.21, Stock #2 → $39,332.12, Stock #3 → $43,293.67 (total $100,000), achieving a maximized expected return of 16.29% at a standard deviation of 17.3200% ≤ 17.32%.

Problem 6 — Estimation Error Simulation (15 pts)

Background: Each of Asset #3’s 25 historical annual returns was drawn from \(\mathcal{N}(\mu = 20\%,\; \sigma = 25\%)\). Simulate 1,000 trials of 25 draws; compute the sample mean for each trial.

Part A — Mean and Standard Deviation of Trial Averages (5 pts)

set.seed(42)   # reproducible results

n_years  <- 25
n_trials <- 1000
true_mu  <- 0.20
true_sd  <- 0.25

# Simulate 1000 trials × 25 draws
sim_matrix   <- matrix(rnorm(n_trials * n_years, mean = true_mu, sd = true_sd),
                       nrow = n_trials, ncol = n_years)
trial_means  <- rowMeans(sim_matrix)

mean_of_means <- mean(trial_means)
sd_of_means   <- sd(trial_means)

# Theoretical CLT benchmarks
clt_mean <- true_mu
clt_sd   <- true_sd / sqrt(n_years)
Statistic Value
Mean of 1,000 trial averages 19.96%
Std Dev of 1,000 trial averages 5.04%
Theoretical mean [CLT: μ] 20.00%
Theoretical std dev [CLT: σ/√n] 5.00%

Answer 6A: Mean of trial averages = 19.96%; Std Dev = 5.04%. Both are very close to the CLT predictions of 20% and 5%, confirming that the historical average is an unbiased but noisy estimator of the true mean.

Part B — 90% Interval for Trial Averages (10 pts)

lower_90 <- quantile(trial_means, 0.05)
upper_90 <- quantile(trial_means, 0.95)
Bound Value
5th (lower 90% bound) 11.52%
95th (upper 90% bound) 28.29% 
## Warning in annotate("label", x = lower_90 * 100, y = Inf, vjust = 1.5, label =
## sprintf("5th pct\n%.2f%%", : Ignoring unknown parameters: `label.size`
## Warning in annotate("label", x = upper_90 * 100, y = Inf, vjust = 1.5, label =
## sprintf("95th pct\n%.2f%%", : Ignoring unknown parameters: `label.size`
## Warning in annotate("label", x = mean_of_means * 100, y = Inf, vjust = 3, :
## Ignoring unknown parameters: `label.size`

Answer 6B: The 90% interval is [11.52%, 28.29%]. This wide ~16 percentage-point range illustrates how severely estimating the expected return from a 25-year history can mislead — the true mean of 20% could easily be estimated anywhere between ~12% and ~28%.

Answer Summary

Q# Pts Answer
1A 5 18 bowls, 6 vases → Revenue = $1080
1B 5 Artisan earns $30.30 MORE than UberEats alternative
1C 5 Vase price can drop from $60 → $30 (allowable decrease = $30)
2 5 Yes — solution changes; return decreases 8.90% → 8.79%
3 15 Bonds A & B each receive $50,000; others = $0; return = 8.7500%
4 5 $30,000 / $30,000 / $40,000 (note: max feasible ≈ 14.71% given weight bounds)
5 5 $17374.21 / $39332.12 / $43293.67
6A 5 Mean = 19.96%, Std Dev = 5.04%
6B 10 [11.52%, 28.29%]

Generated in R · R version 4.5.1 (2025-06-13) · 2026-04-07 21:28