# Load packages
library(tidyverse)
library(tidyquant)
library(moments)
library(ggrepel)
library(scales)
library(stringr)
library(timetk)
library(PerformanceAnalytics)

1. Import Prices

symbols <- c("AAPL", "MSFT", "GOOG", "AMZN", "TSLA", "NFLX")

prices <- tq_get(
  x = symbols,
  get = "stock.prices",
  from = "2012-01-01",
  to   = "2022-12-31"
)

2. Convert to Monthly log Returns

asset_returns_tbl <- prices %>%
  group_by(symbol) %>%
  tq_transmute(select = adjusted,
               mutate_fun = periodReturn,
               period = "monthly",
               type = "log") %>%
  ungroup()
# Create portfolio returns
weights <- rep(1/6, 6)

portfolio_returns_tbl <- asset_returns_tbl %>%
  tq_portfolio(
    assets_col = symbol,
    returns_col = monthly.returns,
    weights = weights,
    col_rename = "returns"
  )
# Step 1: Calculate kurtosis for the portfolio
portfolio_kurt_tbl <- portfolio_returns_tbl %>%
  tq_performance(Ra = returns, performance_fun = table.Stats) %>%
  select(Kurtosis)
# Step 2: Histogram of portfolio returns
portfolio_returns_tbl %>%
  ggplot(aes(x = returns)) +
  geom_histogram(binwidth = 0.01, fill = "cornflowerblue", alpha = 0.5) +
  geom_vline(xintercept = 0.0003, color = "green", size = 1) +
  annotate("text", x = 0.0015, y = 13, label = "Risk-free rate", angle = 90)

Step 3: Scatterplot of Expected Return vs. Kurtosis

mean_kurt_tbl <- asset_returns_tbl %>%
  group_by(symbol) %>%
  summarise(mean = mean(monthly.returns),
            kurt = kurtosis(monthly.returns)) %>%
  ungroup()

portfolio_stats <- portfolio_returns_tbl %>%
  summarise(mean = mean(returns), kurt = kurtosis(returns)) %>%
  mutate(symbol = "portfolio")

mean_kurt_tbl <- bind_rows(mean_kurt_tbl, portfolio_stats)
# Plot
mean_kurt_tbl %>%
  ggplot(aes(x = kurt, y = mean, color = symbol)) +
  geom_point() +
  geom_text_repel(aes(label = symbol)) +
  theme(legend.position = "none") +
  labs(x = "Kurtosis", y = "Expected Return") +
  scale_y_continuous(labels = scales::percent_format(accuracy = 0.1))

Step 4: Rolling 24-Month Kurtosis Line Plot

# Define the rolling window
window <- 24

rolling_kurt_tbl <- portfolio_returns_tbl %>%
  tq_mutate(
    select     = returns,
    mutate_fun = rollapply,
    width      = window,
    FUN        = kurtosis,
    col_rename = "kurt"
  ) %>%
  na.omit() %>%
  select(date, kurt)
# Plot Rolling Kurtosis
rolling_kurt_tbl %>%
  ggplot(aes(x = date, y = kurt)) +
  geom_line(color = "cornflowerblue") +
  scale_y_continuous(breaks = seq(-1, 4, 0.5)) +
  scale_x_date(breaks = scales::pretty_breaks(n = 7)) +
  labs(
    x = NULL,
    y = "Kurtosis",
    title = paste0("Rolling ", window, " Month Kurtosis")
  ) +
  theme(plot.title = element_text(hjust = 0.5)) +
  annotate(
    "text",
    x = as.Date("2017-07-01"),
    y = 3.3,
    label = str_glue("Downside risk\nskyrocketed toward end of 2017"),
    size = 5,
    color = "red"
  )

Step 5: Calculate CAPM Beta

# 5.1 – Get Market Returns (SPY)
market_returns_tbl <- tq_get("SPY",
                             from = "2012-01-01",
                             to   = "2022-12-31") %>%
  tq_transmute(
    select     = adjusted,
    mutate_fun = periodReturn,
    period     = "monthly",
    type       = "log"
  ) %>%
  slice(-1) %>%
  rename(returns = monthly.returns)
# 5.2 – Join Market and Portfolio Returns
portfolio_market_returns_tbl <- left_join(
  market_returns_tbl,
  portfolio_returns_tbl,
  by = "date"
) %>%
  rename(
    market_returns    = returns.x,
    portfolio_returns = returns.y
  )
# Convert to xts format
portfolio_xts <- portfolio_market_returns_tbl %>%
  select(date, portfolio_returns) %>%
  tk_xts(date_var = date)

market_xts <- portfolio_market_returns_tbl %>%
  select(date, market_returns) %>%
  tk_xts(date_var = date)

# Calculate CAPM beta using PerformanceAnalytics
capm_beta <- CAPM.beta(Ra = portfolio_xts, Rb = market_xts)

capm_beta
## [1] 1.274199

#Step 6

# Scatterplot of portfolio vs market returns
portfolio_market_returns_tbl %>%
  ggplot(aes(x = market_returns, y = portfolio_returns)) +
  geom_point(color = "cornflowerblue") +
  geom_smooth(method = "lm", se = FALSE, size = 1.5,
              color = palette_light()[[3]]) +
  labs(
    x = "Market Returns",
    y = "Portfolio Returns",
    title = "CAPM Scatterplot: Portfolio vs. Market Returns"
  )

#Step 7

# Linear regression for fitted values
library(broom)

actual_vs_fitted_tbl <- portfolio_market_returns_tbl %>%
  lm(portfolio_returns ~ market_returns, data = .) %>%
  augment() %>%
  mutate(date = portfolio_market_returns_tbl$date)
# Pivot to long format
actual_fitted_long_tbl <- actual_vs_fitted_tbl %>%
  select(date, portfolio_returns, .fitted) %>%
  pivot_longer(cols = c(portfolio_returns, .fitted),
               names_to = "type",
               values_to = "returns")
# Plot actual vs fitted returns
actual_fitted_long_tbl %>%
  ggplot(aes(x = date, y = returns, color = type)) +
  geom_line() +
  labs(
    x = "Date",
    y = "Returns",
    title = "Actual vs Fitted Portfolio Returns (CAPM)"
  )