Special Projects in Econ v4

Part 0: Environment Setup

0.1: Libraries

0.2: Process Functions

0.3: Latex Functions

# Export summary stats for latex 
latex.summary.stats = function(stats.list,
                             caption = "Return Statistics Summary",
                             label = "tab:return_stats") {
  
  # Convert to data frame
  stats.df = as.data.frame(do.call(rbind, stats.list)) %>%
    tibble::rownames_to_column("Asset") %>%
    mutate(Asset = toupper(Asset))
  
  # Format numbers
  formatted.df = stats.df %>%
    mutate(across(where(is.numeric), ~ {
      if (cur_column() %in% c("Daily.Mean", "Annual.Mean", "Std.Dev", "Annual.Vol")) {
        formatC(.x, format = "f", digits = 5)
      } else if (cur_column() %in% c("VaR.95%", "CVaR.95%")) {
        paste0(formatC(.x*100, format = "f", digits = 2), "%")
      } else if (cur_column() == "Zero.Rate") {
        scales::percent(.x, accuracy = 0.01)
      } else {
        formatC(.x, format = "f", digits = 2)
      }
    }))
  
  # Generate minimal LaTeX output
  kable(formatted.df, 
        format = "latex",
        booktabs = TRUE,
        caption = caption,
        label = label,
        align = c("l", rep("r", ncol(formatted.df)-1))) %>%
    kable_styling(latex_options = c("striped", "hold_position"))
}
# Save asset return scatter plots as images
# Removes titles 
save.return.scatters = function(returns.list, 
                              folder.name = "2.2.1 Return Scatter Plots",
                              width = 8, 
                              height = 5,
                              dpi = 300) {
  
  # Create directory if it doesn't exist
  if (!dir.exists(folder.name)) {
    dir.create(folder.name, recursive = TRUE)
  }
  
  # Generate plots with titles (for viewing in R)
  plots = plot.individual.returns(returns.list)
  
  # Save each plot without titles
  for (asset in names(plots)) {
    # Create version without title
    p.no.title = plots[[asset]] + 
      labs(title = NULL) +
      theme(plot.title = element_blank())
    
    # Construct file path
    file.name = paste0(asset, " Return Scatter.png")
    file.path = file.path(folder.name, file.name)
    
    # Save plot
    ggsave(file.path, 
           plot = p.no.title,
           width = width,
           height = height,
           dpi = dpi)
  }
  
  # Return the original plots with titles for R viewing
  return(plots)
}

# Save asset return AND price histograms as images
# Removes titles
save.histograms = function(price.hists, return.hists,
                          price.folder = "2.2.2 Price Histograms",
                          return.folder = "2.2.2 Return Histograms",
                          width = 8, height = 5, dpi = 300) {
  
  # Create directories if they don't exist
  dir.create(price.folder, showWarnings = FALSE, recursive = TRUE)
  dir.create(return.folder, showWarnings = FALSE, recursive = TRUE)
  
  # Save price histograms
  for (asset in names(price.hists)) {
    # Remove title from price histogram
    p.price = price.hists[[asset]] + 
      labs(title = NULL) +
      theme(plot.title = element_blank())
    
    ggsave(
      file.path(price.folder, paste0(asset, " Price Histogram.png")),
      plot = p.price,
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  # Save return histograms
  for (asset in names(return.hists)) {
    # Remove title from return histogram
    p.return = return.hists[[asset]] + 
      labs(title = NULL) +
      theme(plot.title = element_blank())
    
    ggsave(
      file.path(return.folder, paste0(asset, " Return Histogram.png")),
      plot = p.return,
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  # Return original plots for viewing in R
  return(list(price.hists = price.hists, return.hists = return.hists))
}

# Save asset price line plots as images
# Removes titles, converts to static plots
save.price.line.plots = function(price.plots,
                               folder.name = "2.2.3 Price Line Plots",
                               width = 10,
                               height = 6,
                               dpi = 300) {
  
  # Create directory if it doesn't exist
  if (!dir.exists(folder.name)) {
    dir.create(folder.name, recursive = TRUE)
  }
  
  # Required packages
  if (!requireNamespace("plotly", quietly = TRUE)) install.packages("plotly")
  if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
  
  # Save each plot
  for (asset in names(price.plots)) {
    # Extract data from Plotly object
    plotly_data = plotly::plotly_data(price.plots[[asset]])
    
    # Create base ggplot
    p = ggplot2::ggplot(plotly_data, ggplot2::aes(
      x = if ("Date" %in% names(plotly_data)) Date else if ("x" %in% names(plotly_data)) x else index
    )) +
      ggplot2::labs(x = "Date", y = "Price") +
      ggplot2::theme_minimal()
    
    # Add price line
    if ("Price" %in% names(plotly_data)) {
      p = p + ggplot2::geom_line(ggplot2::aes(y = Price), color = "#2E8B57", linewidth = 0.8)
    } else if ("y" %in% names(plotly_data)) {
      p = p + ggplot2::geom_line(ggplot2::aes(y = y), color = "#2E8B57", linewidth = 0.8)
    }
    
    # Add moving average if present (looking for MA50 column)
    if ("MA50" %in% names(plotly_data)) {
      p = p + ggplot2::geom_line(ggplot2::aes(y = MA50), 
                                color = "#FF8C00", 
                                linetype = "dashed", 
                                linewidth = 0.8)
    } else if (any(grepl("MA", names(plotly_data)))) {
      ma_col = names(plotly_data)[grepl("MA", names(plotly_data))][1]
      p = p + ggplot2::geom_line(ggplot2::aes(y = .data[[ma_col]]), 
                                color = "#FF8C00", 
                                linetype = "dashed", 
                                linewidth = 0.8)
    }
    
    # Save
    ggplot2::ggsave(
      filename = file.path(folder.name, paste0(asset, " Price Line Plot.png")),
      plot = p,
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  # Return original Plotly objects
  return(price.plots)
}

# Helper function to convert Plotly to ggplot
plotly_to_ggplot = function(plotly_obj) {
  plotly_data = plotly_data(plotly_obj)
  
  ggplot(plotly_data, aes(x = x, y = y)) +
    geom_line(color = "#2E8B57", linewidth = 0.8) +
    labs(x = "Date", y = "Price") +
    theme_minimal() +
    scale_y_continuous(labels = scales::comma)
}

# Save time series for returns as images
# Removes titles, converts to static plots
save.return.ts.plots = function(return.plots,
                              folder.name = "2.2.4 Return Time Series",
                              width = 10,
                              height = 6,
                              dpi = 300) {
  
  # Create directory if it doesn't exist
  if (!dir.exists(folder.name)) {
    dir.create(folder.name, recursive = TRUE)
  }
  
  # Required packages
  if (!requireNamespace("plotly", quietly = TRUE)) install.packages("plotly")
  if (!requireNamespace("ggplot2", quietly = TRUE)) install.packages("ggplot2")
  
  # Save each plot
  for (asset in names(return.plots)) {
    # Extract data from Plotly object
    plotly_data = plotly::plotly_data(return.plots[[asset]])
    
    # Create base ggplot
    p = ggplot2::ggplot(plotly_data, ggplot2::aes(x = Date)) +
      ggplot2::geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
      ggplot2::labs(x = "Date", y = "Log Return") +
      ggplot2::theme_minimal() +
      ggplot2::scale_y_continuous(labels = scales::percent_format(accuracy = 0.1))
    
    # Add return line
    if ("Return" %in% names(plotly_data)) {
      p = p + ggplot2::geom_line(ggplot2::aes(y = Return), color = "darkslategray", linewidth = 0.8)
    }
    
    # Add volatility bands if present
    if (all(c("Volatility", "Return") %in% names(plotly_data))) {
      p = p + 
        ggplot2::geom_ribbon(ggplot2::aes(ymin = -Volatility, ymax = Volatility),
                          fill = "gray70", alpha = 0.2) +
        ggplot2::labs(caption = paste(unique(plotly_data$vol.window), "day rolling volatility"))
    }
    
    # Save
    ggplot2::ggsave(
      filename = file.path(folder.name, paste0(asset, " Return TS with Vol Bands.png")),
      plot = p,
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  # Return original Plotly objects
  return(return.plots)
}

# Save ECDF plots as images
# Removes titles
save.ecdf.plots = function(ecdf.plots,
                         folder.name = "2.2.5 ECDF Plots",
                         width = 8,
                         height = 6,
                         dpi = 300) {
  
  # Create directory if it doesn't exist
  if (!dir.exists(folder.name)) {
    dir.create(folder.name, recursive = TRUE)
  }
  
  # Save each plot
  for (asset in names(ecdf.plots)) {
    # Remove title
    p = ecdf.plots[[asset]] +
      ggtitle(NULL) +
      theme(plot.title = element_blank())
    
    # Save
    ggsave(
      filename = file.path(folder.name, paste0(asset, " ECDF Comparison.png")),
      plot = p,
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  # Return original plots for R viewing
  return(ecdf.plots)
}

# Save density plots as images
# Removes titles
save.density.plots = function(density.plots,
                            folder.name = "2.2.6 Density Plots",
                            width = 8,
                            height = 6,
                            dpi = 300) {
  
  if (!dir.exists(folder.name)) {
    dir.create(folder.name, recursive = TRUE)
  }
  
  for (asset in names(density.plots)) {
    ggsave(
      filename = file.path(folder.name, paste0(asset, " Return Density.png")),
      plot = density.plots[[asset]] + 
        ggtitle(NULL) +
        theme(plot.title = element_blank()),
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  return(density.plots)
}

# Save QQ plots as images
# Removes titles
save.qq.plots = function(qq.plots,
                        folder.name = "2.2.7 Q-Q Plots",
                        width = 8,
                        height = 6,
                        dpi = 300) {
  
  # Create directory if it doesn't exist
  if (!dir.exists(folder.name)) {
    dir.create(folder.name, recursive = TRUE)
  }
  
  # Save each plot
  for (asset in names(qq.plots)) {
    # Remove title but keep all other elements
    p = qq.plots[[asset]] + 
      ggtitle(NULL) +
      theme(plot.title = element_blank())
    
    ggsave(
      filename = file.path(folder.name, paste0(asset, " Q-Q Plot.png")),
      plot = p,
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  # Return original plots for R viewing
  return(qq.plots)
}

# Save volatility scaling comparison plots as images
# Removes titles, re-calculates as static images to ease in saving
save.volatility.plots = function(vol.scaled.returns, 
                               vol.scaled.comparison,
                               folder.name = "2.3 Volatility Scaling Comparison",
                               width = 10,
                               height = 6,
                               dpi = 300) {
  
  # Create output directory
  dir.create(folder.name, showWarnings = FALSE, recursive = TRUE)
  
  # Process each asset
  for (asset in names(vol.scaled.returns)) {
    # Get common dates
    common_dates = as.Date(intersect(index(vol.scaled.returns[[asset]]), 
                          index(vol.scaled.comparison[[asset]])))
    
    # Create aligned data frame
    plot_data = data.frame(
      Date = common_dates,
      GARCH = as.numeric(vol.scaled.returns[[asset]][common_dates]),
      Rolling = as.numeric(vol.scaled.comparison[[asset]][common_dates])
    )
    
    # Create static plot
    p = ggplot2::ggplot(plot_data, ggplot2::aes(x = Date)) +
      ggplot2::geom_line(ggplot2::aes(y = GARCH), 
                        color = "#3A0177",  # Purple
                        linewidth = 0.8) +
      ggplot2::geom_line(ggplot2::aes(y = Rolling), 
                        color = "#FFA500",  # Orange
                        linewidth = 0.8) +
      ggplot2::labs(
        x = "Date",
        y = "Scaled Returns",
        caption = paste("GARCH(1,1) vs 63-Day Rolling |", 
                       length(common_dates), "common observations")
      ) +
      ggplot2::theme_minimal()
    
    # Save the plot
    ggplot2::ggsave(
      filename = file.path(folder.name, paste0(asset, " Volatility Scaling.png")),
      plot = p,
      width = width,
      height = height,
      dpi = dpi
    )
    
    message("Saved ", asset, " (", length(common_dates), " observations)")
  }
  
  message("All plots saved to: ", normalizePath(folder.name))
}

# Save excess return plots as images
# Removes titles, re-calculates as static images to ease in saving
save.excess.return.plots = function(log.returns.list,
                                   folder.name = "3.1 Excess Returns",
                                   width = 10,
                                   height = 6,
                                   dpi = 300) {
  
  # Create output directory
  dir.create(folder.name, showWarnings = FALSE, recursive = TRUE)
  
  # Define asset classes
  equity.assets = c("spy", "djci", "grn")
  fixed.income.assets = c("hyg", "lqd", "krbn")
  
  # Calculate excess returns
  excess.returns.list = list()
  excess.returns.list[equity.assets] = lapply(log.returns.list[equity.assets], 
                                            function(x) x - log.returns.list$shv)
  excess.returns.list[fixed.income.assets] = lapply(log.returns.list[fixed.income.assets], 
                                                  function(x) x - log.returns.list$ief)
  excess.returns.list$ust = log.returns.list$ust - log.returns.list$ief
  
  # Create and save plots
  for (asset in names(excess.returns.list)) {
    # Determine benchmark and color
    if (asset %in% equity.assets) {
      benchmark = "shv"
      line.color = "#4B0082"  # Purple
    } else if (asset %in% fixed.income.assets) {
      benchmark = "ief"
      line.color = "#DAA520"  # Gold
    } else {
      benchmark = "ief"
      line.color = "#4682B4"  # Blue for UST
    }
    
    # Create plot using original function
    p = plot.excess.returns(
      excess.returns = excess.returns.list[[asset]],
      asset.name = asset,
      benchmark.name = benchmark,
      add.mean.line = TRUE,
      add.vol.bands = TRUE,
      window.size = 63,
      add.cumulative = TRUE,
      color.scheme = c("shv" = line.color, "ief" = line.color)
    ) +
      theme(plot.title = element_blank(),
            plot.subtitle = element_blank())  # Remove titles for saved version
    
    # Save plot
    ggsave(
      filename = file.path(folder.name, paste0(toupper(asset), " Excess Returns.png")),
      plot = p,
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  message("Successfully saved all excess return plots to:\n", normalizePath(folder.name))
}

# Save correlation matrix as image
# Removes title, but keeps legend showing correlation gradient 
save.correlation.plots = function(log.returns.list,
                                folder.name = "3.2 Correlation Analysis",
                                width = 10,
                                height = 8,
                                dpi = 300) {
  
  # Create output directory
  dir.create(folder.name, showWarnings = FALSE, recursive = TRUE)
  
  # Define asset classes
  equity.assets = c("spy", "djci", "grn")
  fixed.income.assets = c("hyg", "lqd", "krbn")
  
  # Compute all correlations
  cor.analysis = list(
    full_matrix = compute.cor.matrix(log.returns.list, names(log.returns.list)),
    equity_matrix = compute.cor.matrix(log.returns.list, equity.assets),
    fixed_income_matrix = compute.cor.matrix(log.returns.list, fixed.income.assets),
    rolling_examples = list(
      spy_vs_hyg = compute.rolling.cor(log.returns.list, c("spy", "hyg")),
      lqd_vs_ief = compute.rolling.cor(log.returns.list, c("lqd", "ief"))
    )
  )
  
  # Save heatmaps (unchanged)
  heatmap_files = c(
    "Full Asset Correlation.png",
    "Equity Assets Correlation.png",
    "Fixed Income Correlation.png"
  )
  
  heatmap_data = list(
    cor.analysis$full_matrix,
    cor.analysis$equity_matrix,
    cor.analysis$fixed_income_matrix
  )
  
  for (i in 1:3) {
    p = plot.cor.heatmap(heatmap_data[[i]]) + 
      labs(title = NULL)
    ggsave(
      filename = file.path(folder.name, heatmap_files[i]),
      plot = p,
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  # Save rolling correlation plots (FIXED VERSION)
  rolling_pairs = list(
    list(data = cor.analysis$rolling_examples$spy_vs_hyg, 
         name = "SPY vs HYG Rolling Correlation.png",
         color = "#4B0082"),
    list(data = cor.analysis$rolling_examples$lqd_vs_ief,
         name = "LQD vs IEF Rolling Correlation.png",
         color = "#DAA520")
  )
  
  for (plot in rolling_pairs) {
    # Convert xts to data frame
    plot_data = data.frame(
      Date = index(plot$data),
      Correlation = as.numeric(plot$data)
    )
    
    p = ggplot(plot_data, aes(x = Date, y = Correlation)) +
      geom_line(color = plot$color, linewidth = 0.8) +
      geom_hline(yintercept = 0, linetype = "dashed", color = "gray40") +
      labs(x = NULL, y = "Correlation") +
      ylim(-1, 1) +
      theme_minimal()
    
    ggsave(
      filename = file.path(folder.name, plot$name),
      plot = p,
      width = width,
      height = height/2,
      dpi = dpi
    )
  }
  
  message("Successfully saved all correlation plots to:\n", normalizePath(folder.name))
  return(cor.analysis)
}

# Save beta plots + beta regression with confidence bands as image
# Removes title
save.beta.plots = function(log.returns.list,
                         folder.name = "3.3 Beta Analysis",
                         width = 10,
                         height = 6,
                         dpi = 300) {
  
  # Create output directory
  dir.create(folder.name, showWarnings = FALSE, recursive = TRUE)
  
  # 1. Save static beta table
  beta.table = compute.static.betas(log.returns.list)
  write.csv(beta.table, file.path(folder.name, "Static Beta Table.csv"), row.names = FALSE)
  
  # 2. Save rolling beta charts for all assets
  rolling.beta.dir = file.path(folder.name, "Rolling Betas")
  dir.create(rolling.beta.dir, showWarnings = FALSE)
  
  for (asset in names(log.returns.list)) {
    # Generate plot without title
    p = plot.rolling.beta(log.returns.list, asset) +
      labs(title = NULL)
    
    # Save plot
    ggsave(
      filename = file.path(rolling.beta.dir, paste0(toupper(asset), " Rolling Beta.png")),
      plot = p,
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  # 3. Save regression plots for all assets
  regression.dir = file.path(folder.name, "Regression Plots")
  dir.create(regression.dir, showWarnings = FALSE)
  
  for (asset in setdiff(names(log.returns.list), c("shv", "ief"))) {  # Skip benchmarks
    # Generate plot without title
    p = plot.beta.regression(log.returns.list, asset) +
      labs(title = NULL)
    
    # Save plot
    ggsave(
      filename = file.path(regression.dir, paste0(toupper(asset), " Beta Regression.png")),
      plot = p,
      width = width,
      height = height,
      dpi = dpi
    )
  }
  
  message("Beta analysis saved to:\n", normalizePath(folder.name))
  return(beta.table)  # Return the table for reference
}

Part 1: Asset Identification

# 1. Load assets by their symbol name #########################################
symbol.names = c("ust", "spy", "djci", "grn", "krbn", "hyg", "lqd", "shv", "ief") # Adjust as necessary

# Structural Break dates
break.dates = c("2022-02-24", "2022-12-18", "2023-04-23", "2023-08-29", "2024-04-18", "2020-11-04", "2021-07-14", "2022-08-16")
# Recession data
recessions = getSymbols("USREC", src = "FRED", auto.assign = FALSE)


# 2. Loop to load the data ####################################################
data.list = lapply(symbol.names, function(symbol) {
  raw.data = load.data(symbol)
  renamed.data = rename.columns(raw.data)
  cleaned.data = clean.na(renamed.data)
  return(cleaned.data)
})
names(data.list) = symbol.names
# Data is now stored in data.list


# 3. Align all time series by date ############################################
data.list.aligned = align.data(data.list)
# Remove rows with NA values across any column
data.list.aligned = lapply(data.list.aligned, function(x) na.omit(x))
# Find the common dates
common.dates = Reduce(intersect, lapply(data.list, index))
common.dates = as.Date(common.dates, origin = "1970-01-01")
# range(common.dates) - Check range of common dates
data.list.aligned = lapply(data.list, function(x) x[common.dates])
# Data now is only for days in which all assets have data


# 4. Finding log returns for each asset #######################################
log.returns.list = compute.log.returns(data.list.aligned)

# Scaled returns
vol.scaled.returns = compute.scaled.returns(log.returns.list)

Part 2: Exploratory Analysis

2.1: Simple Summary Statistics

# 2.1.1 Find mean, median, standard deviation, skewness, kurtosis ############
summary.stats.list = simple.summary.statistics(log.returns.list)
summary.stats.list # view in console

## $ust
##      Daily Mean     Annual Mean         Std Dev      Annual Vol        Skewness 
##   -5.062620e-04   -1.275780e-01    9.819116e-03    1.558736e-01    1.852816e-01 
## Excess Kurtosis      VaR 95%.5%        CVaR 95%       Zero Rate    Observations 
##   -1.976198e+00   -1.631311e-02   -2.121645e-02    1.043478e-02    1.150000e+03 
## 
## $spy
##      Daily Mean     Annual Mean         Std Dev      Annual Vol        Skewness 
##    5.206010e-04    1.311914e-01    1.040119e-02    1.651139e-01   -3.172307e-01 
## Excess Kurtosis      VaR 95%.5%        CVaR 95%       Zero Rate    Observations 
##   -1.225584e+00   -1.683330e-02   -2.435606e-02    1.739130e-03    1.150000e+03 
## 
## $djci
##      Daily Mean     Annual Mean         Std Dev      Annual Vol        Skewness 
##    6.996524e-05    1.763124e-02    3.305672e-03    5.247592e-02    3.099838e+00 
## Excess Kurtosis      VaR 95%.5%        CVaR 95%       Zero Rate    Observations 
##    9.214906e+01    0.000000e+00   -2.886261e-04    9.234783e-01    1.150000e+03 
## 
## $grn
##      Daily Mean     Annual Mean         Std Dev      Annual Vol        Skewness 
##    8.469648e-04    2.134351e-01    2.683234e-02    4.259501e-01   -2.804260e-01 
## Excess Kurtosis      VaR 95%.5%        CVaR 95%       Zero Rate    Observations 
##    1.542315e+00   -3.928052e-02   -6.046667e-02    1.304348e-02    1.150000e+03 
## 
## $krbn
##      Daily Mean     Annual Mean         Std Dev      Annual Vol        Skewness 
##    2.732644e-04    6.886263e-02    2.057251e-02    3.265784e-01   -1.539841e+00 
## Excess Kurtosis      VaR 95%.5%        CVaR 95%       Zero Rate    Observations 
##    1.237164e+01   -2.952526e-02   -4.997694e-02    4.347826e-03    1.150000e+03 
## 
## $hyg
##      Daily Mean     Annual Mean         Std Dev      Annual Vol        Skewness 
##   -5.528583e-05   -1.393203e-02    4.874542e-03    7.738095e-02   -1.133693e-01 
## Excess Kurtosis      VaR 95%.5%        CVaR 95%       Zero Rate    Observations 
##    1.995521e+00   -8.005141e-03   -1.147046e-02    7.826087e-03    1.150000e+03 
## 
## $lqd
##      Daily Mean     Annual Mean         Std Dev      Annual Vol        Skewness 
##   -2.024853e-04   -5.102630e-02    5.555306e-03    8.818775e-02    9.329215e-03 
## Excess Kurtosis      VaR 95%.5%        CVaR 95%       Zero Rate    Observations 
##   -1.369065e+00   -9.750599e-03   -1.249420e-02    1.739130e-03    1.150000e+03 
## 
## $shv
##      Daily Mean     Annual Mean         Std Dev      Annual Vol        Skewness 
##   -2.280165e-06   -5.746016e-04    6.165331e-04    9.787160e-03   -5.330251e+00 
## Excess Kurtosis      VaR 95%.5%        CVaR 95%       Zero Rate    Observations 
##    2.758984e+01   -1.818147e-04   -2.030120e-03    2.617391e-01    1.150000e+03 
## 
## $ief
##      Daily Mean     Annual Mean         Std Dev      Annual Vol        Skewness 
##   -2.202760e-04   -5.550955e-02    4.911023e-03    7.796007e-02    1.376115e-01 
## Excess Kurtosis      VaR 95%.5%        CVaR 95%       Zero Rate    Observations 
##   -1.765564e+00   -8.399981e-03   -1.084614e-02    6.086957e-03    1.150000e+03

format.summary.stats(summary.stats.list) # view in viewer

Return Statistics Summary

Asset	Daily Mean	Annual Mean	Std Dev	Annual Vol	Skewness	Excess Kurtosis	VaR 95%.5%	CVaR 95%	Zero Rate	Observations
UST	-0.00	-0.13	0.01	0.16	0.19	-1.98	-0.02	-0.02	0.01	1150.00
SPY	0.00	0.13	0.01	0.17	-0.32	-1.23	-0.02	-0.02	0.00	1150.00
DJCI	0.00	0.02	0.00	0.05	3.10	92.15	0.00	-0.00	0.92	1150.00
GRN	0.00	0.21	0.03	0.43	-0.28	1.54	-0.04	-0.06	0.01	1150.00
KRBN	0.00	0.07	0.02	0.33	-1.54	12.37	-0.03	-0.05	0.00	1150.00
HYG	-0.00	-0.01	0.00	0.08	-0.11	2.00	-0.01	-0.01	0.01	1150.00
LQD	-0.00	-0.05	0.01	0.09	0.01	-1.37	-0.01	-0.01	0.00	1150.00
SHV	-0.00	-0.00	0.00	0.01	-5.33	27.59	-0.00	-0.00	0.26	1150.00
IEF	-0.00	-0.06	0.00	0.08	0.14	-1.77	-0.01	-0.01	0.01	1150.00

# 2.1.2 Export for Latex #####################################################
latex.summary.stats(summary.stats.list)

2.2: Data Visualization

# 2.2.1 Scatter Plot (for prices and returns) ################################
price.plots = plot.individual.prices(data.list.aligned)
# To view price scatterplot:
# price.plots$spy (change as necessary for desired asset)

return.plots = plot.individual.returns(log.returns.list)
# To view returns scatterplot:
return.plots$spy # change as necessary

return.plots$grn

return.plots$djci

# Run to save return scatters
save.return.scatters(log.returns.list)

## $ust

## 
## $spy

## 
## $djci

## 
## $grn

## 
## $krbn

## 
## $hyg

## 
## $lqd

## 
## $shv

## 
## $ief

# 2.2.2 Histograms (for prices and returns) ##################################
# Generate histograms
price.hists = plot.price.histograms(data.list.aligned)
return.hists = plot.return.histograms(log.returns.list)

# To view specific histograms (change as necessary)
price.hists$spy

return.hists$spy

# Run to save price and return histograms
save.histograms(price.hists, return.hists)

## $price.hists
## $price.hists$ust

## 
## $price.hists$spy

## 
## $price.hists$djci

## 
## $price.hists$grn

## 
## $price.hists$krbn

## 
## $price.hists$hyg

## 
## $price.hists$lqd

## 
## $price.hists$shv

## 
## $price.hists$ief

## 
## 
## $return.hists
## $return.hists$ust

## 
## $return.hists$spy

## 
## $return.hists$djci

## 
## $return.hists$grn

## 
## $return.hists$krbn

## 
## $return.hists$hyg

## 
## $return.hists$lqd

## 
## $return.hists$shv

## 
## $return.hists$ief

# 2.2.3 Line Plots for Price (like in v3) ####################################
all.price.plots = plot.all.prices.over.time(
  data.list.aligned,
  show.ma = TRUE,
  ma.period = 50
)
# To view specific plots (change as necessary)
all.price.plots$spy    # View SPY plot

all.price.plots$grn   # View UST plot

# Arrange 2 plots horizontally
# subplot(all.price.plots$spy, all.price.plots$ust, nrows = 1)

# Run to save price and return histograms
save.price.line.plots(all.price.plots)

## $ust
## 
## $spy
## 
## $djci
## 
## $grn
## 
## $krbn
## 
## $hyg
## 
## $lqd
## 
## $shv
## 
## $ief

# Arrange two assets in the same plot
# With structural breaks
carbon.asset.prices = plot.interactive.dual.assets(
  data.list = data.list.aligned,
  break.dates = break.dates,
  colors = c("darkgreen", "purple")
)

# Display the plot
carbon.asset.prices

# Save as html object for presentation
htmlwidgets::saveWidget(
  widget = carbon.asset.prices,  # Your plotly object
  file = "carbon_assets_comparison.html",
  selfcontained = TRUE  # All dependencies included
)

# 2.2.4 Time Series for Returns ##############################################
return.plots = plot.all.returns.over.time(log.returns.list)
return.plots$spy # change as necessary

# With volatility bands
return.plots.vol = plot.all.returns.over.time(
  log.returns.list,
  show.volatility = TRUE,
  vol.window = 21  # 21-day rolling volatility
)
return.plots.vol$spy

# Run to save return time series
save.return.ts.plots(return.plots.vol)

## $ust
## 
## $spy
## 
## $djci
## 
## $grn
## 
## $krbn
## 
## $hyg
## 
## $lqd
## 
## $shv
## 
## $ief

# 2.2.5 ECDF Plots for Returns ################################################
ecdf.plots = plot.ecdf.comparison(log.returns.list)
ecdf.plots$spy  # change as necessary

# Run to save ECDF plots
save.ecdf.plots(ecdf.plots)

## $ust

## 
## $spy

## 
## $djci

## 
## $grn

## 
## $krbn

## 
## $hyg

## 
## $lqd

## 
## $shv

## 
## $ief

# 2.2.6 Kernel Density Plots #################################################
density.plots = plot.kernel.density(log.returns.list)
density.plots$krbn

# Run to save density plots
save.density.plots(density.plots)

## $ust

## 
## $spy

## 
## $djci

## 
## $grn

## 
## $krbn

## 
## $hyg

## 
## $lqd

## 
## $shv

## 
## $ief

# 2.2.7 Q-Q Plots ############################################################
qq.plots = plot.qq.comparison(log.returns.list)
qq.plots$ust # adjust as necessary

# Run to save QQ plots
save.qq.plots(qq.plots)

## $ust

## 
## $spy

## 
## $djci

## 
## $grn

## 
## $krbn

## 
## $hyg

## 
## $lqd

## 
## $shv

## 
## $ief

2.3: Dynamic Volatility Thresholds

# 2.3.1 Scaled Returns Calculation ###########################################

# GARCH-scaled returns for ALL assets
vol.scaled.returns = compute.scaled.returns(
  log.returns.list, 
  method = "GARCH"  # Uses GARCH(1,1) model
)

# Rolling-window scaled returns for COMPARISON 
vol.scaled.comparison = compute.scaled.returns(
  log.returns.list,
  method = "rolling",
  window = 63  # 3-month window
)


# 2.3.2 Comparison Plot ######################################################
# Create comparison data for selected assets
selected.assets = symbol.names  # Modify as needed

# Generate plots for all selected assets
interactive.plots = lapply(selected.assets, create.interactive.plot)
names(interactive.plots) = selected.assets

# View SPY plot
interactive.plots$grn # Change as necessary

# 2.3.3 Key Technical Notes ##################################################

# Color Scheme:
#   Purple: GARCH(1,1)
#   Orange: Rolling 63-Day

# Critical Parameters:
#   GARCH.order = c(1,1)
#   window = 63: Quarterly volatility
#   min.obs = 50


# 2.3.4 Save Plots ###########################################################
# Run to save plots
save.volatility.plots(vol.scaled.returns, vol.scaled.comparison)

Part 3: Financial Analysis

3.1: Excess Returns

# 3.1.1 Excess returns for each asset ######################################################################################
# Importance of Asset-Class Matching 
# Matching benchmarks to the same asset class is critical because:
# 1. Equities and fixed-income assets have fundamentally different risk profiles
# 2. Using mismatched benchmarks (e.g., comparing bonds to equity indices) would:
#    - Overstate "excess returns" for fixed-income during equity bull markets
#    - Understate true credit/spread risk in bonds
# 3. Proper matching isolates the specific risk premia to measure:
#    - Equity risk premium (vs risk-free rate)
#    - Credit spread (vs duration-matched Treasuries)

# 3.1.2 Group Assets by Class ################################################
# NEED TO MODIFY IF ADDING NEW ASSETS
equity.assets = c("spy", "djci", "grn", "krbn")  # Equity-like: stocks, commodities, REITs
fixed.income.assets = c("hyg", "lqd")  # Bonds: high-yield, IG corporates, climate bond
benchmarks = c("shv", "ief") # Risk-free & duration-matched

# Calculate Excess Returns 
excess.returns.list = list()

# 1. Equity-like assets vs SHV (risk-free)
for (asset in equity.assets) {
  excess.returns.list[[asset]] = log.returns.list[[asset]] - log.returns.list$shv
}

# 2. Fixed-income assets vs IEF (duration-matched)
for (asset in fixed.income.assets) {
  excess.returns.list[[asset]] = log.returns.list[[asset]] - log.returns.list$ief
}

# 3. Special case: UST vs IEF (leverage analysis)
excess.returns.list$ust = log.returns.list$ust - log.returns.list$ief

plot.settings = list(
  equity = list(
    benchmark = "shv",
    add.vol.bands = TRUE,
    add.cumulative = TRUE,
    window.size = 63,  # 3-month rolling window
    color = "#4B0082"  # Purple for equity
  ),
  fixed.income = list(
    benchmark = "ief",
    add.vol.bands = TRUE,
    add.cumulative = TRUE,
    window.size = 63,  # 3-month rolling window
    color = "#DAA520"  # Gold for fixed income
  )
)
# Generate all plots
excess.plots = list()

# Equity assets vs SHV
for (asset in equity.assets) {
  excess.plots[[asset]] = plot.excess.returns(
    excess.returns.list[[asset]],
    asset,
    benchmark.name = plot.settings$equity$benchmark,
    add.vol.bands = plot.settings$equity$add.vol.bands,
    add.cumulative = plot.settings$equity$add.cumulative,
    window.size = plot.settings$equity$window.size,
    color.scheme = c("shv" = plot.settings$equity$color)
  )
}

# Fixed income vs IEF
for (asset in fixed.income.assets) {
  excess.plots[[asset]] = plot.excess.returns(
    excess.returns.list[[asset]],
    asset,
    benchmark.name = plot.settings$fixed.income$benchmark,
    add.vol.bands = plot.settings$fixed.income$add.vol.bands,
    add.cumulative = plot.settings$fixed.income$add.cumulative,
    window.size = plot.settings$fixed.income$window.size,
    color.scheme = c("ief" = plot.settings$fixed.income$color)
  )
}

# Special case for UST (leveraged bonds)
excess.plots$ust = plot.excess.returns(
  excess.returns.list$ust,
  "ust",
  benchmark.name = "ief",
  add.vol.bands = TRUE,
  add.cumulative = TRUE,
  window.size = 63,
  color.scheme = c("ief" = "#4682B4")  # Blue for leverage analysis
)

# 3.1.3 View Plots ###########################################################
excess.plots$spy  # SPY vs SHV

excess.plots$hyg  # HYG vs IEF

excess.plots$grn

# 3.1.4 Save Plots #############################################################
# Run to save excess return plots
save.excess.return.plots(log.returns.list)


# 3.1.5 Interactive Plots ######################################################
# Generate all interactive plots with structural breaks
interactive.excess.plots = list()

# Equity assets vs SHV
for (asset in equity.assets) {
  interactive.excess.plots[[asset]] = plot.interactive.excess.returns(
    excess.returns.list[[asset]],
    asset,
    benchmark.name = "shv",
    break.dates = break.dates,
    color = "#4B0082"  # Purple for equity
  )
}

# Fixed income assets vs IEF
for (asset in fixed.income.assets) {
  interactive.excess.plots[[asset]] = plot.interactive.excess.returns(
    excess.returns.list[[asset]],
    asset,
    benchmark.name = "ief",
    break.dates = break.dates,
    color = "#DAA520"  # Gold for fixed income
  )
}

# Special case for UST (leveraged bonds)
interactive.excess.plots$ust = plot.interactive.excess.returns(
  excess.returns.list$ust,
  "ust",
  benchmark.name = "ief",
  break.dates = break.dates,
  color = "#4682B4"  # Blue for leverage analysis
)

# View sample plots
interactive.excess.plots$spy  # View SPY interactive plot

interactive.excess.plots$hyg  # View HYG interactive plot 

interactive.excess.plots$grn  # View GRN interactive plot

# Save all as HTML files
output.dir = "interactive_excess_returns"
if (!dir.exists(output.dir)) dir.create(output.dir)

for (asset in names(interactive.excess.plots)) {
  htmlwidgets::saveWidget(
    widget = interactive.excess.plots[[asset]],
    file = file.path(output.dir, paste0("excess_returns_", asset, ".html")),
    selfcontained = TRUE,
    title = paste(toupper(asset), "Excess Returns")
  )
}

# Quarto embedding example (paste this in your .qmd file)
# ```{=html}
# <iframe src="interactive_excess_returns/excess_returns_spy.html" 
#         width="100%" height="550px" style="border:none">
# </iframe>
# ```

# Save as HTML files
for (asset in names(interactive.excess.plots)) {
  htmlwidgets::saveWidget(
    interactive.excess.plots[[asset]],
    file = paste0("excess_returns_", asset, ".html"),
    selfcontained = TRUE
  )
}

3.2: Correlation Analysis

# 3.2.1 Generate all correlations (all assets, equities, and fixed-income assets) #######################################################################
cor.analysis = list(
  full_matrix = compute.cor.matrix(log.returns.list, names(log.returns.list)),
  equity_matrix = compute.cor.matrix(log.returns.list, equity.assets),
  fixed_income_matrix = compute.cor.matrix(log.returns.list, fixed.income.assets),
  
  rolling_examples = list(
    spy_vs_hyg = compute.rolling.cor(log.returns.list, c("spy", "hyg")),
    lqd_vs_ief = compute.rolling.cor(log.returns.list, c("lqd", "ief"))
  )
)

cor.analysis

## $full_matrix
##              UST        SPY          DJCI          GRN        KRBN        HYG
## UST   1.00000000 0.08391995 -0.0117354890 -0.079478041 -0.03035372 0.42756732
## SPY   0.08391995 1.00000000  0.0961820258  0.155008376  0.23796875 0.73737106
## DJCI -0.01173549 0.09618203  1.0000000000  0.077496451  0.08689711 0.06590468
## GRN  -0.07947804 0.15500838  0.0774964508  1.000000000  0.87711652 0.11457016
## KRBN -0.03035372 0.23796875  0.0868971126  0.877116522  1.00000000 0.21838560
## HYG   0.42756732 0.73737106  0.0659046817  0.114570158  0.21838560 1.00000000
## LQD   0.86337394 0.35418825  0.0001822377  0.006212548  0.07161501 0.67185322
## SHV   0.05509452 0.03904444 -0.0006680272  0.036044080  0.02481405 0.18939433
## IEF   0.97476338 0.09073090 -0.0122067894 -0.071633334 -0.02015878 0.45257520
##               LQD           SHV         IEF
## UST  0.8633739368  0.0550945203  0.97476338
## SPY  0.3541882547  0.0390444430  0.09073090
## DJCI 0.0001822377 -0.0006680272 -0.01220679
## GRN  0.0062125481  0.0360440802 -0.07163333
## KRBN 0.0716150150  0.0248140511 -0.02015878
## HYG  0.6718532163  0.1893943294  0.45257520
## LQD  1.0000000000  0.1485443881  0.88540168
## SHV  0.1485443881  1.0000000000  0.15570246
## IEF  0.8854016838  0.1557024569  1.00000000
## 
## $equity_matrix
##             SPY       DJCI        GRN       KRBN
## SPY  1.00000000 0.09618203 0.15500838 0.23796875
## DJCI 0.09618203 1.00000000 0.07749645 0.08689711
## GRN  0.15500838 0.07749645 1.00000000 0.87711652
## KRBN 0.23796875 0.08689711 0.87711652 1.00000000
## 
## $fixed_income_matrix
##           HYG       LQD
## HYG 1.0000000 0.6718532
## LQD 0.6718532 1.0000000
## 
## $rolling_examples
## $rolling_examples$spy_vs_hyg
##            m.c.seq.row..seq.n...seq.col..drop...FALSE.
## 2021-01-29                                   0.7455191
## 2021-02-01                                   0.7373467
## 2021-02-02                                   0.7473795
## 2021-02-03                                   0.7478035
## 2021-02-04                                   0.7479127
## 2021-02-05                                   0.7474876
## 2021-02-08                                   0.7496898
## 2021-02-09                                   0.7493800
## 2021-02-10                                   0.7548242
## 2021-02-11                                   0.7591040
##        ...                                            
## 2025-02-14                                   0.6498239
## 2025-02-18                                   0.6555032
## 2025-02-19                                   0.6570469
## 2025-02-20                                   0.6542478
## 2025-02-21                                   0.6514648
## 2025-02-24                                   0.6475066
## 2025-02-25                                   0.6394352
## 2025-02-26                                   0.6329395
## 2025-02-27                                   0.6316327
## 2025-02-28                                   0.6343618
## 
## $rolling_examples$lqd_vs_ief
##            m.c.seq.row..seq.n...seq.col..drop...FALSE.
## 2021-01-29                                   0.7322314
## 2021-02-01                                   0.7319750
## 2021-02-02                                   0.7332279
## 2021-02-03                                   0.7311797
## 2021-02-04                                   0.7320267
## 2021-02-05                                   0.7326836
## 2021-02-08                                   0.7296902
## 2021-02-09                                   0.7270867
## 2021-02-10                                   0.7236393
## 2021-02-11                                   0.7222727
##        ...                                            
## 2025-02-14                                   0.9439260
## 2025-02-18                                   0.9463398
## 2025-02-19                                   0.9458633
## 2025-02-20                                   0.9460082
## 2025-02-21                                   0.9450068
## 2025-02-24                                   0.9450417
## 2025-02-25                                   0.9460573
## 2025-02-26                                   0.9453145
## 2025-02-27                                   0.9436629
## 2025-02-28                                   0.9437162

# 3.2.2 Generate Plots #######################################################
cor.plots = list(
  full_heatmap = plot.cor.heatmap(cor.analysis$full_matrix, "Full Asset Correlation"),
  equity_heatmap = plot.cor.heatmap(cor.analysis$equity_matrix, "Equity Assets Correlation"),
  fixed_income_heatmap = plot.cor.heatmap(cor.analysis$fixed_income_matrix, "Fixed Income Correlation"),
  
  rolling_plots = list(
    spy_hyg = ggplot(fortify(cor.analysis$rolling_examples$spy_vs_hyg), 
                     aes(x = Index, y = Value)) +
      geom_line(color = "#4B0082") +
      labs(title = "6-Month Rolling: SPY vs HYG", x = NULL, y = "Correlation") +
      ylim(-1, 1) +
      theme_minimal(),
    
    lqd_ief = ggplot(fortify(cor.analysis$rolling_examples$lqd_vs_ief), 
                    aes(x = Index, y = Value)) +
      geom_line(color = "#DAA520") +
      labs(title = "6-Month Rolling: LQD vs IEF", x = NULL, y = "Correlation") +
      ylim(-1, 1) +
      theme_minimal()
  )
)

# 3.2.3 Access Matrix ########################################################
# cor.analysis$full_matrix # Numerical 
plot.cor.heatmap(cor.analysis$full_matrix) # can modify argument after $ if necessary

# 3.2.4 Save Matrix ##########################################################
# Run to save correlation matrix
save.correlation.plots(log.returns.list)

## $full_matrix
##              UST        SPY          DJCI          GRN        KRBN        HYG
## UST   1.00000000 0.08391995 -0.0117354890 -0.079478041 -0.03035372 0.42756732
## SPY   0.08391995 1.00000000  0.0961820258  0.155008376  0.23796875 0.73737106
## DJCI -0.01173549 0.09618203  1.0000000000  0.077496451  0.08689711 0.06590468
## GRN  -0.07947804 0.15500838  0.0774964508  1.000000000  0.87711652 0.11457016
## KRBN -0.03035372 0.23796875  0.0868971126  0.877116522  1.00000000 0.21838560
## HYG   0.42756732 0.73737106  0.0659046817  0.114570158  0.21838560 1.00000000
## LQD   0.86337394 0.35418825  0.0001822377  0.006212548  0.07161501 0.67185322
## SHV   0.05509452 0.03904444 -0.0006680272  0.036044080  0.02481405 0.18939433
## IEF   0.97476338 0.09073090 -0.0122067894 -0.071633334 -0.02015878 0.45257520
##               LQD           SHV         IEF
## UST  0.8633739368  0.0550945203  0.97476338
## SPY  0.3541882547  0.0390444430  0.09073090
## DJCI 0.0001822377 -0.0006680272 -0.01220679
## GRN  0.0062125481  0.0360440802 -0.07163333
## KRBN 0.0716150150  0.0248140511 -0.02015878
## HYG  0.6718532163  0.1893943294  0.45257520
## LQD  1.0000000000  0.1485443881  0.88540168
## SHV  0.1485443881  1.0000000000  0.15570246
## IEF  0.8854016838  0.1557024569  1.00000000
## 
## $equity_matrix
##             SPY       DJCI        GRN
## SPY  1.00000000 0.09618203 0.15500838
## DJCI 0.09618203 1.00000000 0.07749645
## GRN  0.15500838 0.07749645 1.00000000
## 
## $fixed_income_matrix
##            HYG        LQD       KRBN
## HYG  1.0000000 0.67185322 0.21838560
## LQD  0.6718532 1.00000000 0.07161501
## KRBN 0.2183856 0.07161501 1.00000000
## 
## $rolling_examples
## $rolling_examples$spy_vs_hyg
##            m.c.seq.row..seq.n...seq.col..drop...FALSE.
## 2021-01-29                                   0.7455191
## 2021-02-01                                   0.7373467
## 2021-02-02                                   0.7473795
## 2021-02-03                                   0.7478035
## 2021-02-04                                   0.7479127
## 2021-02-05                                   0.7474876
## 2021-02-08                                   0.7496898
## 2021-02-09                                   0.7493800
## 2021-02-10                                   0.7548242
## 2021-02-11                                   0.7591040
##        ...                                            
## 2025-02-14                                   0.6498239
## 2025-02-18                                   0.6555032
## 2025-02-19                                   0.6570469
## 2025-02-20                                   0.6542478
## 2025-02-21                                   0.6514648
## 2025-02-24                                   0.6475066
## 2025-02-25                                   0.6394352
## 2025-02-26                                   0.6329395
## 2025-02-27                                   0.6316327
## 2025-02-28                                   0.6343618
## 
## $rolling_examples$lqd_vs_ief
##            m.c.seq.row..seq.n...seq.col..drop...FALSE.
## 2021-01-29                                   0.7322314
## 2021-02-01                                   0.7319750
## 2021-02-02                                   0.7332279
## 2021-02-03                                   0.7311797
## 2021-02-04                                   0.7320267
## 2021-02-05                                   0.7326836
## 2021-02-08                                   0.7296902
## 2021-02-09                                   0.7270867
## 2021-02-10                                   0.7236393
## 2021-02-11                                   0.7222727
##        ...                                            
## 2025-02-14                                   0.9439260
## 2025-02-18                                   0.9463398
## 2025-02-19                                   0.9458633
## 2025-02-20                                   0.9460082
## 2025-02-21                                   0.9450068
## 2025-02-24                                   0.9450417
## 2025-02-25                                   0.9460573
## 2025-02-26                                   0.9453145
## 2025-02-27                                   0.9436629
## 2025-02-28                                   0.9437162

# 3.2.5 Interactive Matrix (for presentation) ################################
# Interactive matrix with significance
p.cor = cor.matrix.interactive(
  log.returns.list,
  assets.to.include = symbol.names,
  p.threshold = 0.05,
  r.threshold = 0.3
)

# Access matrix
p.cor

# Save as html for presentation
htmlwidgets::saveWidget(p.cor, "correlation_matrix.html")

3.3: Beta (Market Sensitivity)

# 3.3.1 Static Beta Table #####################################################
beta.table = compute.static.betas(log.returns.list)
beta.table

##                           Asset       Beta         Alpha          R2
## market.excess[valid.obs]1   SPY 1.00000000 -1.170318e-16 1.000000000
## market.excess[valid.obs]4  KRBN 0.46949155  7.574186e-03 0.056314338
## market.excess[valid.obs]3   GRN 0.39600982  1.618291e-01 0.023564986
## market.excess[valid.obs]5   HYG 0.34188968 -5.840688e-02 0.549125268
## market.excess[valid.obs]6   LQD 0.18587658 -7.494392e-02 0.123531960
## market.excess[valid.obs]    UST 0.07742611 -1.372056e-01 0.006739171
## market.excess[valid.obs]7   IEF 0.03972524 -6.016939e-02 0.007239770
## market.excess[valid.obs]2  DJCI 0.03181549  1.401364e-02 0.009671300

# 3.3.2 Rolling Beta Chart ##############################################################################
# Example use:
plot.rolling.beta(log.returns.list, "lqd")

plot.rolling.beta(log.returns.list, "grn")

# 3.3.3 Beta Regression Plot with Confidence Bands ###########################
# Example use:
plot.beta.regression(log.returns.list, "spy")

plot.beta.regression(log.returns.list, "lqd")

# 3.3.4 Save Beta Plots ######################################################
save.beta.plots(log.returns.list)

##                           Asset       Beta         Alpha          R2
## market.excess[valid.obs]1   SPY 1.00000000 -1.170318e-16 1.000000000
## market.excess[valid.obs]4  KRBN 0.46949155  7.574186e-03 0.056314338
## market.excess[valid.obs]3   GRN 0.39600982  1.618291e-01 0.023564986
## market.excess[valid.obs]5   HYG 0.34188968 -5.840688e-02 0.549125268
## market.excess[valid.obs]6   LQD 0.18587658 -7.494392e-02 0.123531960
## market.excess[valid.obs]    UST 0.07742611 -1.372056e-01 0.006739171
## market.excess[valid.obs]7   IEF 0.03972524 -6.016939e-02 0.007239770
## market.excess[valid.obs]2  DJCI 0.03181549  1.401364e-02 0.009671300

# 3.3.5 Interactive Beta (for presentation) ##################################
# Generate plot
p.betas = plot.interactive.betas(
  returns.list = log.returns.list,
  break.dates = break.dates,
  window = 63  # 3-month rolling window
)

# View plot
p.betas

# Save for Quarto
htmlwidgets::saveWidget(
  p.betas, 
  "carbon_betas_comparison.html",
  selfcontained = TRUE
)

3.4: Tail Risk Analysis

# 3.4.1 Determine tail risk ##################################################
# 95%
tail.risk.table = compute.tail.risk(log.returns.list)


# 3.4.2 Visualize tail risk ##################################################
kable(tail.risk.table, align = "c") %>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>%
  row_spec(which(tail.risk.table$`Asset Type` == "Carbon"),
           background = "#E8F5E9", bold = TRUE)

Asset	VaR	ES	Worst Day
ust	-1.63%	-2.12%	-3.64%
spy	-1.68%	-2.44%	-4.45%
djci	0.00%	-0.80%	-3.53%
grn	-3.93%	-6.05%	-16.51%
krbn	-2.95%	-5.00%	-21.97%
hyg	-0.80%	-1.15%	-3.40%
lqd	-0.98%	-1.25%	-2.33%
shv	-0.02%	-0.20%	-0.43%
ief	-0.84%	-1.08%	-1.79%

Part 4: Portfolio Analysis

4.0: Functions & Libraries

# 4.0.0 Load libraries used in portfolio optimization ##########################
library(ggrepel)
library(Matrix)
# Rest of the libraries are in previous section


# 4.0.1 Portfolio Optimization Functions #######################################
# Markowitz optimization
markowitz.optim = function(mu, sigma, target.return = NULL) {
  n = length(mu)
  Dmat = 2 * sigma
  dvec = rep(0, n)
  
  if(is.null(target.return)) {
    Amat = cbind(1, diag(n))
    bvec = c(1, rep(0, n))
    meq = 1
  } else {
    Amat = cbind(1, mu, diag(n))
    bvec = c(1, target.return, rep(0, n))
    meq = 2
  }
  
  Dmat = Dmat + diag(1e-8, n)  # Add small value to diagonal for numerical stability
  
  result = tryCatch({
    sol = solve.QP(Dmat = Dmat, dvec = dvec, Amat = Amat, bvec = bvec, meq = meq)
    list(weights = sol$solution, 
         return = sum(sol$solution * mu),
         risk = sqrt(t(sol$solution) %*% sigma %*% sol$solution),
         status = "success")
  }, error = function(e) {
    list(weights = rep(NA, n),
         return = NA,
         risk = NA,
         status = "failed")
  })
  return(result)
}

# Generate the efficient frontier after optimization - updated for fourth portfolio
generate.frontier = function(returns, portfolio.name) {
  mu = colMeans(returns) * 252  # Annualize returns
  sigma = cov(returns) * 252    # Annualize covariance
  n.assets = ncol(returns)
  risk.free.rate = -0.08
  
  # Ensure positive definite matrix
  if(any(eigen(sigma)$values < 1e-8)) {
    sigma = as.matrix(Matrix::nearPD(sigma)$mat)
  }
  
  target.returns = seq(min(mu)*0.9, max(mu)*1.1, length.out = 50)
  
  frontier = lapply(target.returns, function(ret) {
    pf = markowitz.optim(mu, sigma, ret)
    if(pf$status == "success") {
      data.frame(Return = pf$return, Risk = pf$risk, Target = ret, t(pf$weights))
    } else NULL
  }) %>% bind_rows() %>% arrange(Risk)
  
  frontier$Efficient = cummax(frontier$Return) == frontier$Return
  efficient.frontier = subset(frontier, Efficient)
  
  # Calculate GMV portfolio
  gmv = markowitz.optim(mu, sigma, NULL)
  gmv.df = data.frame(Risk = gmv$risk, Return = gmv$return,
                      Sharpe = (gmv$return - risk.free.rate)/gmv$risk)
  
  # Calculate Tangency portfolio
  sharpe.ratios = (efficient.frontier$Return - risk.free.rate) / efficient.frontier$Risk
  tangency = efficient.frontier[which.max(sharpe.ratios), ]
  tangency$Sharpe = max(sharpe.ratios)
  
  list(
    frontier = efficient.frontier,
    gmv = gmv.df,
    tangency = tangency,
    assets = data.frame(Risk = sqrt(diag(sigma)), Return = mu, Asset = colnames(returns)),
    name = portfolio.name
  )
}


# 4.0.2 Visualization Functions ################################################
# Create static plot of efficient frontier
create.static.plot = function(portfolio.results) {
  colors = c("Base" = "#1f77b4", 
             "Base+GRN" = "#2ca02c",
             "Base+KRBN" = "#d62728",
             "Base+GRN+KRBN" = "#9467bd")
  
  # Ensure all portfolios have plottable data
  plot.data = lapply(names(portfolio.results), function(name) {
    if (nrow(portfolio.results[[name]]$frontier) > 0) {
      data.frame(Risk = portfolio.results[[name]]$frontier$Risk,
                 Return = portfolio.results[[name]]$frontier$Return,
                 Portfolio = name,
                 Type = "Frontier")
    }
  }) %>% bind_rows()
  
  ggplot(plot.data, aes(x = Risk, y = Return, color = Portfolio)) +
    geom_line(linewidth = 1) +
    scale_color_manual(values = colors) +
    labs(title = "Portfolio Efficient Frontiers with Carbon Assets",
         x = "Annualized Risk",
         y = "Annualized Return") +
    theme_minimal()
}

# Create dynamic plot of efficient frontier
create.interactive.plot = function(portfolio.results) {
  colors = c(
    "Base" = "#000000",
    "Base+GRN" = "#ff8da1",
    "Base+KRBN" = "#00ffff",
    "Base+GRN+KRBN" = "#9467bd"  # Orange for combined portfolio
  )
  
  # Initialize plot
  p = plot_ly()
  
  # Add traces for each portfolio
  for (name in names(portfolio.results)) {
    if (nrow(portfolio.results[[name]]$frontier) > 0) {
      p = p %>% add_trace(
        data = portfolio.results[[name]]$frontier,
        x = ~Risk, y = ~Return,
        type = "scatter", mode = "lines",
        line = list(color = colors[name], width = 2),
        name = name,
        hoverinfo = "text",
        text = ~paste("Portfolio:", name,
                      "<br>Return:", round(Return, 4),
                      "<br>Risk:", round(Risk, 4))
      )
      
      # Add GMV point
      p = p %>% add_trace(
        x = portfolio.results[[name]]$gmv$Risk,
        y = portfolio.results[[name]]$gmv$Return,
        type = "scatter", mode = "markers",
        marker = list(
          color = colors[name],
          size = 12,
          symbol = "diamond",
          line = list(width = 1, color = "black")
        ),
        name = paste(name, "GMV"),
        hoverinfo = "text",
        text = paste("Portfolio:", name,
                     "<br>Type: GMV",
                     "<br>Return:", round(portfolio.results[[name]]$gmv$Return, 4),
                     "<br>Risk:", round(portfolio.results[[name]]$gmv$Risk, 4),
                     "<br>Sharpe:", round(portfolio.results[[name]]$gmv$Sharpe, 4))
      )
      
      # Add Tangency point
      p = p %>% add_trace(
        x = portfolio.results[[name]]$tangency$Risk,
        y = portfolio.results[[name]]$tangency$Return,
        type = "scatter", mode = "markers",
        marker = list(
          color = colors[name],
          size = 12,
          symbol = "heart",
          line = list(width = 1, color = "black")
        ),
        name = paste(name, "Tangency"),
        hoverinfo = "text",
        text = paste("Portfolio:", name,
                     "<br>Type: Tangency",
                     "<br>Return:", round(portfolio.results[[name]]$tangency$Return, 4),
                     "<br>Risk:", round(portfolio.results[[name]]$tangency$Risk, 4),
                     "<br>Sharpe:", round(portfolio.results[[name]]$tangency$Sharpe, 4))
      )
    }
  }
  
  # Add assets (only from Base portfolio)
  base.assets = portfolio.results[["Base"]]$assets
  p = p %>% add_trace(
    x = base.assets$Risk,
    y = base.assets$Return,
    type = "scatter", mode = "markers",
    marker = list(color = "black", size = 8),
    text = base.assets$Asset,
    hoverinfo = "text",
    name = "Assets"
  )
  
  # Layout with adjusted title position
  p %>% layout(
    title = list(
      text = "<b>Portfolio Comparison with Carbon Assets</b>",
      x = 0.05,
      y = 0.97,  # Lowered from default (0.98)
      xanchor = "left",
      yanchor = "top",
      font = list(size = 20)
    ),
    margin = list(t = 80),  # Increased top margin
    xaxis = list(title = "Risk (Standard Deviation)"),
    yaxis = list(title = "Expected Return"),
    hoverlabel = list(bgcolor = "white", font = list(size = 12)),
    legend = list(orientation = "h", y = -0.2)
  )
}

4.1: Data Preparation

# 4.1.1 Define the components ##################################################
# Tickers used in construction of base portfolio
base.tickers = c("AAPL", "MSFT", "AMZN", "GOOG", "META", "TSLA", "JPM", 
                 "V", "WMT", "UNH", "SPY") 
carbon.tickers = c("GRN", "KRBN")
all.tickers = c(base.tickers, carbon.tickers)


# 4.1.2 Download data ##########################################################
# Handles if data cannot be found for any symbol
successful.tickers = c()
for(ticker in all.tickers) {
  tryCatch({
    getSymbols(ticker, from = Sys.Date() - 365*3, auto.assign = TRUE)
    successful.tickers = c(successful.tickers, ticker)
  }, error = function(e) {
    message(paste("Could not download data for", ticker))
  })
}


# 4.1.3 Create portfolio (only includes successful tickers) ####################
available.carbon = intersect(carbon.tickers, successful.tickers)
all.portfolios = list(
  "Base" = base.tickers,
  "Base+GRN" = c(base.tickers, "GRN")["GRN" %in% successful.tickers],
  "Base+KRBN" = c(base.tickers, "KRBN")["KRBN" %in% successful.tickers],
  "Base+GRN+KRBN" = c(base.tickers, available.carbon)
)

# Remove any portfolios that ended up with no carbon assets
all.portfolios = all.portfolios[
  sapply(all.portfolios, function(x) {
    is.base.portfolio = identical(names(all.portfolios)[which(sapply(all.portfolios, identical, x))], "Base")
    has.carbon.assets = length(x) > length(base.tickers)
    is.base.portfolio | has.carbon.assets
  })
]

# 4.1.4 Prepare returns function ###############################################
prepare.returns = function(tickers) {
  prices = do.call(merge, lapply(tickers, function(x) Ad(get(x))))
  colnames(prices) = tickers
  na.omit(Return.calculate(prices, method = "arithmetic"))
}

# 4.1.5 Calculate returns ######################################################
returns.list = lapply(all.portfolios, prepare.returns)

4.2 Efficient Frontier Visualization

# 4.2.1 Generate frontiers #####################################################
# Generate frontiers
portfolio.results = lapply(names(all.portfolios), function(name) {
  generate.frontier(returns.list[[name]], name)
})
names(portfolio.results) = names(all.portfolios)


# 4.2.2 Display plots ##########################################################
portfolio.colors = c(
  "Base" = "#1f77b4",
  "Base+GRN" = "#2ca02c", 
  "Base+KRBN" = "#d62728",
  "Base+GRN+KRBN" = "#9467bd"
)

static.plot = create.static.plot(portfolio.results)
interactive.plot = create.interactive.plot(portfolio.results)

interactive.plot

# 4.2.3 Save interactive plot ##################################################
htmlwidgets::saveWidget(
    widget = interactive.plot,
    file = "efficient_frontier.html",
    selfcontained = TRUE
  )

Part 5: Presentation Outputs

Part 6: Essay Appendix Figures

6.0: Essay Export Functions

# Price plot (without title)
plot.dual.assets.paper = function(data.list, asset1 = "GRN", asset2 = "KRBN", 
                                  ma.period = 50, break.dates = NULL,
                                  colors = c("#1f77b4", "#ff7f0e")) {
  
  # Extract price data
  asset1.data = Ad(data.list[[tolower(asset1)]])
  asset2.data = Ad(data.list[[tolower(asset2)]])
  
  # Calculate moving averages
  asset1.ma = rollmean(asset1.data, ma.period, fill = NA, align = "right")
  asset2.ma = rollmean(asset2.data, ma.period, fill = NA, align = "right")
  
  # Create data frame with all series
  plot.data = data.frame(
    Date = index(asset1.data),
    Asset1.Price = as.numeric(asset1.data),
    Asset1.MA = as.numeric(asset1.ma),
    Asset2.Price = as.numeric(asset2.data),
    Asset2.MA = as.numeric(asset2.ma)
  )
  
  # Calculate scaling factor for secondary axis
  scale.factor = mean(asset1.data, na.rm = TRUE) / mean(asset2.data, na.rm = TRUE)
  
  # Create plotly plot with legend
  p = plot_ly(plot.data, x = ~Date) %>%
    add_lines(y = ~Asset1.Price, name = paste(asset1, "Price"),
              line = list(color = colors[1], width = 2),
              hoverinfo = "text",
              text = ~paste(asset1, "<br>Date:", Date,
                            "<br>Price:", round(Asset1.Price, 2),
                            "<br>MA:", round(Asset1.MA, 2))) %>%
    add_lines(y = ~Asset1.MA, name = paste(asset1, "MA"),
              line = list(color = colors[1], width = 2, dash = "dash"),
              hoverinfo = "text",
              text = ~paste(asset1, "MA<br>Date:", Date,
                            "<br>MA:", round(Asset1.MA, 2))) %>%
    add_lines(y = ~Asset2.Price * scale.factor, name = paste(asset2, "Price"),
              line = list(color = colors[2], width = 2),
              yaxis = "y2",
              hoverinfo = "text",
              text = ~paste(asset2, "<br>Date:", Date,
                            "<br>Price:", round(Asset2.Price, 2),
                            "<br>MA:", round(Asset2.MA, 2))) %>%
    add_lines(y = ~Asset2.MA * scale.factor, name = paste(asset2, "MA"),
              line = list(color = colors[2], width = 2, dash = "dash"),
              yaxis = "y2",
              hoverinfo = "text",
              text = ~paste(asset2, "MA<br>Date:", Date,
                            "<br>MA:", round(Asset2.MA, 2)))
  
  # Add layout with axis titles and compact legend
  p = layout(p,
             margin = list(l = 60, r = 80, b = 60, t = 40, pad = 4),  # Adjusted margins
             xaxis = list(title = "Date"),
             yaxis = list(
               title = paste(asset1, "Price ($)"),
               color = colors[1],
               automargin = TRUE
             ),
             yaxis2 = list(
               title = paste(asset2, "Price ($)"),
               overlaying = "y",
               side = "right",
               color = colors[2],
               automargin = TRUE
             ),
             legend = list(
               x = 0.5,
               y = -0.25,  # Position below plot
               orientation = "h",
               xanchor = "center",
               font = list(size = 10)  # Smaller font for compactness
             ),
             hovermode = "x unified"
  )
  
  # Add structural breaks if provided
  if(!is.null(break.dates)) {
    breaks = data.frame(Date = as.Date(break.dates),
                        Y = max(plot.data$Asset1.Price, na.rm = TRUE))
    p = p %>%
      add_segments(data = breaks,
                   x = ~Date, xend = ~Date,
                   y = 0, yend = ~Y,
                   line = list(color = "red", dash = "dot", width = 1),
                   name = "Break Point",
                   hoverinfo = "text",
                   text = ~paste("Structural Break<br>Date:", Date),
                   showlegend = FALSE)
  }
  
  return(p)
}

# Excess return plot (without title)
plot.paper.excess.returns = function(excess.returns, 
                                     asset.name, 
                                     benchmark.name,
                                     break.dates = NULL,
                                     window.size = 63,
                                     color = "#4B0082") {
  
  # Prepare data
  plot.data = data.frame(
    Date = index(excess.returns),
    Excess.Return = as.numeric(coredata(excess.returns))
  ) %>% 
    filter(!is.na(Excess.Return)) %>%
    mutate(
      Rolling.SD = zoo::rollapply(Excess.Return, window.size, sd, fill = NA, align = "right"),
      Cumulative = cumsum(Excess.Return)
    )
  
  # Create plot without title
  p = plot_ly(plot.data, x = ~Date) %>%
    
    # Main excess return line
    add_lines(y = ~Excess.Return, name = "Daily Excess Return",
              line = list(color = color, width = 1.5),
              hoverinfo = "text",
              text = ~paste0("<b>", toupper(asset.name), " vs ", toupper(benchmark.name), "</b>",
                             "<br>Date: ", Date,
                             "<br>Excess Return: ", round(Excess.Return*100, 2), "%",
                             "<br>Cumulative: ", round(Cumulative*100, 2), "%")) %>%
    
    # Zero line
    add_lines(y = 0, name = "Zero", 
              line = list(color = "gray", dash = "dot", width = 1),
              hoverinfo = "none") %>%
    
    # Volatility bands
    add_ribbons(ymin = ~-Rolling.SD, ymax = ~Rolling.SD,
                name = paste0(window.size, "-day Volatility"),
                fillcolor = "rgba(200,200,200,0.2)",
                line = list(color = "rgba(0,0,0,0)"),
                hoverinfo = "text",
                text = ~paste0("Volatility Band: ±", round(Rolling.SD*100, 2), "%")) %>%
    
    # Cumulative return (secondary axis)
    add_lines(y = ~Cumulative, name = "Cumulative Excess",
              line = list(color = "darkgreen", dash = "dot", width = 1),
              yaxis = "y2",
              hoverinfo = "text",
              text = ~paste0("Cumulative: ", round(Cumulative*100, 2), "%"))
  
  # Add structural breaks if provided
  if(!is.null(break.dates)) {
    break.df = data.frame(Date = as.Date(break.dates))
    p = p %>%
      add_segments(
        data = break.df,
        x = ~Date, xend = ~Date,
        y = min(plot.data$Excess.Return, na.rm = TRUE), 
        yend = max(plot.data$Excess.Return, na.rm = TRUE),
        name = "Break Point",
        line = list(color = "red", dash = "dot", width = 1.2),
        hoverinfo = "text",
        text = ~paste("Structural Break<br>Date:", Date),
        showlegend = FALSE
      )
  }
  
  # Layout without title
  p = p %>%
    layout(
      yaxis = list(title = "Daily Excess Return",
                   tickformat = ".2%",
                   hoverformat = ".2%"),
      yaxis2 = list(title = "Cumulative Return",
                    overlaying = "y",
                    side = "right",
                    tickformat = ".2%",
                    hoverformat = ".2%",
                    showgrid = FALSE),
      legend = list(orientation = "h", y = -0.2),
      hovermode = "x unified",
      margin = list(t = 40, r = 80)  # Reduced top margin (no title), kept right margin
    )
  
  return(p)
}

# 3 Month Rolling Betas (without title)
plot.interactive.betas.paper = function(returns.list, 
                                  asset1 = "GRN", 
                                  asset2 = "KRBN",
                                  market.asset = "SPY",
                                  risk.free.asset = "SHV",
                                  break.dates = NULL,
                                  window = 63) {
  
  # Calculate excess returns
  market.excess = returns.list[[tolower(market.asset)]] - returns.list[[tolower(risk.free.asset)]]
  asset1.excess = returns.list[[tolower(asset1)]] - returns.list[[tolower(risk.free.asset)]]
  asset2.excess = returns.list[[tolower(asset2)]] - returns.list[[tolower(risk.free.asset)]]
  
  # Compute rolling betas
  calc.rolling.beta = function(asset.excess) {
    merged = na.omit(merge.xts(asset.excess, market.excess))
    rollapply(merged, width = window,
              FUN = function(x) coef(lm(x[,1] ~ x[,2]))[2],
              by.column = FALSE)
  }
  
  beta1 = calc.rolling.beta(asset1.excess)
  beta2 = calc.rolling.beta(asset2.excess)
  
  # Prepare plot data
  plot.data = data.frame(
    Date = index(beta1),
    GRN.beta = as.numeric(beta1),
    KRBN.beta = as.numeric(beta2)
  )
  
  # Create plot
  p = plot_ly(plot.data, x = ~Date) %>%
    add_lines(y = ~GRN.beta, name = paste(asset1, "Beta"),
              line = list(color = "#4B0082", width = 2),
              hoverinfo = "text",
              text = ~paste(asset1, "Beta<br>Date:", Date, 
                            "<br>Beta:", round(GRN.beta, 3))) %>%
    add_lines(y = ~KRBN.beta, name = paste(asset2, "Beta"),
              line = list(color = "#DAA520", width = 2),
              hoverinfo = "text",
              text = ~paste(asset2, "Beta<br>Date:", Date,
                            "<br>Beta:", round(KRBN.beta, 3))) %>%
    add_lines(y = ~1, name = "Market Beta",
              line = list(color = "gray", dash = "dot"))
  
  # Add structural breaks
  if(!is.null(break.dates)) {
    breaks.df = data.frame(Date = as.Date(break.dates))
    p = p %>%
      add_segments(
        data = breaks.df,
        x = ~Date, xend = ~Date,
        y = min(plot.data$GRN.beta, plot.data$KRBN.beta, na.rm = TRUE),
        yend = max(plot.data$GRN.beta, plot.data$KRBN.beta, na.rm = TRUE),
        name = "Break Point",
        line = list(color = "red", dash = "dot", width = 1),
        hoverinfo = "text",
        text = ~paste("Structural Break<br>Date:", Date),
        showlegend = FALSE
      )
  }
  
  # Final layout without title
  p %>% layout(
    yaxis = list(title = "Beta Coefficient"),
    xaxis = list(title = ""),
    hovermode = "x unified",
    legend = list(orientation = "h", y = -0.2),
    margin = list(t = 30, r = 50, b = 60)  # Reduced top margin (t) since no title
  )
}

6.1: Essay Outputs

# 6.1.1 Price Dynamics ####
carbon.asset.prices.paper = plot.dual.assets.paper(
  data.list = data.list.aligned,
  break.dates = break.dates,
  colors = c("darkgreen", "purple")
)

carbon.asset.prices.paper

# 6.1.2 Excess Returns ####
paper.excess.plots = list()
for (asset in equity.assets) {
  paper.excess.plots[[asset]] = plot.paper.excess.returns(
    excess.returns.list[[asset]],
    asset,
    benchmark.name = "shv",
    break.dates = break.dates,
    color = "#191970"  # Purple for equity
  )
}

# Save as images
paper.excess.plots$grn

paper.excess.plots$krbn

# 6.1.3 Correlation Plot ####
plot.cor.heatmap(cor.analysis$full_matrix)

# 6.1.4 Market Betas ####
# Beta table
beta.table

##                           Asset       Beta         Alpha          R2
## market.excess[valid.obs]1   SPY 1.00000000 -1.170318e-16 1.000000000
## market.excess[valid.obs]4  KRBN 0.46949155  7.574186e-03 0.056314338
## market.excess[valid.obs]3   GRN 0.39600982  1.618291e-01 0.023564986
## market.excess[valid.obs]5   HYG 0.34188968 -5.840688e-02 0.549125268
## market.excess[valid.obs]6   LQD 0.18587658 -7.494392e-02 0.123531960
## market.excess[valid.obs]    UST 0.07742611 -1.372056e-01 0.006739171
## market.excess[valid.obs]7   IEF 0.03972524 -6.016939e-02 0.007239770
## market.excess[valid.obs]2  DJCI 0.03181549  1.401364e-02 0.009671300

# Rolling betas
p.betas.paper = plot.interactive.betas.paper(
  returns.list = log.returns.list,
  break.dates = break.dates,
  window = 63  # 3-month rolling window
)
# View plot
p.betas.paper

# 6.1.5 Efficient Frontier ####
static.plot = create.static.plot(portfolio.results)
interactive.plot = create.interactive.plot(portfolio.results)