#####################################################
# Bass Diffusion Model                              #
# Chapter 13 - Chestnut Ridge Smart Watch Example   #
#####################################################

# This script reproduces the full Chestnut Ridge smart watch Bass diffusion example in R:
#   Step 1. Data setup (cumulative and lagged cumulative sales)
#   Step 2. Estimate the Bass diffusion model (N, p, and q)
#   Step 3. Forecast sales for 50 quarters
#   Step 4. Visualize actual vs. predicted sales
#   Step 5. Use the forecasts for sales planning
# Run it from top to bottom. Each step prints its result to the console or Plots pane.

#####################################################
# Step 1: Data Setup                                #
#####################################################

## Install Packages (if needed - only installs packages that are missing)
packages <- c("ggplot2")
new_packages <- packages[!(packages %in% installed.packages()[, "Package"])]
if (length(new_packages) > 0) install.packages(new_packages)

## Load Packages and Set Seed
library(ggplot2)
set.seed(1)

## Read in the sales data (see Table 13.4)
sales <- read.csv(file.choose()) ## Choose the file sales_data.csv
sales$date <- as.Date(sales$date, format = "%m/%d/%Y")
head(sales)

## Create cumulative sales and lag cumulative sales variables
## N(t-1) is cumulative sales up to the prior quarter, and N(t-1)^2 is its square
## The lag for the first quarter is 0 since there are no sales before launch
sales$cumsales <- cumsum(sales$sales)
sales$cumsales_1 <- c(0, head(sales$cumsales, -1))
sales$cumsales2 <- sales$cumsales^2
sales$cumsales2_1 <- c(0, head(sales$cumsales2, -1))
print(sales)

#####################################################
# Step 2: Estimate the Bass Diffusion Model         #
#####################################################

## Run the Bass Diffusion Model regression (equation 13.5)
## n(t) = a + b*N(t-1) + c*N(t-1)^2
bass <- lm(sales ~ cumsales_1 + cumsales2_1, data = sales)
summary(bass) ## R-squared = 0.7752

## Determine N, p, and q from the regression coefficients
a <- coef(bass)[[1]] ## Intercept                   ## 2.3721
b <- coef(bass)[[2]] ## Coefficient on cumsales_1   ## 0.1179
c <- coef(bass)[[3]] ## Coefficient on cumsales2_1  ## -0.0002692
N1 <- (-b + sqrt(b^2 - 4*a*c)) / (2*c)
N2 <- (-b - sqrt(b^2 - 4*a*c)) / (2*c)
N <- max(N1, N2) ## Market potential (millions of units)
p <- a / N       ## Coefficient of innovation
q <- b + p       ## Coefficient of imitation
print(c(N = N, p = p, q = q)) ## N = 457.34, p = 0.0052, q = 0.1231

#####################################################
# Step 3: Forecast Sales                            #
#####################################################

## Forecast sales for 50 quarters using the Bass diffusion model
## n(t) = p*N + (q - p)*N(t-1) - (q/N)*N(t-1)^2
periods <- 50
forecasts <- data.frame(t = 1:periods,
                        date = seq(min(sales$date), by = "quarter", length.out = periods))
psales <- double(periods); pcumsales <- double(periods + 1)
for (i in 1:periods) {
  psales[i] <- p*N + (q - p)*pcumsales[i] - (q/N)*pcumsales[i]^2
  pcumsales[i+1] <- pcumsales[i] + psales[i]
}
forecasts$psales <- psales
forecasts$pcumsales <- pcumsales[-1] ## Remove the first value which is 0

## Add the actual sales to the forecasts (NA for quarters with no data yet)
forecasts$sales <- sales$sales[match(forecasts$t, sales$period)]
forecasts$cumsales <- sales$cumsales[match(forecasts$t, sales$period)]
print(forecasts)

#####################################################
# Step 4: Visualize the Sales Forecasts             #
#####################################################

## Function to plot actual vs. predicted sales by quarter
## Actual values only exist for the first 24 quarters
plot.forecast <- function(actual, predicted, title, ylab) {
  df <- data.frame(t = rep(forecasts$t, 2),
                   type = rep(c("Actual", "Predicted"), each = periods),
                   value = c(actual, predicted))
  ggplot(df, aes(x = t, y = value, color = type)) +
    geom_line(linewidth = 1, na.rm = TRUE) +
    scale_color_manual(values = c(Actual = "darkorange", Predicted = "steelblue")) +
    scale_x_continuous(breaks = seq(2, periods, by = 2)) +
    labs(title = title, x = "Quarter", y = ylab, color = NULL) +
    theme_minimal() +
    theme(legend.position = "bottom")
}

## Figure 13.4: Actual vs. Predicted Sales by Quarter
plot.forecast(forecasts$sales, forecasts$psales,
              "Actual vs. Predicted Sales Per Quarter", "Sales (millions of units)")

## Figure 13.5: Actual vs. Predicted Cumulative Sales by Quarter
plot.forecast(forecasts$cumsales, forecasts$pcumsales,
              "Actual vs. Predicted Cumulative Sales Per Quarter", "Cumulative Sales (millions of units)")

## Quarter in which predicted sales peak
forecasts[which.max(forecasts$psales), c("t", "date", "psales")] ## Quarter 27 (10/1/2020), 15.28 million

#####################################################
# Step 5: Using Forecasts for Sales Planning        #
#####################################################

## Predicted market share of the Chestnut Ridge smart watch (first choice rule, Chapter 12)
cr.share <- 0.30

## Expected Chestnut Ridge sales in the next quarter (quarter 25)
next.qtr <- forecasts[forecasts$t == max(sales$period) + 1, ]
next.qtr$psales            ## 15.13 million smart watches in the market
next.qtr$psales * cr.share ## 4.54 million Chestnut Ridge smart watches

## Expected Chestnut Ridge sales for each remaining forecast quarter
future <- forecasts[forecasts$t > max(sales$period), c("t", "date", "psales")]
future$cr.sales <- future$psales * cr.share
print(future)
sum(future$cr.sales) ## 73.51 million over quarters 25 to 50