##############################################################
# Chapter 15 - Using Marketing Mix Models to Optimize the    #
#              Marketing Mix                                 #
# Chestnut Ridge HDTV Marketing Mix Model (R)                #
##############################################################

# This script reproduces the full Chestnut Ridge HDTV marketing mix model example in R:
#   Step 1. Install packages (if needed), load packages, and set seed
#   Step 2. Read in the marketing mix data
#   Step 3. Look at means of variables
#   Step 4. Create natural log, lag, and weekday variables
#   Step 5. Check for unit root (augmented Dickey-Fuller test)
#   Step 6. Check for multicollinearity
#   Step 7. Run the marketing mix regression
#   Step 8. Create predicted values for quantity and revenue
#   Step 9. Convert date to a Date so it plots on a time axis
#   Step 10. Plot marketing spend and pricing over time (Figure 15.6)
#   Step 11. Plot predicted vs. actual quantity sold (Figure 15.7)
#   Step 12. Plot predicted vs. actual revenue (Figure 15.8)
# Run it from top to bottom. Each step prints its result to the console or Plots pane.


##############################################################
# Step 1: Install packages (if needed), load packages, and   #
#         set seed                                           #
##############################################################

packages <- c("tseries", "ggplot2", "tidyr")
new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new_packages) > 0) install.packages(new_packages)

library(tseries)
library(ggplot2)
library(tidyr)
set.seed(1)


##############################################################
# Step 2: Read in the marketing mix data                     #
##############################################################

mmix <- read.csv(file.choose())  ## Choose the file mmix_data.csv


##############################################################
# Step 3: Look at means of variables                         #
##############################################################

## (date is text, so its mean is NA and R prints a warning)
lapply(sapply(mmix, mean), round, 2)


##############################################################
# Step 4: Create natural log, lag, and weekday variables     #
##############################################################

mmix$ln_quantity <- log(mmix$quantity)
mmix$ln_price <- log(mmix$price)
mmix$ln_digital_ad <- log(mmix$digital_ad)
mmix$ln_digital_search <- log(mmix$digital_search)
mmix$ln_print <- log(mmix$print + 1)  ## Some time periods have 0 spend
mmix$ln_tv <- log(mmix$tv)

## Lag of ln_quantity (day 1 has no prior day, so it is set to 0)
mmix$lln_quantity <- c(0, mmix$ln_quantity[1:length(mmix$ln_quantity) - 1])

## Day-of-week variable (used as dummy variables in the regression)
## Note: weekdays() returns day names in your computer's language setting.
## On an English-language system the days sort alphabetically, so Friday
## becomes the baseline day in the regression. On a non-English system the
## day names (and therefore the baseline day and the dummy variable labels)
## will differ. To get English day names, first run:
## Sys.setlocale("LC_TIME", "English")  ## Windows
## Sys.setlocale("LC_TIME", "en_US.UTF-8")  ## Mac/Linux
mmix$weekdays <- weekdays(as.Date(mmix$date))


##############################################################
# Step 5: Check for unit root (augmented Dickey-Fuller test) #
##############################################################

adf.test(mmix$ln_quantity)


##############################################################
# Step 6: Check for multicollinearity                        #
##############################################################

cor_vars <- c("ln_quantity", "lln_quantity", "ln_price", "ln_digital_ad",
              "ln_digital_search", "ln_print", "ln_tv")
cor_data <- mmix[cor_vars]
cor_table <- cor(cor_data)
round(cor_table, 2)

## Combine ln_digital_ad and ln_digital_search
mmix$ln_digital <- log(mmix$digital_ad + mmix$digital_search)


##############################################################
# Step 7: Run the marketing mix regression                   #
##############################################################

mmix_reg <- lm(ln_quantity ~ lln_quantity + ln_price + ln_digital + ln_print +
               ln_tv + factor(weekdays), data = mmix)
summary(mmix_reg)


##############################################################
# Step 8: Create predicted values for quantity and revenue   #
##############################################################

mmix$pred_quantity <- exp(predict(mmix_reg))
mmix$pred_revenue <- mmix$price * mmix$pred_quantity


##############################################################
# Step 9: Convert date to a Date so it plots on a time axis  #
##############################################################

mmix$date <- as.Date(mmix$date)


##############################################################
# Step 10: Plot marketing spend and pricing over time        #
#          (Figure 15.6)                                     #
##############################################################

mix_long <- pivot_longer(mmix,
                         cols = c(digital_ad, digital_search, print, tv, price),
                         names_to = "variable", values_to = "value")
mix_long$variable <- factor(mix_long$variable,
                            levels = c("digital_ad", "digital_search", "print",
                                       "tv", "price"),
                            labels = c("Digital Ad ($)", "Digital Search ($)",
                                       "Print ($)", "TV ($)", "Price ($)"))

ggplot(mix_long, aes(x = date, y = value)) +
  geom_line(color = "#2a78d6", linewidth = 0.4) +
  facet_grid(variable ~ ., scales = "free_y", switch = "y") +
  scale_x_date(date_breaks = "2 months", date_labels = "%b %Y") +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Marketing Spend and Pricing Over Time",
       x = "Date", y = NULL) +
  theme_minimal() +
  theme(strip.placement = "outside",
        strip.text.y.left = element_text(angle = 0, hjust = 1),
        axis.text.x = element_text(angle = 45, hjust = 1),
        panel.grid.minor = element_blank())


##############################################################
# Step 11: Plot predicted vs. actual quantity sold (Figure   #
#          15.7)                                             #
##############################################################

qty_long <- pivot_longer(mmix, cols = c(quantity, pred_quantity),
                         names_to = "series", values_to = "value")
qty_long$series <- factor(qty_long$series,
                          levels = c("quantity", "pred_quantity"),
                          labels = c("Actual Quantity", "Predicted Quantity"))

ggplot(qty_long, aes(x = date, y = value, color = series, linetype = series)) +
  geom_line(linewidth = 0.5) +
  scale_color_manual(values = c("#2a78d6", "#eb6834")) +
  scale_linetype_manual(values = c("solid", "dashed")) +
  scale_x_date(date_breaks = "2 months", date_labels = "%b %Y") +
  scale_y_continuous(labels = scales::comma) +
  labs(title = "Predicted vs. Actual Quantity Sold",
       x = "Date", y = "Quantity Sold", color = NULL, linetype = NULL) +
  theme_minimal() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45, hjust = 1),
        panel.grid.minor = element_blank())


##############################################################
# Step 12: Plot predicted vs. actual revenue (Figure 15.8)   #
##############################################################

rev_long <- pivot_longer(mmix, cols = c(revenue, pred_revenue),
                         names_to = "series", values_to = "value")
rev_long$series <- factor(rev_long$series,
                          levels = c("revenue", "pred_revenue"),
                          labels = c("Actual Revenue", "Predicted Revenue"))

ggplot(rev_long, aes(x = date, y = value, color = series, linetype = series)) +
  geom_line(linewidth = 0.5) +
  scale_color_manual(values = c("#2a78d6", "#eb6834")) +
  scale_linetype_manual(values = c("solid", "dashed")) +
  scale_x_date(date_breaks = "2 months", date_labels = "%b %Y") +
  scale_y_continuous(labels = scales::dollar) +
  labs(title = "Predicted vs. Actual Revenue",
       x = "Date", y = "Revenue", color = NULL, linetype = NULL) +
  theme_minimal() +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45, hjust = 1),
        panel.grid.minor = element_blank())