Simple Stock Over/Undervaluation Analyzer

################################################################################
### Fetch data on the selected stock and analyze for under/overvaluation.  This
### is a simple method which involves determining the long-term trend of the
### stock, and seeing if it is currently above/below that trend.   It should
### work decently on stocks that are "well behaved" (long track record of
### continuous gains, and a high r^2)
################################################################################

### Stock ticker for the stock of interest
stock <- "DG"

### What are the start/end dates you are interested in
date_start <- as.Date("1980-01-01")
date_end   <- as.Date(Sys.Date())

################################################################################
### No changes should be necessary below here.
################################################################################

### Load necessary R packages
library(tidyverse)
library(lubridate)
library(quantmod)
library(broom)
library(scales)
library(cowplot)

### Function to add recession bars to graphs
geom_recession_bars <- function(date_start, date_end, fill = "darkseagreen4") {

  date_start <- as.Date(date_start, origin = "1970-01-01")
  date_end   <- as.Date(date_end,   origin = "1970-01-01")

  recessions_tibble <-
    tribble(
      ~peak,     ~trough,
      "1953-07-01", "1954-05-01",
      "1957-08-01", "1958-04-01",
      "1960-04-01", "1961-02-01",
      "1969-12-01", "1970-11-01",
      "1973-11-01", "1975-03-01",
      "1980-01-01", "1980-07-01",
      "1981-07-01", "1982-11-01",
      "1990-07-01", "1991-03-01",
      "2001-03-01", "2001-11-01",
      "2007-12-01", "2009-06-01",
      "2020-02-01", "2020-04-01"
    ) %>%
    mutate(peak   = as.Date(peak), trough = as.Date(trough))

  recessions_trim <- recessions_tibble %>%
    filter(peak   >= min(date_start) & trough <= max(date_end))

  if (nrow(recessions_trim) > 0) {
    recession_bars <- geom_rect(data = recessions_trim, inherit.aes = F,
                                fill = fill, alpha = 0.25,
                                aes(xmin = as.Date(peak,   origin = "1970-01-01"),
                                    xmax = as.Date(trough, origin = "1970-01-01"),
                                    ymin = -Inf, ymax = +Inf))
  } else {
    recession_bars <- geom_blank()
  }
}

### Fetch stock prices from Yahoo Finance
prices <-
  getSymbols(stock, env = NULL, auto.assign = FALSE, src = "yahoo",
             from = date_start, to = date_end) %>%
  fortify.zoo() %>%
  as_tibble() %>%
  select(Index, 2) %>%
  rename(date = Index,
         value = 2) %>%
  arrange(date) %>%
  mutate(log_value = log(value)) %>%
  mutate(decimal_date = decimal_date(date))

### Do a linear regression and pluck out the coefficients
fit_prices <- lm(log_value ~ decimal_date, data = prices %>% filter(value != 0))

coeff1 <- fit_prices %>% tidy() %>% filter(term == "(Intercept)") %>% pull("estimate")
coeff2 <- fit_prices %>% tidy() %>% filter(term == "decimal_date") %>% pull("estimate")
adj_r2 <- fit_prices %>% glance() %>% pull("adj.r.squared")

growth <- paste((coeff2 * 100) %>% round(digits = 1), "%/yr", sep = "")
doubling_time <- paste((log(2) / log(1 + coeff2)) %>% round(digits = 1), "yrs")

### Use the regression coefficients to compute the trend
prices <-
  prices %>%
  mutate(fit_log_value = coeff1 + coeff2 * decimal_date) %>%
  mutate(fit_value = exp(fit_log_value)) %>%
  mutate(value_per = 1 - (fit_value / value))

valuation <- prices %>% arrange(date) %>% tail(n = 1) %>% pull(value_per)

# Determine number of years of trend growth represented by over/undervaluation
years <- (abs(valuation) / coeff2) %>% round(digits = 1)

# Format string to display valuation in useful form
valuation <- (100 * abs(valuation)) %>% round(digits = 1)

if (valuation > 0) {
   valuation_txt <- paste("overvalued by ",  valuation, "% (", years, " years growth)", sep = "")
} else {
   valuation_txt <- paste("undervalued by ", valuation, "% (", years, " years growth)", sep = "")
}


##########################################################################
### Graph 1:  Price of stock
##########################################################################

p_norm <-
  ggplot(data = prices) +
  theme_bw() +
  theme(legend.title = element_blank()) +
  geom_line(aes(x = date, y = value)) +
  geom_line(aes(x = date, y = fit_value)) +
  geom_recession_bars(min(prices$date), max(prices$date)) +
  scale_x_date(breaks = pretty_breaks(10)) +
  scale_y_continuous(breaks = pretty_breaks(10)) +
  labs(title = paste(stock, " Stock Price, trend growth = ", growth, ", doubling time = ", doubling_time, sep = ""),
       caption = "", x = "", y = "US $")

##########################################################################
### Graph 2:  log(Price of stock)
##########################################################################

caption  <- paste("Data retrieved from Yahoo Finance on", format(Sys.time(), "%B %d, %Y at %I:%M %p %Z"))

p_log <-
  ggplot(data = prices) +
  theme_bw() +
  theme(legend.title = element_blank()) +
  geom_line(aes(x = date, y = log_value)) +
  geom_line(aes(x = date, y = fit_log_value)) +
  geom_recession_bars(min(prices$date), max(prices$date)) +
  scale_x_date(breaks = pretty_breaks(10)) +
  scale_y_continuous(breaks = pretty_breaks(10)) +
  labs(title = paste("log(", stock, " Stock Price), adj R^2 = ", adj_r2 %>% round(digits = 3), sep = ""),
       caption = "", x = "", y = "US $")

##########################################################################
### Graph 3:  Standard deviation above/below trend
##########################################################################

# Compute mean / standard deviation of % above/below trend
mean_stock <- prices %>% pull(value_per) %>% mean()
sd_stock   <- prices %>% pull(value_per) %>% sd()

p_trend <-
  ggplot(data = prices) +
  theme_bw() +
  theme(legend.title = element_blank()) +
  geom_line(aes(x = date, y = 1 - (fit_value / value))) +
  geom_hline(yintercept = 0, linetype = "dotted") +
  geom_recession_bars(min(prices$date), max(prices$date)) +
  scale_x_date(breaks = pretty_breaks(10)) +
  scale_y_continuous(breaks = pretty_breaks(10),
                     labels = scales::percent_format(accuracy = 1),
                     sec.axis = sec_axis(~ (. - mean_stock) / sd_stock,
                                           name = "Z-Score",
                                           breaks = pretty_breaks(6))
                     ) +
  labs(title = paste(stock, " Stock Price % above/below trend, ", valuation_txt, sep = ""),
       caption = caption, x = "", y = "%")

### Print final graph
print(plot_grid(p_norm, p_log, p_trend, nrow = 3, ncol = 1, align = "hv"))