Load Required Libraries

# Core data manipulation and visualization
library(tidyverse)      # Data manipulation and visualization
library(lubridate)      # Date handling
library(scales)         # Formatting for plots

# Financial data and analysis
library(quantmod)       # Download financial data
library(PerformanceAnalytics)  # Portfolio analytics
library(quadprog)       # Quadratic programming for optimization

# Tables and presentation
library(knitr)          # Tables
library(kableExtra)     # Enhanced tables
library(corrplot)       # Correlation visualization
library(rmdformats)     # Document formatting

# Set random seed for reproducibility
set.seed(42)

# Define color palette for dark theme plots
dark_palette <- c(
  "#00d4ff",  # Cyan

  "#ff6b6b",  # Coral

  "#4ecdc4",  # Teal
  "#ffe66d",  # Yellow

  "#95e1d3",  # Mint

  "#f38181",  # Salmon

  "#aa96da",  # Lavender
  "#fcbad3",  # Pink
  "#a8d8ea",  # Light blue

  "#ff9f43",  # Orange
  "#74b9ff"   # Sky blue
)

Download Real Market Data

We’ll work with a diversified set of US stocks representing different sectors, plus the S&P 500 as our market benchmark.

# Define our investment universe
# We're selecting stocks from different sectors for diversification
tickers <- c(
  "AAPL",   # Technology - Apple
  "MSFT",   # Technology - Microsoft
  "JPM",    # Financials - JPMorgan Chase
  "JNJ",    # Healthcare - Johnson & Johnson
  "XOM",    # Energy - Exxon Mobil
  "PG",     # Consumer Staples - Procter & Gamble
  "AMZN",   # Consumer Discretionary - Amazon
  "DIS",    # Communication Services - Disney
  "CAT",    # Industrials - Caterpillar
  "NEE"     # Utilities - NextEra Energy
)

# Market benchmark
market_ticker <- "SPY"  # S&P 500 ETF

# Risk-free rate proxy
rf_ticker <- "^IRX"  # 13-week Treasury Bill rate

# Download data period (5 years of monthly data)
start_date <- "2019-01-01"
end_date <- "2024-01-01"

# Function to safely download and extract adjusted prices
get_prices <- function(ticker, start, end) {
  tryCatch({
    data <- getSymbols(ticker, src = "yahoo", from = start, to = end, 
                       auto.assign = FALSE, periodicity = "monthly")
    Ad(data)  # Get adjusted closing prices
  }, error = function(e) {
    message(paste("Error downloading", ticker, ":", e$message))
    return(NULL)
  })
}

# Download all stock data
cat("Downloading stock data...\n")

## Downloading stock data...

stock_prices <- do.call(merge, lapply(tickers, get_prices, start_date, end_date))
colnames(stock_prices) <- tickers

# Download market data
cat("Downloading market data...\n")

## Downloading market data...

market_prices <- get_prices(market_ticker, start_date, end_date)
colnames(market_prices) <- "Market"

# Download risk-free rate (and convert from annual % to monthly decimal)
cat("Downloading risk-free rate data...\n")

## Downloading risk-free rate data...

rf_data <- get_prices(rf_ticker, start_date, end_date)
rf_rate <- mean(as.numeric(rf_data), na.rm = TRUE) / 100 / 12  # Monthly risk-free rate

cat(paste("Monthly risk-free rate:", round(rf_rate * 100, 4), "%\n"))

## Monthly risk-free rate: 0.1583 %

Calculate Returns

# Calculate monthly returns
stock_returns <- Return.calculate(stock_prices, method = "log")
market_returns <- Return.calculate(market_prices, method = "log")

# Remove first row (NA from return calculation)
stock_returns <- stock_returns[-1, ]
market_returns <- market_returns[-1, ]

# Merge to align dates automatically (xts handles this well)
all_returns <- merge(stock_returns, market_returns, join = "inner")
colnames(all_returns) <- c(tickers, "Market")

# Remove any rows with NA values
all_returns <- na.omit(all_returns)

# Separate back into stock and market returns
stock_returns <- all_returns[, tickers]
market_returns <- all_returns[, "Market"]

# Convert to data frames for easier manipulation
returns_df <- data.frame(
  Date = index(all_returns),
  coredata(all_returns)
)

# Display summary statistics
cat(paste0(
  "\n=== RETURN SUMMARY STATISTICS ===\n",
  "Number of observations: ", nrow(returns_df), "\n",
  "Date range: ", min(returns_df$Date), " to ", max(returns_df$Date), "\n\n"
))

## 
## === RETURN SUMMARY STATISTICS ===
## Number of observations: 59
## Date range: 2019-02-01 to 2023-12-01

# Summary statistics table
summary_stats <- returns_df %>%
  select(-Date) %>%
  summarise(across(everything(), list(
    Mean = ~mean(., na.rm = TRUE) * 12,  # Annualized
    StdDev = ~sd(., na.rm = TRUE) * sqrt(12),  # Annualized
    Min = ~min(., na.rm = TRUE),
    Max = ~max(., na.rm = TRUE)
  ))) %>%
  pivot_longer(everything(), 
               names_to = c("Asset", "Stat"), 
               names_pattern = "(.*)_(.*)") %>%
  pivot_wider(names_from = Stat, values_from = value)

kable(summary_stats, digits = 4, 
      caption = "Annualized Return Statistics (2019-2024)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Visualize Return Distributions

# Reshape for plotting
returns_long <- returns_df %>%
  pivot_longer(-Date, names_to = "Asset", values_to = "Return")

# Histogram of returns by asset
ggplot(returns_long, aes(x = Return, fill = Asset)) +
  geom_histogram(bins = 30, alpha = 0.8, color = "#1a1a2e") +
  facet_wrap(~Asset, scales = "free_y", ncol = 4) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "#ff6b6b", linewidth = 0.8) +
  scale_fill_manual(values = dark_palette) +
  labs(title = "Distribution of Monthly Returns by Asset",
       subtitle = "Red dashed line indicates zero return",
       x = "Monthly Return",
       y = "Frequency")

Commentary: The histograms show the return distributions for each asset. Notice that:

• Most returns cluster around zero (normal market behavior)
• Some assets show wider distributions (higher volatility)
• The red line at zero helps visualize whether returns skew positive or negative

The Core Issue: Small Input Changes -> Large Output Changes

Let’s demonstrate why portfolio optimization is so sensitive to expected return estimates.

# Calculate sample statistics
# Convert xts to matrix for calculations
stock_mat <- as.matrix(coredata(stock_returns))

# Expected returns (annualized)
exp_returns <- colMeans(stock_mat, na.rm = TRUE) * 12

# Variance-covariance matrix (annualized)
cov_matrix <- cov(stock_mat, use = "complete.obs") * 12

# Standard deviations
std_devs <- sqrt(diag(cov_matrix))

# Display expected returns
cat("=== ESTIMATED EXPECTED RETURNS (Annualized) ===\n\n")

## === ESTIMATED EXPECTED RETURNS (Annualized) ===

exp_return_df <- data.frame(
  Asset = names(exp_returns),
  Expected_Return = exp_returns,
  Std_Dev = std_devs
) %>%
  arrange(desc(Expected_Return))

kable(exp_return_df, digits = 4,
      col.names = c("Asset", "Expected Return", "Std Dev")) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Markowitz Mean-Variance Optimization

#' Solve for the optimal portfolio weights using quadratic programming
#' 
#' This function implements the Markowitz mean-variance optimization
#' to find weights that maximize the Sharpe ratio
#' 
#' @param exp_returns Vector of expected returns
#' @param cov_matrix Covariance matrix
#' @param rf Risk-free rate
#' @param allow_short Whether to allow short selling
#' @return Vector of optimal weights
optimize_portfolio <- function(exp_returns, cov_matrix, rf, allow_short = TRUE) {
  n <- length(exp_returns)
  
  # Excess returns
  excess_returns <- exp_returns - rf
  
  if (allow_short) {
    # Unconstrained optimization (allows short selling)
    # Maximize: w'mu - rf / sqrt(w'Sumw)
    # Equivalent to: minimize w'Sumw subject to w'mu = 1
    
    # Using quadprog: minimize 0.5 * w' D w - d' w
    # Subject to: A' w >= b
    
    Dmat <- 2 * cov_matrix
    dvec <- rep(0, n)
    
    # Constraints: sum of weights = 1
    Amat <- cbind(rep(1, n), excess_returns)
    bvec <- c(1, 0.01)  # Require positive expected excess return
    
    result <- solve.QP(Dmat, dvec, Amat, bvec, meq = 1)
    
    # Scale to sum to 1 and have positive Sharpe ratio
    weights <- result$solution
    
    # Rescale weights for target excess return
    weights <- weights / sum(weights)
    
  } else {
    # Constrained optimization (no short selling, max 100% in any asset)
    
    # Add non-negativity constraints
    Dmat <- 2 * cov_matrix
    dvec <- rep(0, n)
    
    # Constraints matrix: 
    # 1. Sum of weights = 1
    # 2. Each weight >= 0
    # 3. Each weight <= 1
    Amat <- cbind(
      rep(1, n),           # Sum = 1
      diag(n),             # w_i >= 0
      -diag(n)             # w_i <= 1 (i.e., -w_i >= -1)
    )
    bvec <- c(1, rep(0, n), rep(-1, n))
    
    result <- solve.QP(Dmat, dvec, Amat, bvec, meq = 1)
    weights <- result$solution
  }
  
  names(weights) <- names(exp_returns)
  return(weights)
}

# Calculate optimal portfolio with short selling allowed
weights_short <- optimize_portfolio(exp_returns, cov_matrix, rf_rate * 12, allow_short = TRUE)

# Calculate optimal portfolio without short selling
weights_long <- optimize_portfolio(exp_returns, cov_matrix, rf_rate * 12, allow_short = FALSE)

# Display comparison
weight_comparison <- data.frame(
  Asset = names(weights_short),
  With_Short_Selling = weights_short,
  Long_Only = weights_long
) %>%
  mutate(
    Difference = With_Short_Selling - Long_Only
  )

cat("\n=== OPTIMAL PORTFOLIO WEIGHTS COMPARISON ===\n\n")

## 
## === OPTIMAL PORTFOLIO WEIGHTS COMPARISON ===

kable(weight_comparison, digits = 4,
      col.names = c("Asset", "With Short Selling", "Long Only", "Difference")) %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(2, color = ifelse(weight_comparison$With_Short_Selling < 0, "red", "green")) %>%
  column_spec(3, color = "blue")

Commentary: Notice the dramatic differences:

• With short selling: The optimizer may suggest extreme positions (large positive and negative weights)
• Long only: Weights are constrained to be between 0% and 100%
• The “Difference” column shows how constraints fundamentally change the portfolio

Demonstrating Sensitivity to Expected Return Changes

#' Demonstrate how small changes in expected returns affect optimal weights
#' 
#' We'll perturb one stock's expected return by small amounts and observe
#' how the entire optimal portfolio changes
sensitivity_analysis <- function(exp_returns, cov_matrix, rf, 
                                  perturb_asset, perturbations) {
  
  results <- map_dfr(perturbations, function(delta) {
    # Create perturbed expected returns
    perturbed_returns <- exp_returns
    perturbed_returns[perturb_asset] <- exp_returns[perturb_asset] + delta
    
    # Optimize
    weights <- optimize_portfolio(perturbed_returns, cov_matrix, rf, allow_short = TRUE)
    
    # Calculate portfolio statistics
    port_return <- sum(weights * perturbed_returns)
    port_vol <- sqrt(t(weights) %*% cov_matrix %*% weights)
    sharpe <- (port_return - rf) / port_vol
    
    data.frame(
      Perturbation = delta,
      t(weights),
      Portfolio_Return = port_return,
      Portfolio_Vol = port_vol,
      Sharpe_Ratio = sharpe
    )
  })
  
  return(results)
}

# Perturbations: -5% to +5% change in expected return
perturbations <- seq(-0.05, 0.05, by = 0.01)

# Run sensitivity analysis on AAPL
sensitivity_results <- sensitivity_analysis(
  exp_returns, cov_matrix, rf_rate * 12,
  perturb_asset = "AAPL",
  perturbations = perturbations
)

cat("\n=== SENSITIVITY ANALYSIS: Perturbing AAPL Expected Return ===\n\nHow do optimal weights change when we adjust AAPL's expected return?\n\n")

## 
## === SENSITIVITY ANALYSIS: Perturbing AAPL Expected Return ===
## 
## How do optimal weights change when we adjust AAPL's expected return?

# Display key results
kable(sensitivity_results %>% 
        select(Perturbation, AAPL, MSFT, JPM, Portfolio_Return, Sharpe_Ratio) %>%
        filter(Perturbation %in% c(-0.05, -0.02, 0, 0.02, 0.05)),
      digits = 4,
      caption = "Portfolio Weights Under Different AAPL Expected Return Assumptions") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

# Visualize sensitivity
sensitivity_long <- sensitivity_results %>%
  select(Perturbation, all_of(tickers)) %>%
  pivot_longer(-Perturbation, names_to = "Asset", values_to = "Weight")

ggplot(sensitivity_long, aes(x = Perturbation, y = Weight, color = Asset)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "#ff6b6b", linewidth = 0.8) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "#ff6b6b", linewidth = 0.8) +
  scale_color_manual(values = dark_palette) +
  labs(
    title = "Sensitivity of Optimal Weights to AAPL Expected Return",
    subtitle = "Small changes in expected return cause large changes in optimal allocation",
    x = "Perturbation to AAPL Expected Return",
    y = "Optimal Weight",
    color = "Asset"
  ) +
  scale_x_continuous(labels = percent_format()) +
  scale_y_continuous(labels = percent_format())

Commentary: This chart is crucial for understanding the sensitivity problem:

• The x-axis shows changes to AAPL’s expected return (from -5% to +5%)
• The y-axis shows how the optimal weights change
• Notice how even a 1-2% change in expected return can cause weights to swing dramatically
• This is why practitioners are skeptical of pure mean-variance optimization

Overview of the Treynor-Black Approach

The Treynor-Black model provides a practical framework for active portfolio management:

1. Start with a passive benchmark (the market portfolio)
2. Identify mispriced securities (those with non-zero alpha)
3. Construct an active portfolio from mispriced securities
4. Combine active and passive to maximize Sharpe ratio

Step 1: Estimate Alpha and Beta for Each Security

#' Run CAPM regressions for each security
#' 
#' For each stock, we estimate: R_i - R_f = alpha_i + beta_i * (R_m - R_f) + e_i
#' 
#' alpha_i: Security's abnormal return (what we're looking for!)
#' beta_i: Sensitivity to market movements
#' sigma_e_i: Idiosyncratic (firm-specific) risk
run_capm_regressions <- function(stock_returns, market_returns, rf) {
  # Convert to numeric vectors/matrix
  stock_mat <- as.matrix(coredata(stock_returns))
  market_vec <- as.numeric(coredata(market_returns))
  
  results <- map_dfr(colnames(stock_mat), function(ticker) {
    # Get excess returns
    stock_excess <- stock_mat[, ticker] - rf
    market_excess <- market_vec - rf
    
    # Run regression
    model <- lm(stock_excess ~ market_excess)
    
    # Extract results
    data.frame(
      Asset = ticker,
      Alpha = coef(model)[1] * 12,  # Annualized
      Alpha_SE = summary(model)$coefficients[1, 2] * sqrt(12),
      Beta = coef(model)[2],
      Beta_SE = summary(model)$coefficients[2, 2],
      Sigma_e = sd(residuals(model)) * sqrt(12),  # Annualized
      R_squared = summary(model)$r.squared,
      T_stat_alpha = coef(model)[1] / summary(model)$coefficients[1, 2]
    )
  })
  
  return(results)
}

# Run CAPM regressions
capm_results <- run_capm_regressions(stock_returns, market_returns, rf_rate)

cat("\n=== CAPM REGRESSION RESULTS ===\n\nModel: R_i - R_f = alpha + beta x (R_m - R_f) + epsilon\n\n")

## 
## === CAPM REGRESSION RESULTS ===
## 
## Model: R_i - R_f = alpha + beta x (R_m - R_f) + epsilon

kable(capm_results, digits = 4,
      col.names = c("Asset", "Alpha (Ann.)", "Alpha SE", "Beta", 
                    "Beta SE", "sigma(epsilon) (Ann.)", "R^2", "t-stat (alpha)"),
      caption = "CAPM Regression Results for Each Security") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(2, color = ifelse(capm_results$Alpha > 0, "green", "red"))

Commentary: Understanding the regression results:

• Alpha (alpha): The annualized abnormal return. Positive alpha means the stock outperforms what CAPM predicts. Negative alpha means underperformance.
• Beta (beta): Market sensitivity. beta > 1 means more volatile than market; beta < 1 means less volatile.
• sigma(epsilon): Idiosyncratic volatility - the firm-specific risk that can be diversified away.
• R^2: How much of the stock’s movement is explained by market movements.
• t-stat (alpha): Statistical significance of alpha. |t| > 2 suggests significant alpha.

# Visualize alpha vs beta
ggplot(capm_results, aes(x = Beta, y = Alpha, label = Asset)) +
  geom_point(aes(size = Sigma_e, color = Alpha), alpha = 0.8) +
  geom_text(vjust = -1, hjust = 0.5, size = 3, color = "#e0e0e0") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "#ff6b6b", linewidth = 0.8) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "#00d4ff", linewidth = 0.8) +
  scale_color_gradient2(low = "#ff6b6b", mid = "#808080", high = "#4ecdc4", midpoint = 0) +
  labs(
    title = "Security Alpha vs Beta",
    subtitle = "Point size represents idiosyncratic risk (sigma_epsilon)",
    x = "Beta (Market Sensitivity)",
    y = "Alpha (Annualized Abnormal Return)",
    size = "sigma(epsilon)",
    color = "Alpha"
  )

Commentary: This scatter plot shows:

• Vertical axis (Alpha): Stocks above zero have positive abnormal returns; below zero have negative
• Horizontal axis (Beta): Stocks to the right of 1 are more aggressive; to the left are more defensive
• Point size: Larger points have higher idiosyncratic risk
• Ideally, we want stocks in the upper portion (positive alpha) with smaller points (low idiosyncratic risk)

Step 2: Calculate the Information Ratio for Each Security

The information ratio tells us how much alpha we get per unit of idiosyncratic risk:

\[IR_i = \frac{\alpha_i}{\sigma_{\varepsilon,i}}\]

# Calculate information ratios
capm_results <- capm_results %>%
  mutate(
    # Information ratio = alpha / idiosyncratic volatility
    Info_Ratio = Alpha / Sigma_e,
    
    # Alpha-to-variance ratio (used for weighting)
    Alpha_to_Var = Alpha / (Sigma_e^2),
    
    # Squared information ratio (contribution to portfolio Sharpe ratio)
    IR_Squared = Info_Ratio^2
  )

cat("\n=== INFORMATION RATIOS ===\n\nInformation Ratio = alpha / sigma(epsilon)\nThis measures abnormal return per unit of diversifiable risk\n\n")

## 
## === INFORMATION RATIOS ===
## 
## Information Ratio = alpha / sigma(epsilon)
## This measures abnormal return per unit of diversifiable risk

kable(capm_results %>% 
        select(Asset, Alpha, Sigma_e, Info_Ratio, IR_Squared) %>%
        arrange(desc(Info_Ratio)),
      digits = 4,
      col.names = c("Asset", "Alpha", "sigma(epsilon)", "Info Ratio", "IR^2"),
      caption = "Information Ratios Ranked from Best to Worst") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(4, color = ifelse(capm_results %>% 
                                   arrange(desc(Info_Ratio)) %>% 
                                   pull(Info_Ratio) > 0, "green", "red"))

Commentary: The information ratio is crucial in Treynor-Black:

• High positive IR: The stock offers good alpha relative to its idiosyncratic risk -> overweight
• Negative IR: The stock underperforms for its risk -> underweight or short
• The squared IR tells us how much each security contributes to the portfolio’s enhanced Sharpe ratio

Step 3: Construct the Active Portfolio

In the Treynor-Black model, we weight each security in the active portfolio proportionally to its alpha-to-variance ratio:

\[w_i^A = \frac{\alpha_i / \sigma_{\varepsilon,i}^2}{\sum_{j=1}^n \alpha_j / \sigma_{\varepsilon,j}^2}\]

#' Construct the Treynor-Black active portfolio
#' 
#' Weights are proportional to alpha/variance ratio
construct_active_portfolio <- function(capm_results) {
  # Calculate alpha-to-variance ratios
  alpha_var_ratios <- capm_results$Alpha / (capm_results$Sigma_e^2)
  names(alpha_var_ratios) <- capm_results$Asset
  
  # Normalize to sum to 1
  active_weights <- alpha_var_ratios / sum(alpha_var_ratios)
  
  # Calculate active portfolio characteristics
  alpha_A <- sum(active_weights * capm_results$Alpha)
  beta_A <- sum(active_weights * capm_results$Beta)
  sigma_e_A <- sqrt(sum(active_weights^2 * capm_results$Sigma_e^2))
  info_ratio_A <- alpha_A / sigma_e_A
  
  list(
    weights = active_weights,
    alpha = alpha_A,
    beta = beta_A,
    sigma_e = sigma_e_A,
    info_ratio = info_ratio_A
  )
}

# Construct active portfolio
active_portfolio <- construct_active_portfolio(capm_results)

cat("\n=== ACTIVE PORTFOLIO CONSTRUCTION ===\n\nWeights are proportional to alpha_i / sigma^2(epsilon_i)\n\n")

## 
## === ACTIVE PORTFOLIO CONSTRUCTION ===
## 
## Weights are proportional to alpha_i / sigma^2(epsilon_i)

# Display active portfolio weights
active_weights_df <- data.frame(
  Asset = names(active_portfolio$weights),
  Weight = active_portfolio$weights,
  Alpha = capm_results$Alpha,
  Alpha_to_Var = capm_results$Alpha_to_Var
) %>%
  arrange(desc(Weight))

kable(active_weights_df, digits = 4,
      col.names = c("Asset", "Weight in Active Portfolio", "Alpha", "alpha/sigma^2(epsilon)"),
      caption = "Active Portfolio Weights (Treynor-Black)") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(2, color = ifelse(active_weights_df$Weight > 0, "green", "red"))

cat(paste0(
  "\n=== ACTIVE PORTFOLIO CHARACTERISTICS ===\n\n",
  sprintf("Alpha (alpha_A):              %0.4f (%0.2f%% annualized)\n", 
          active_portfolio$alpha, active_portfolio$alpha * 100),
  sprintf("Beta (beta_A):               %0.4f\n", active_portfolio$beta),
  sprintf("Idiosyncratic Vol (sigma_epsilonA): %0.4f (%0.2f%% annualized)\n", 
          active_portfolio$sigma_e, active_portfolio$sigma_e * 100),
  sprintf("Information Ratio (IR_A): %0.4f\n", active_portfolio$info_ratio)
))

## 
## === ACTIVE PORTFOLIO CHARACTERISTICS ===
## 
## Alpha (alpha_A):              0.2919 (29.19% annualized)
## Beta (beta_A):               0.7952
## Idiosyncratic Vol (sigma_epsilonA): 0.1702 (17.02% annualized)
## Information Ratio (IR_A): 1.7148

Commentary: Understanding the active portfolio:

• Positive weights: Securities with positive alpha-to-variance ratios (we want more of these)
• Negative weights: Securities with negative alpha (we want to short these, or avoid them)
• The active portfolio alpha is the weighted average of individual alphas
• The active portfolio beta tells us how much market exposure the active portfolio has

# Visualize active portfolio weights
ggplot(active_weights_df, aes(x = reorder(Asset, Weight), y = Weight, fill = Weight > 0)) +
  geom_col(alpha = 0.9) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "#e0e0e0") +
  coord_flip() +
  scale_fill_manual(values = c("#ff6b6b", "#4ecdc4"), guide = "none") +
  labs(
    title = "Active Portfolio Weights (Treynor-Black Model)",
    subtitle = "Positive = Long position, Negative = Short position",
    x = "Asset",
    y = "Weight in Active Portfolio"
  ) +
  scale_y_continuous(labels = percent_format())

Step 4: Combine Active Portfolio with Market Index

Now we determine the optimal mix between the active portfolio (A) and the market index (M).

#' Calculate optimal weight in active portfolio
#' 
#' Formula: w_A* = w_A^0 / [1 + (1 - beta_A) * w_A^0]
#' where w_A^0 = (alpha_A / sigma^2_A) / (E[R_M] / sigma^2_M)
calculate_optimal_active_weight <- function(active_portfolio, market_returns, rf) {
  # Market statistics (annualized)
  market_vec <- as.numeric(coredata(market_returns))
  market_excess <- market_vec - rf
  market_return <- mean(market_excess) * 12
  market_vol <- sd(market_excess) * sqrt(12)
  
  # Active portfolio statistics
  alpha_A <- active_portfolio$alpha
  beta_A <- active_portfolio$beta
  sigma_A <- sqrt(beta_A^2 * market_vol^2 + active_portfolio$sigma_e^2)
  
  # Initial weight (before beta adjustment)
  w_A_0 <- (alpha_A / sigma_A^2) / (market_return / market_vol^2)
  
  # Adjusted weight (accounting for beta)
  w_A_star <- w_A_0 / (1 + (1 - beta_A) * w_A_0)
  
  # Market weight
  w_M_star <- 1 - w_A_star
  
  list(
    w_A_0 = w_A_0,
    w_A_star = w_A_star,
    w_M_star = w_M_star,
    market_return = market_return,
    market_vol = market_vol
  )
}

# Calculate optimal weights
optimal_weights <- calculate_optimal_active_weight(
  active_portfolio, 
  market_returns, 
  rf_rate
)

cat(paste0(
  "\n=== OPTIMAL PORTFOLIO ALLOCATION ===\n\n",
  "Step 1: Initial Active Weight (before beta adjustment)\n",
  sprintf("  w_A^0 = (alpha_A / sigma^2_A) / (E[R_M] / sigma^2_M) = %0.4f (%0.2f%%)\n\n", 
          optimal_weights$w_A_0, optimal_weights$w_A_0 * 100),
  "Step 2: Adjusted Active Weight (accounting for beta_A)\n",
  sprintf("  w_A* = w_A^0 / [1 + (1 - beta_A) x w_A^0] = %0.4f (%0.2f%%)\n\n", 
          optimal_weights$w_A_star, optimal_weights$w_A_star * 100),
  "Step 3: Market Weight\n",
  sprintf("  w_M* = 1 - w_A* = %0.4f (%0.2f%%)\n\n", 
          optimal_weights$w_M_star, optimal_weights$w_M_star * 100)
))

## 
## === OPTIMAL PORTFOLIO ALLOCATION ===
## 
## Step 1: Initial Active Weight (before beta adjustment)
##   w_A^0 = (alpha_A / sigma^2_A) / (E[R_M] / sigma^2_M) = 1.7921 (179.21%)
## 
## Step 2: Adjusted Active Weight (accounting for beta_A)
##   w_A* = w_A^0 / [1 + (1 - beta_A) x w_A^0] = 1.3110 (131.10%)
## 
## Step 3: Market Weight
##   w_M* = 1 - w_A* = -0.3110 (-31.10%)

# Calculate final weights in each individual security
final_weights <- active_portfolio$weights * optimal_weights$w_A_star
final_weights["Market"] <- optimal_weights$w_M_star

cat("\n=== FINAL PORTFOLIO WEIGHTS (All Securities + Market Index) ===\n\n")

## 
## === FINAL PORTFOLIO WEIGHTS (All Securities + Market Index) ===

final_weights_df <- data.frame(
  Asset = names(final_weights),
  Final_Weight = final_weights
) %>%
  arrange(desc(Final_Weight))

kable(final_weights_df, digits = 4,
      col.names = c("Asset", "Final Portfolio Weight"),
      caption = "Complete Treynor-Black Optimal Portfolio") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  column_spec(2, color = ifelse(final_weights_df$Final_Weight > 0, "green", "red"))

Commentary: The final Treynor-Black portfolio shows:

• Market index weight: This is our passive, diversified exposure
• Individual stock weights: These are our active bets based on alpha estimates
• The sum of all weights equals 100% (or close to it)

Step 5: Calculate Portfolio Performance Metrics

#' Calculate the Sharpe ratio enhancement from active management
#' 
#' Key formula: SR_P^2 = SR_M^2 + Sum(alpha_i / sigma_epsiloni)^2
calculate_sharpe_enhancement <- function(capm_results, market_returns, rf) {
  # Market Sharpe ratio
  market_vec <- as.numeric(coredata(market_returns))
  market_excess <- market_vec - rf
  market_sr <- mean(market_excess) * 12 / (sd(market_excess) * sqrt(12))
  
  # Sum of squared information ratios
  sum_ir_squared <- sum(capm_results$IR_Squared)
  
  # Portfolio Sharpe ratio
  portfolio_sr_squared <- market_sr^2 + sum_ir_squared
  portfolio_sr <- sqrt(portfolio_sr_squared)
  
  list(
    market_sr = market_sr,
    sum_ir_squared = sum_ir_squared,
    portfolio_sr = portfolio_sr,
    improvement = portfolio_sr - market_sr,
    improvement_pct = (portfolio_sr / market_sr - 1) * 100
  )
}

sharpe_enhancement <- calculate_sharpe_enhancement(capm_results, market_returns, rf_rate)

cat(paste0(
  "\n=== SHARPE RATIO ENHANCEMENT FROM ACTIVE MANAGEMENT ===\n\n",
  "Formula: SR_P^2 = SR_M^2 + Sum(alpha_i / sigma_epsiloni)^2\n\n",
  sprintf("Market Sharpe Ratio (SR_M):     %0.4f\n", sharpe_enhancement$market_sr),
  sprintf("Sum of IR^2 from active stocks:  %0.4f\n", sharpe_enhancement$sum_ir_squared),
  sprintf("Portfolio Sharpe Ratio (SR_P):  %0.4f\n", sharpe_enhancement$portfolio_sr),
  sprintf("\nSharpe Ratio Improvement:       %0.4f (%0.2f%% increase)\n", 
          sharpe_enhancement$improvement, sharpe_enhancement$improvement_pct)
))

## 
## === SHARPE RATIO ENHANCEMENT FROM ACTIVE MANAGEMENT ===
## 
## Formula: SR_P^2 = SR_M^2 + Sum(alpha_i / sigma_epsiloni)^2
## 
## Market Sharpe Ratio (SR_M):     0.5964
## Sum of IR^2 from active stocks:  2.9404
## Portfolio Sharpe Ratio (SR_P):  1.8155
## 
## Sharpe Ratio Improvement:       1.2191 (204.41% increase)

Commentary: The Sharpe ratio enhancement demonstrates the value of active management:

• The market Sharpe ratio is what you get from passive investing
• The sum of squared IRs represents the contribution from security selection
• The portfolio Sharpe ratio shows total risk-adjusted performance
• This formula guarantees that optimal active management can only improve (never hurt) the Sharpe ratio

# Visualize contribution to Sharpe ratio
sr_contribution <- capm_results %>%
  select(Asset, IR_Squared) %>%
  add_row(Asset = "Market (Base)", IR_Squared = sharpe_enhancement$market_sr^2) %>%
  mutate(
    Cumulative_SR = sqrt(cumsum(IR_Squared)),
    Contribution_Type = ifelse(Asset == "Market (Base)", "Passive", "Active")
  )

ggplot(sr_contribution, aes(x = reorder(Asset, -IR_Squared), y = IR_Squared, fill = Contribution_Type)) +
  geom_col(alpha = 0.9) +
  labs(
    title = "Contribution to Portfolio Sharpe Ratio^2",
    subtitle = "Active management adds to the base market Sharpe ratio",
    x = "Asset",
    y = "Contribution to SR^2 (Information Ratio^2)",
    fill = "Type"
  ) +
  scale_fill_manual(values = c("Active" = "#00d4ff", "Passive" = "#4ecdc4")) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Understanding Tracking Error

Tracking error measures how much your portfolio deviates from the benchmark:

\[TE = R_P - R_M\]

The tracking risk is the standard deviation of tracking error:

\[\sigma(TE) = w_A^* \sqrt{(1 - \beta_A)^2 \sigma_M^2 + \sigma_{\varepsilon A}^2}\]

#' Calculate tracking risk
#' 
#' Tracking risk measures how much the portfolio deviates from the benchmark
calculate_tracking_risk <- function(active_portfolio, market_returns, 
                                     active_weight, rf) {
  # Market volatility (annualized)
  market_vec <- as.numeric(coredata(market_returns))
  market_vol <- sd(market_vec - rf) * sqrt(12)
  
  # Components of tracking risk
  beta_A <- active_portfolio$beta
  sigma_e_A <- active_portfolio$sigma_e
  
  # Tracking risk formula
  tracking_risk <- active_weight * sqrt((1 - beta_A)^2 * market_vol^2 + sigma_e_A^2)
  
  # If we put 100% in active portfolio
  tracking_risk_full <- sqrt((1 - beta_A)^2 * market_vol^2 + sigma_e_A^2)
  
  list(
    tracking_risk = tracking_risk,
    tracking_risk_full = tracking_risk_full,
    beta_component = (1 - beta_A)^2 * market_vol^2,
    idiosyncratic_component = sigma_e_A^2
  )
}

tracking <- calculate_tracking_risk(
  active_portfolio, 
  market_returns, 
  optimal_weights$w_A_star, 
  rf_rate
)

cat(paste0(
  "\n=== TRACKING RISK ANALYSIS ===\n\n",
  "Tracking Error: TE = R_P - R_M\n",
  "Tracking Risk: sigma(TE) = w_A x sqrt[(1-beta_A)^2 x sigma_M^2 + sigma^2_epsilonA]\n\n",
  sprintf("Active Portfolio Beta (beta_A):         %0.4f\n", active_portfolio$beta),
  sprintf("Active Portfolio sigma_epsilonA:               %0.4f (%0.2f%%)\n", 
          active_portfolio$sigma_e, active_portfolio$sigma_e * 100),
  sprintf("Optimal Active Weight (w_A*):        %0.4f (%0.2f%%)\n\n", 
          optimal_weights$w_A_star, optimal_weights$w_A_star * 100),
  sprintf("Beta Contribution to Variance:       %0.6f\n", tracking$beta_component),
  sprintf("Idiosyncratic Contribution:          %0.6f\n", tracking$idiosyncratic_component),
  sprintf("\nTracking Risk (Current Allocation):  %0.4f (%0.2f%% annualized)\n", 
          tracking$tracking_risk, tracking$tracking_risk * 100),
  sprintf("Tracking Risk (100%% Active):         %0.4f (%0.2f%% annualized)\n", 
          tracking$tracking_risk_full, tracking$tracking_risk_full * 100)
))

## 
## === TRACKING RISK ANALYSIS ===
## 
## Tracking Error: TE = R_P - R_M
## Tracking Risk: sigma(TE) = w_A x sqrt[(1-beta_A)^2 x sigma_M^2 + sigma^2_epsilonA]
## 
## Active Portfolio Beta (beta_A):         0.7952
## Active Portfolio sigma_epsilonA:               0.1702 (17.02%)
## Optimal Active Weight (w_A*):        1.3110 (131.10%)
## 
## Beta Contribution to Variance:       0.001475
## Idiosyncratic Contribution:          0.028973
## 
## Tracking Risk (Current Allocation):  0.2288 (22.88% annualized)
## Tracking Risk (100% Active):         0.1745 (17.45% annualized)

Commentary: Tracking risk tells us:

• Beta component: If beta_A != 1, market movements themselves create tracking error
• Idiosyncratic component: Firm-specific risks that cause deviation from benchmark
• Current tracking risk: Based on our optimal allocation
• Full active tracking risk: What we’d have if we put 100% in the active portfolio

Setting Tracking Risk Constraints

In practice, portfolio managers often face tracking risk limits:

#' Calculate maximum active weight given tracking risk constraint
#' 
#' w_A(TE) = sigma*(TE) / sqrt[(1-beta_A)^2 x sigma_M^2 + sigma^2_epsilonA]
max_active_weight <- function(target_tracking_risk, active_portfolio, 
                               market_returns, rf) {
  market_vec <- as.numeric(coredata(market_returns))
  market_vol <- sd(market_vec - rf) * sqrt(12)
  beta_A <- active_portfolio$beta
  sigma_e_A <- active_portfolio$sigma_e
  
  denominator <- sqrt((1 - beta_A)^2 * market_vol^2 + sigma_e_A^2)
  target_tracking_risk / denominator
}

# Calculate for various tracking risk targets
tracking_targets <- seq(0.01, 0.10, by = 0.01)

tracking_analysis <- map_dfr(tracking_targets, function(target) {
  max_w_A <- max_active_weight(target, active_portfolio, market_returns, rf_rate)
  
  data.frame(
    Target_Tracking_Risk = target,
    Max_Active_Weight = min(max_w_A, 1),  # Cap at 100%
    Market_Weight = max(1 - max_w_A, 0)
  )
})

cat("\n=== TRACKING RISK CONSTRAINT ANALYSIS ===\n\nMaximum Active Weight = Target sigma(TE) / sqrt[(1-beta_A)^2 x sigma_M^2 + sigma^2_epsilonA]\n\n")

## 
## === TRACKING RISK CONSTRAINT ANALYSIS ===
## 
## Maximum Active Weight = Target sigma(TE) / sqrt[(1-beta_A)^2 x sigma_M^2 + sigma^2_epsilonA]

kable(tracking_analysis, digits = 4,
      col.names = c("Target Tracking Risk", "Max Active Weight", "Min Market Weight"),
      caption = "Maximum Active Portfolio Weight Under Different Tracking Risk Limits") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

# Visualize the trade-off
ggplot(tracking_analysis, aes(x = Target_Tracking_Risk)) +
  geom_line(aes(y = Max_Active_Weight, color = "Active Portfolio"), linewidth = 1.2) +
  geom_line(aes(y = Market_Weight, color = "Market Index"), linewidth = 1.2) +
  geom_point(aes(y = Max_Active_Weight), size = 2, color = "#00d4ff") +
  geom_point(aes(y = Market_Weight), size = 2, color = "#4ecdc4") +
  geom_vline(xintercept = tracking$tracking_risk, linetype = "dashed", color = "#ff6b6b") +
  annotate("text", x = tracking$tracking_risk + 0.005, y = 0.5, 
           label = "Current\nOptimal", color = "#ff6b6b", size = 3) +
  labs(
    title = "Portfolio Allocation Under Tracking Risk Constraints",
    subtitle = "Higher tracking risk tolerance allows more active management",
    x = "Target Tracking Risk (Annualized)",
    y = "Portfolio Weight",
    color = "Allocation"
  ) +
  scale_x_continuous(labels = percent_format()) +
  scale_y_continuous(labels = percent_format()) +
  scale_color_manual(values = c("Active Portfolio" = "#00d4ff", 
                                "Market Index" = "#4ecdc4"))

Commentary: This chart shows the fundamental trade-off:

• Tight tracking risk constraint (left): Forces more market index allocation (closet indexing)
• Loose tracking risk constraint (right): Allows more active portfolio allocation
• The red dashed line shows our current optimal allocation
• Fund managers must balance potential outperformance against career risk of significant underperformance

The Problem of Overconfident Forecasts

Analysts often overestimate the magnitude of their alpha predictions. Historical backtesting can help adjust forecasts.

#' Simulate an analyst's forecasting track record
#' 
#' We'll create synthetic data showing how forecasts compare to actual outcomes
set.seed(123)

# Simulate 40 quarters (10 years) of forecasting data
n_quarters <- 40

# True relationship: Actual = 0.005 + 0.4 x Forecast + noise
# This means the analyst's forecasts are inflated (only 40% materializes)
true_intercept <- 0.005  # Small positive bias
true_slope <- 0.4        # Forecasts are overconfident (60% shrinkage needed)
noise_sd <- 0.02         # Estimation error

# Analyst's raw forecasts (they tend to be overconfident)
analyst_forecasts <- rnorm(n_quarters, mean = 0.08, sd = 0.05)

# Actual realized abnormal returns
actual_returns <- true_intercept + true_slope * analyst_forecasts + 
                  rnorm(n_quarters, mean = 0, sd = noise_sd)

# Create data frame
forecast_data <- data.frame(
  Quarter = 1:n_quarters,
  Year = rep(2014:2023, each = 4),
  Q = rep(paste0("Q", 1:4), 10),
  Forecast = analyst_forecasts,
  Actual = actual_returns
)

cat(paste0(
  "\n=== ANALYST FORECASTING TRACK RECORD ===\n\n",
  "We've simulated 40 quarters of alpha forecasts and actual outcomes.\n",
  "Let's assess how accurate this analyst's forecasts have been.\n\n",
  "Forecast Statistics:\n",
  sprintf("  Mean Forecast:    %0.4f (%0.2f%%)\n", 
          mean(forecast_data$Forecast), mean(forecast_data$Forecast) * 100),
  sprintf("  Mean Actual:      %0.4f (%0.2f%%)\n", 
          mean(forecast_data$Actual), mean(forecast_data$Actual) * 100),
  sprintf("  Forecast Std Dev: %0.4f\n", sd(forecast_data$Forecast)),
  sprintf("  Actual Std Dev:   %0.4f\n\n", sd(forecast_data$Actual))
))

## 
## === ANALYST FORECASTING TRACK RECORD ===
## 
## We've simulated 40 quarters of alpha forecasts and actual outcomes.
## Let's assess how accurate this analyst's forecasts have been.
## 
## Forecast Statistics:
##   Mean Forecast:    0.0823 (8.23%)
##   Mean Actual:      0.0378 (3.78%)
##   Forecast Std Dev: 0.0449
##   Actual Std Dev:   0.0270

# Display recent data
kable(tail(forecast_data, 12), digits = 4,
      caption = "Last 12 Quarters of Forecasting Data") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Evaluating Forecast Accuracy

We regress actual outcomes on forecasts to assess the analyst’s skill:

\[u_t = a + b \times \alpha_t^f + \varepsilon_t\]

# Run the accuracy regression
accuracy_model <- lm(Actual ~ Forecast, data = forecast_data)

cat("\n=== FORECAST ACCURACY REGRESSION ===\n\nModel: Actual Return = a + b x Forecast + epsilon\n\n")

## 
## === FORECAST ACCURACY REGRESSION ===
## 
## Model: Actual Return = a + b x Forecast + epsilon

summary(accuracy_model)

## 
## Call:
## lm(formula = Actual ~ Forecast, data = forecast_data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.045669 -0.010855 -0.001053  0.010246  0.043485 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.003028   0.006472   0.468    0.643    
## Forecast    0.422342   0.069267   6.097 4.19e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01942 on 38 degrees of freedom
## Multiple R-squared:  0.4945, Adjusted R-squared:  0.4812 
## F-statistic: 37.18 on 1 and 38 DF,  p-value: 4.186e-07

# Extract coefficients
intercept <- coef(accuracy_model)[1]
slope <- coef(accuracy_model)[2]
r_squared <- summary(accuracy_model)$r.squared

bias_type <- ifelse(intercept > 0, "positive", "negative")

cat(paste0(
  "\n=== INTERPRETATION ===\n\n",
  sprintf("Intercept (a):  %0.4f\n", intercept),
  sprintf("  -> The analyst has a slight %s bias\n\n", bias_type),
  sprintf("Slope (b):      %0.4f\n", slope),
  sprintf("  -> Only %.0f%% of the analyst's forecast typically materializes\n", slope * 100),
  sprintf("  -> The analyst is OVERCONFIDENT by %.0f%%\n\n", (1 - slope) * 100),
  sprintf("R-squared:      %0.4f\n", r_squared),
  sprintf("  -> The forecast explains %.0f%% of variation in actual outcomes\n", r_squared * 100)
))

## 
## === INTERPRETATION ===
## 
## Intercept (a):  0.0030
##   -> The analyst has a slight positive bias
## 
## Slope (b):      0.4223
##   -> Only 42% of the analyst's forecast typically materializes
##   -> The analyst is OVERCONFIDENT by 58%
## 
## R-squared:      0.4945
##   -> The forecast explains 49% of variation in actual outcomes

# Visualize forecast vs actual
# Using base geom_smooth without message
forecast_plot <- ggplot(forecast_data, aes(x = Forecast, y = Actual)) +
  geom_point(alpha = 0.7, size = 3, color = "#00d4ff") +
  geom_abline(intercept = intercept, slope = slope, color = "#ff6b6b", linewidth = 1) +
  geom_ribbon(aes(ymin = predict(accuracy_model, interval = "confidence")[,2],
                  ymax = predict(accuracy_model, interval = "confidence")[,3]),
              alpha = 0.2, fill = "#ff6b6b") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "#4ecdc4", linewidth = 0.8) +
  annotate("text", x = 0.15, y = 0.15, label = "Perfect Forecast\n(45 degree line)", 
           color = "#4ecdc4", size = 3) +
  annotate("text", x = 0.15, y = slope * 0.15 + intercept + 0.02, 
           label = "Actual Relationship", color = "#ff6b6b", size = 3) +
  labs(
    title = "Analyst Forecast Accuracy Assessment",
    subtitle = paste0("Actual = ", round(intercept, 3), " + ", round(slope, 3), 
                      " x Forecast (R-squared = ", round(r_squared, 3), ")"),
    x = "Analyst's Alpha Forecast",
    y = "Actual Realized Abnormal Return"
  ) +
  scale_x_continuous(labels = percent_format()) +
  scale_y_continuous(labels = percent_format())

print(forecast_plot)

Commentary: This chart is crucial for understanding forecast adjustment:

• Blue dashed line (45 degree): What perfect forecasting would look like (Forecast = Actual)
• Red line: The actual relationship between forecasts and outcomes
• The gap between the lines shows the analyst’s overconfidence
• Points should scatter around the red line with some noise

Adjusting Future Forecasts

Now when the analyst issues a new forecast, we shrink it based on historical accuracy:

\[\alpha_{adjusted} = a + b \times \alpha_{forecast}\]

#' Adjust a new forecast based on historical accuracy
#' 
#' This is the BLUE (Best Linear Unbiased Estimator) of the true alpha
adjust_forecast <- function(raw_forecast, accuracy_model) {
  intercept <- coef(accuracy_model)[1]
  slope <- coef(accuracy_model)[2]
  
  adjusted <- intercept + slope * raw_forecast
  
  list(
    raw = raw_forecast,
    adjusted = adjusted,
    shrinkage = (raw_forecast - adjusted) / raw_forecast * 100
  )
}

# Example: Analyst forecasts 20% alpha for a new stock
new_forecast <- 0.20

adjusted <- adjust_forecast(new_forecast, accuracy_model)

cat(paste0(
  "\n=== ADJUSTING A NEW FORECAST ===\n\n",
  sprintf("Analyst's Raw Forecast: %0.4f (%0.2f%%)\n", 
          adjusted$raw, adjusted$raw * 100),
  sprintf("Adjusted Forecast:      %0.4f (%0.2f%%)\n", 
          adjusted$adjusted, adjusted$adjusted * 100),
  sprintf("Shrinkage Applied:      %0.1f%%\n\n", adjusted$shrinkage),
  "The adjusted forecast accounts for the analyst's historical tendency\n",
  "to overestimate alphas. This is statistically optimal (BLUE).\n"
))

## 
## === ADJUSTING A NEW FORECAST ===
## 
## Analyst's Raw Forecast: 0.2000 (20.00%)
## Adjusted Forecast:      0.0875 (8.75%)
## Shrinkage Applied:      56.3%
## 
## The adjusted forecast accounts for the analyst's historical tendency
## to overestimate alphas. This is statistically optimal (BLUE).

# Show shrinkage for various forecast levels
forecast_levels <- seq(0.05, 0.30, by = 0.05)

shrinkage_table <- map_dfr(forecast_levels, function(f) {
  adj <- adjust_forecast(f, accuracy_model)
  data.frame(
    Raw_Forecast = f,
    Adjusted_Forecast = adj$adjusted,
    Shrinkage_Pct = adj$shrinkage
  )
})

cat("\n=== SHRINKAGE TABLE ===\nHow different forecast levels are adjusted:\n\n")

## 
## === SHRINKAGE TABLE ===
## How different forecast levels are adjusted:

kable(shrinkage_table, digits = 4,
      col.names = c("Raw Forecast", "Adjusted Forecast", "Shrinkage %"),
      caption = "Forecast Adjustment at Different Levels") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Commentary: The shrinkage adjustment is important because:

• Higher forecasts are shrunk more (in absolute terms)
• The adjustment is based on the analyst’s actual track record
• Using adjusted forecasts leads to more realistic portfolio weights
• This prevents the optimizer from overreacting to overconfident predictions

Complete Implementation with Adjusted Alphas

#' Complete Treynor-Black implementation with all adjustments
#' 
#' This function implements the full model including:
#' 1. CAPM regression for alpha/beta estimation
#' 2. Alpha adjustment for forecast precision
#' 3. Active portfolio construction
#' 4. Optimal combination with market index
#' 5. Tracking risk calculation
complete_treynor_black <- function(stock_returns, market_returns, rf, 
                                    alpha_shrinkage = 0.5,  # Default 50% shrinkage
                                    tracking_risk_limit = NULL) {
  
  # Step 1: Run CAPM regressions
  capm_results <- run_capm_regressions(stock_returns, market_returns, rf)
  
  # Step 2: Apply alpha shrinkage (simulating forecast adjustment)
  capm_results$Alpha_Adjusted <- capm_results$Alpha * alpha_shrinkage
  capm_results$Info_Ratio_Adj <- capm_results$Alpha_Adjusted / capm_results$Sigma_e
  capm_results$Alpha_to_Var_Adj <- capm_results$Alpha_Adjusted / (capm_results$Sigma_e^2)
  
  # Step 3: Construct active portfolio with adjusted alphas
  alpha_var_ratios <- capm_results$Alpha_to_Var_Adj
  names(alpha_var_ratios) <- capm_results$Asset
  active_weights <- alpha_var_ratios / sum(alpha_var_ratios)
  
  alpha_A <- sum(active_weights * capm_results$Alpha_Adjusted)
  beta_A <- sum(active_weights * capm_results$Beta)
  sigma_e_A <- sqrt(sum(active_weights^2 * capm_results$Sigma_e^2))
  
  # Step 4: Calculate optimal active weight
  market_vec <- as.numeric(coredata(market_returns))
  market_excess <- market_vec - rf
  market_return <- mean(market_excess) * 12
  market_vol <- sd(market_excess) * sqrt(12)
  
  sigma_A <- sqrt(beta_A^2 * market_vol^2 + sigma_e_A^2)
  w_A_0 <- (alpha_A / sigma_A^2) / (market_return / market_vol^2)
  w_A_star <- w_A_0 / (1 + (1 - beta_A) * w_A_0)
  
  # Step 5: Apply tracking risk constraint if specified
  if (!is.null(tracking_risk_limit)) {
    tracking_denominator <- sqrt((1 - beta_A)^2 * market_vol^2 + sigma_e_A^2)
    max_w_A <- tracking_risk_limit / tracking_denominator
    w_A_star <- min(w_A_star, max_w_A)
  }
  
  w_M_star <- 1 - w_A_star
  
  # Final portfolio weights
  final_weights <- c(active_weights * w_A_star, Market = w_M_star)
  
  # Portfolio statistics
  tracking_risk <- w_A_star * sqrt((1 - beta_A)^2 * market_vol^2 + sigma_e_A^2)
  
  market_sr <- market_return / market_vol
  sum_ir_sq <- sum(capm_results$Info_Ratio_Adj^2)
  portfolio_sr <- sqrt(market_sr^2 + sum_ir_sq)
  
  list(
    capm_results = capm_results,
    active_weights = active_weights,
    active_alpha = alpha_A,
    active_beta = beta_A,
    active_sigma_e = sigma_e_A,
    w_A_star = w_A_star,
    w_M_star = w_M_star,
    final_weights = final_weights,
    tracking_risk = tracking_risk,
    market_sharpe = market_sr,
    portfolio_sharpe = portfolio_sr,
    sharpe_improvement = portfolio_sr - market_sr
  )
}

Compare Different Scenarios

# Scenario 1: No alpha adjustment, no tracking risk constraint
scenario_1 <- complete_treynor_black(
  stock_returns, market_returns, rf_rate,
  alpha_shrinkage = 1.0,  # Use raw alphas
  tracking_risk_limit = NULL
)

# Scenario 2: 50% alpha shrinkage, no tracking risk constraint  
scenario_2 <- complete_treynor_black(
  stock_returns, market_returns, rf_rate,
  alpha_shrinkage = 0.5,
  tracking_risk_limit = NULL
)

# Scenario 3: 50% alpha shrinkage with 5% tracking risk limit
scenario_3 <- complete_treynor_black(
  stock_returns, market_returns, rf_rate,
  alpha_shrinkage = 0.5,
  tracking_risk_limit = 0.05
)

# Scenario 4: 50% alpha shrinkage with 2% tracking risk limit (very constrained)
scenario_4 <- complete_treynor_black(
  stock_returns, market_returns, rf_rate,
  alpha_shrinkage = 0.5,
  tracking_risk_limit = 0.02
)

# Compare scenarios
scenario_comparison <- data.frame(
  Scenario = c("Raw Alphas, No TE Limit",
               "50% Shrinkage, No TE Limit",
               "50% Shrinkage, 5% TE Limit",
               "50% Shrinkage, 2% TE Limit"),
  Active_Weight = c(scenario_1$w_A_star, scenario_2$w_A_star, 
                    scenario_3$w_A_star, scenario_4$w_A_star),
  Market_Weight = c(scenario_1$w_M_star, scenario_2$w_M_star,
                    scenario_3$w_M_star, scenario_4$w_M_star),
  Tracking_Risk = c(scenario_1$tracking_risk, scenario_2$tracking_risk,
                    scenario_3$tracking_risk, scenario_4$tracking_risk),
  Portfolio_SR = c(scenario_1$portfolio_sharpe, scenario_2$portfolio_sharpe,
                   scenario_3$portfolio_sharpe, scenario_4$portfolio_sharpe),
  SR_Improvement = c(scenario_1$sharpe_improvement, scenario_2$sharpe_improvement,
                     scenario_3$sharpe_improvement, scenario_4$sharpe_improvement)
)

cat("\n=== SCENARIO COMPARISON ===\n\n")

## 
## === SCENARIO COMPARISON ===

kable(scenario_comparison, digits = 4,
      col.names = c("Scenario", "Active Weight", "Market Weight", 
                    "Tracking Risk", "Portfolio SR", "SR Improvement"),
      caption = "Comparison of Different Treynor-Black Implementation Scenarios") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

Commentary: The scenario analysis reveals important trade-offs:

• Raw alphas lead to more aggressive active positions (potentially unrealistic)
• Alpha shrinkage reduces active weights and tracking risk
• Tracking risk constraints force more conservative allocations
• Sharpe ratio improvement decreases as constraints tighten (trade-off between risk control and performance)

# Visualize scenario comparison
scenario_long <- scenario_comparison %>%
  pivot_longer(cols = c(Active_Weight, Market_Weight),
               names_to = "Allocation_Type",
               values_to = "Weight") %>%
  mutate(Allocation_Type = ifelse(Allocation_Type == "Active_Weight", 
                                   "Active Portfolio", "Market Index"))

ggplot(scenario_long, aes(x = Scenario, y = Weight, fill = Allocation_Type)) +
  geom_col(position = "stack", alpha = 0.9) +
  geom_text(aes(label = sprintf("%.1f%%", Weight * 100)), 
            position = position_stack(vjust = 0.5), size = 3, color = "#1a1a2e") +
  coord_flip() +
  labs(
    title = "Portfolio Allocation Under Different Scenarios",
    subtitle = "How alpha adjustment and tracking risk limits affect allocation",
    x = "",
    y = "Portfolio Weight",
    fill = "Allocation"
  ) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_manual(values = c("Active Portfolio" = "#00d4ff", 
                               "Market Index" = "#4ecdc4"))

Final Portfolio Summary

# Use Scenario 2 as our recommended portfolio
recommended <- scenario_2

# Build the output as a single string
output <- paste0(
  "\n",
  "```\n",
  "+================================================================+\n",
  "|         RECOMMENDED TREYNOR-BLACK PORTFOLIO                   |\n",
  "|         (50% Alpha Shrinkage, No Tracking Risk Limit)         |\n",
  "+================================================================+\n\n",
  "ALLOCATION SUMMARY\n",
  "------------------------------------------------------------------\n",
  sprintf("Active Portfolio Weight:    %6.2f%%\n", recommended$w_A_star * 100),
  sprintf("Market Index Weight:        %6.2f%%\n", recommended$w_M_star * 100),
  "------------------------------------------------------------------\n\n",
  "PERFORMANCE METRICS\n",
  "------------------------------------------------------------------\n",
  sprintf("Market Sharpe Ratio:        %6.4f\n", recommended$market_sharpe),
  sprintf("Portfolio Sharpe Ratio:     %6.4f\n", recommended$portfolio_sharpe),
  sprintf("Sharpe Ratio Improvement:   %6.4f\n", recommended$sharpe_improvement),
  sprintf("Tracking Risk:              %6.2f%%\n", recommended$tracking_risk * 100),
  "------------------------------------------------------------------\n\n",
  "ACTIVE PORTFOLIO CHARACTERISTICS\n",
  "------------------------------------------------------------------\n",
  sprintf("Alpha (Adjusted):           %6.4f (%0.2f%%)\n", 
          recommended$active_alpha, recommended$active_alpha * 100),
  sprintf("Beta:                       %6.4f\n", recommended$active_beta),
  sprintf("Idiosyncratic Vol:          %6.4f (%0.2f%%)\n", 
          recommended$active_sigma_e, recommended$active_sigma_e * 100),
  "------------------------------------------------------------------\n\n",
  "FINAL PORTFOLIO WEIGHTS\n",
  "------------------------------------------------------------------\n"
)

# Add portfolio weights
final_weights_sorted <- sort(recommended$final_weights, decreasing = TRUE)
for (i in seq_along(final_weights_sorted)) {
  asset <- names(final_weights_sorted)[i]
  weight <- final_weights_sorted[i]
  bar <- paste(rep("#", max(1, round(abs(weight) * 50))), collapse = "")
  if (weight >= 0) {
    output <- paste0(output, sprintf("%-8s %7.2f%% %s\n", asset, weight * 100, bar))
  } else {
    output <- paste0(output, sprintf("%-8s %7.2f%% %s (SHORT)\n", asset, weight * 100, bar))
  }
}
output <- paste0(output, "------------------------------------------------------------------\n", "```\n")

cat(output)

+================================================================+
|         RECOMMENDED TREYNOR-BLACK PORTFOLIO                   |
|         (50% Alpha Shrinkage, No Tracking Risk Limit)         |
+================================================================+

ALLOCATION SUMMARY
------------------------------------------------------------------
Active Portfolio Weight:     75.71%
Market Index Weight:         24.29%
------------------------------------------------------------------

PERFORMANCE METRICS
------------------------------------------------------------------
Market Sharpe Ratio:        0.5964
Portfolio Sharpe Ratio:     1.0444
Sharpe Ratio Improvement:   0.4480
Tracking Risk:               13.21%
------------------------------------------------------------------

ACTIVE PORTFOLIO CHARACTERISTICS
------------------------------------------------------------------
Alpha (Adjusted):           0.1459 (14.59%)
Beta:                       0.7952
Idiosyncratic Vol:          0.1702 (17.02%)
------------------------------------------------------------------

FINAL PORTFOLIO WEIGHTS
------------------------------------------------------------------
MSFT       56.55% ############################
AAPL       37.90% ###################
Market     24.29% ############
PG         14.39% #######
CAT         5.80% ###
NEE         1.16% #
XOM        -0.83% # (SHORT)
JPM        -3.03% ## (SHORT)
AMZN       -3.20% ## (SHORT)
JNJ        -4.25% ## (SHORT)
DIS       -28.78% ############## (SHORT)
------------------------------------------------------------------

Summary of Concepts Covered

output <- paste0(
  "\n```\n",
  "+================================================================+\n",
  "|              KEY TAKEAWAYS FROM CHAPTER 9                      |\n",
  "+================================================================+\n\n",
  "1. SENSITIVITY PROBLEM\n",
  "   - Mean-variance optimization is highly sensitive to expected returns\n",
  "   - Small changes in inputs cause large changes in optimal weights\n",
  "   - This makes pure Markowitz optimization impractical\n\n",
  "2. TREYNOR-BLACK MODEL\n",
  "   - Combines active stock picking with passive index investing\n",
  "   - Weight each security by alpha/sigma^2(epsilon) ratio\n",
  "   - Optimal mix maximizes Sharpe ratio: SR^2_P = SR^2_M + Sum(alpha/sigma_epsilon)^2\n\n",
  "3. TRACKING RISK\n",
  "   - Measures deviation from benchmark (TE = R_P - R_M)\n",
  "   - Managers often constrained by tracking risk limits\n",
  "   - Trade-off: Lower tracking risk means less active management\n\n",
  "4. FORECAST PRECISION\n",
  "   - Analysts tend to be overconfident in alpha forecasts\n",
  "   - Adjust forecasts using historical accuracy: alpha_adj = a + b x alpha_raw\n",
  "   - Shrinkage leads to more realistic portfolio allocations\n\n",
  "5. PRACTICAL IMPLEMENTATION\n",
  "   - Always apply alpha shrinkage (50-60% is common)\n",
  "   - Consider tracking risk constraints if benchmark-relative\n",
  "   - The market index provides diversification and reduces extreme positions\n",
  "```\n"
)

cat(output)

+================================================================+
|              KEY TAKEAWAYS FROM CHAPTER 9                      |
+================================================================+

1. SENSITIVITY PROBLEM
   - Mean-variance optimization is highly sensitive to expected returns
   - Small changes in inputs cause large changes in optimal weights
   - This makes pure Markowitz optimization impractical

2. TREYNOR-BLACK MODEL
   - Combines active stock picking with passive index investing
   - Weight each security by alpha/sigma^2(epsilon) ratio
   - Optimal mix maximizes Sharpe ratio: SR^2_P = SR^2_M + Sum(alpha/sigma_epsilon)^2

3. TRACKING RISK
   - Measures deviation from benchmark (TE = R_P - R_M)
   - Managers often constrained by tracking risk limits
   - Trade-off: Lower tracking risk means less active management

4. FORECAST PRECISION
   - Analysts tend to be overconfident in alpha forecasts
   - Adjust forecasts using historical accuracy: alpha_adj = a + b x alpha_raw
   - Shrinkage leads to more realistic portfolio allocations

5. PRACTICAL IMPLEMENTATION
   - Always apply alpha shrinkage (50-60% is common)
   - Consider tracking risk constraints if benchmark-relative
   - The market index provides diversification and reduces extreme positions

Session Information

sessionInfo()

## R version 4.4.3 (2025-02-28 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] C
## system code page: 65001
## 
## time zone: Asia/Shanghai
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices datasets  utils     methods   base     
## 
## other attached packages:
##  [1] rmdformats_1.0.4           corrplot_0.95             
##  [3] kableExtra_1.4.0           knitr_1.50                
##  [5] quadprog_1.5-8             PerformanceAnalytics_2.0.8
##  [7] quantmod_0.4.28            TTR_0.24.4                
##  [9] xts_0.14.1                 zoo_1.8-14                
## [11] scales_1.4.0               lubridate_1.9.4           
## [13] forcats_1.0.1              stringr_1.6.0             
## [15] dplyr_1.1.4                purrr_1.2.0               
## [17] readr_2.1.6                tidyr_1.3.1               
## [19] tibble_3.3.0               tidyverse_2.0.0           
## [21] ggplot2_4.0.1             
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        generics_0.1.4     renv_1.0.7         xml2_1.5.1        
##  [5] stringi_1.8.7      lattice_0.22-6     hms_1.1.4          digest_0.6.39     
##  [9] magrittr_2.0.4     evaluate_1.0.5     grid_4.4.3         timechange_0.3.0  
## [13] RColorBrewer_1.1-3 bookdown_0.46      fastmap_1.2.0      jsonlite_2.0.0    
## [17] viridisLite_0.4.2  textshaping_1.0.4  jquerylib_0.1.4    cli_3.6.5         
## [21] rlang_1.1.6        withr_3.0.2        cachem_1.1.0       yaml_2.3.11       
## [25] tools_4.4.3        tzdb_0.5.0         curl_7.0.0         vctrs_0.6.5       
## [29] R6_2.6.1           lifecycle_1.0.4    pkgconfig_2.0.3    pillar_1.11.1     
## [33] bslib_0.9.0        gtable_0.3.6       glue_1.8.0         systemfonts_1.3.1 
## [37] xfun_0.54          tidyselect_1.2.1   rstudioapi_0.17.1  farver_2.1.2      
## [41] htmltools_0.5.9    labeling_0.4.3     svglite_2.2.2      rmarkdown_2.30    
## [45] compiler_4.4.3     S7_0.2.1

This workbook was created as a practical companion to Chapter 9: Optimal Portfolio Construction. All examples use real market data and demonstrate the full implementation of the Treynor-Black model with practical adjustments for forecast precision and tracking risk constraints.