---
title: "ACCT 420: Advanced linear regression"
author: 
  - name: Dr. Richard M. Crowley
    email: rcrowley@smu.edu.sg
    affiliation: https://rmc.link/
title-slide-attributes:
  data-background-video: "../Backgrounds/title.mp4"
  data-background-video-loop: "true"
  data-background-size: stretch
format:
  revealjs:
    theme: default
    css: quarto_white.css
    navigation-mode: vertical
    controls: true
    transition: slide
    background-transition: fade
    transition-speed: default
    width: 1536
    height: 864
    logo: ../Logos/logo2.png
    code-link: true
    code-line-numbers: false
    chalkboard:
      theme: whiteboard
    auto-stretch: false
revealjs-plugins:
  - attribution
execute:
  echo: true
  warning: false
  message: false
editor: source
---

# Front Matter {data-background-video="../Backgrounds/section.mp4" data-background-video-loop="true"}

<!--# Note that the downlit and xml2 pacakges are required to compile -->
<!--# Installed extensions include `attribution`, `qrcode`, and `material-icons` -->

```{r helpers, echo=FALSE, warning=FALSE, message=F}
library(knitr)
library(kableExtra)

#' Pretty formatted HTML tables
#' 
#' @param df The dataframe to format as a table
#' @param names Names to use for columns if not using the df column names. Vector
#' @param highlight_cols Columns to highlight. Default is c(1) for highlighting the first row. Use c() to remove.
#' @param highlight_rows Rows to highlight. Default is c() to not highlight any rows.
#' @param color_rows Rows to highlight in color. Default is c() to not highlight any rows.
#' @param highlight_end_rows Rows to highlight. Uses a reverse indexing where 1 is the last row, 2 is second to last, etc.
#' @param full Use to span the full width of a slide. FALSE by default.
#' @param fixed_header Whether the first row should be treated as a header row. TRUE by default.
html_df <- function(df, names=NULL, highlight_cols=c(1), highlight_rows=c(), color_rows=c(), highlight_end_rows=c(), full=F, fixed_header=T) {
  if(!length(names)) {
    names=colnames(df)
  }
  kbl(df,"html", col.names = names, align = c("l",rep('c',length(cols)-1))) %>%
    kable_paper(c("striped","hover"), full_width=full, fixed_thead=F, html_font = "\"Source Sans Pro\", Helvetica, sans-serif") %>%
    {if(length(highlight_cols)) column_spec(., highlight_cols,bold=T, background = "#ffffff99") else .} %>%
    {if(length(highlight_rows)) row_spec(., highlight_rows,bold=T, background = "#ffffff99") else .} %>%
    {if(length(color_rows)) row_spec(., color_rows,bold=T, background = "#aaaaff99") else .} %>%
    {if(length(highlight_end_rows)) row_spec(., nrow(df) + 1 - highlight_end_rows,bold=T, background = "#ffffff99") else .}
}

library(tidyverse)
```

<!--# List packages here that are not being correctly highlighted -->
<!--# The usual cause is calling the library from a echo=FALSE cell -->

```{r}
#| include: false

library(tidyverse)
df <- read.csv("../../Data/Session_3-1.csv", stringsAsFactors=FALSE)

# load in relevant data from Session 2
load("../../Data/Session_2_export.RData")
```

<!--# code to create the exercise data for shiny -->

```{r}
#| include: false
#| eval: false

wmt <- filter(df, tic == "WMT") %>%
  select(tic, fyearq, fqtr, revtq) %>%
  mutate(revtq_lag=lag(revtq))
write_csv(wmt, '../shiny/S3/data/S3.csv')



wmt_training <- filter(wmt, fyearq < 2015)
wmt_testing <- filter(wmt, fyearq >= 2015)
mod1 <- lm(revtq ~ revtq_lag, data=wmt_training)
mod2 <- lm(revtq ~ revtq_lag:factor(fqtr), data=wmt_training)

saveRDS(mod1, "../shiny/S3/data/S3_mod1.rds")
saveRDS(mod2, "../shiny/S3/data/S3_mod2.rds")

# For session 4:
# Data is from: https://fred.stlouisfed.org/series/JHDUSRGDPBR
df_fed <- read_csv('../../Data/Session_4_Fed.csv') %>%
  rename(recession = 2)
  
library(lubridate)
df_fed <- df_fed %>%
  mutate(date = ymd(DATE),
         fyearq = year(date),
         fqtr = quarter(date)) %>%
  select(fyearq, fqtr, recession)
wmt <- wmt %>% left_join(df_fed)
write_csv(wmt, '../shiny/S4/data/S4.csv')

# For session 5
wmt <- wmt %>% select(-revtq_lag) %>%
  mutate(revtq_yoy = (revtq - lag(revtq, 4)) / lag(revtq, 4),
         recession_lag = lag(recession),
         revtq_down = factor(ifelse(revtq_yoy<0,1,0), levels=c(0,1)),
         revtq_yoy_lag = lag(revtq_yoy),
         date = fyearq + (fqtr-1)/4)

s5mod1 <- glm(revtq_down ~ revtq_yoy_lag, data=df, family=binomial)

wmt2 <- wmt
wmt2$pred <- predict(s5mod1, wmt2, type="response")
wmt2 <- wmt2 %>% filter(!is.na(revtq_down), !is.na(pred))

write_csv(wmt, '../shiny/S5/data/S5_1.csv')
write_csv(wmt2, '../shiny/S5/data/S5_2.csv')
saveRDS(s5mod1, "../shiny/S5/data/S5_mod1.rds")

```

## Learning objectives {data-background="../Backgrounds/default.jpg"}

:::: columns

::: {.column width="50%"}

![](../Course_Diagrams/Map-3.png){fig-alt="Roadmap for the course"}

:::

::: {.column width="50%"}
- <b>Theory:</b>
    - Further understand stats treatments
        - Panel data
        - Fixed effects
        - Time (seasonality)
- <b>Application:</b>
    - Predicting revenue quarterly and weekly
- <b>Methodology:</b>
    - Univariate
    - Linear regression (OLS)
    - Visualization
:::

::::

# Application: Quarterly retail revenue {data-background-video="../Backgrounds/section.mp4"}

## The question {data-background="../Backgrounds/default.jpg"}

> How can we predict quarterly revenue for retail companies, leveraging our knowledge of such companies?

 More specifically...

- Consider time dimensions
    - What matters:
        - Last quarter?
        - Last year?
        - Other time frames?
    - Cyclicality

# Time and OLS {data-background-video="../Backgrounds/section.mp4"}

## Time matters a lot for retail {data-background="../Backgrounds/default.jpg"}

![](../Figures/SG_GSS.jpg){width="500px" fig-align="center"} ![](../Figures/SG_1111.jpg){width="500px" fig-align="center"}
![](../Figures/SG_black_friday.jpg){width="500px" fig-align="center"} ![](../Figures/SG_xmas_2016.jpg){width="500px" fig-align="center"}

## How to capture time effects? {data-background="../Backgrounds/default.jpg"}

:::: {.columns}

::: {.column width="50%"}
- Autoregression: our main focus today
    - Regress $y_t$ on earlier value(s) of itself
        - Last quarter, last year, etc.
- Controlling for time directly in the model
    - We can do this using something called a *fixed effect*
        - This absorbs the *noise* in a model due to a particular grouping in the data
:::

::: {.column width="50%"}
![](../Figures/calendar.png){width="400px" fig-align="center"}
:::

::::

## Noise in data {data-background="../Backgrounds/default.jpg"}

> Statistical noise is random error in the data

- Many sources of noise:
    - Factors we neglected to included in the model
    - Error in measurement
        - Accounting measurement!
    - Unexpected events / shocks
    - Groups within the data that behave differently

> Noise is OK, but the more we remove, the better!

## Removing noise: Revisiting SG real estate revenue {data-background="../Backgrounds/case.jpg"}

- Different companies may behave slightly differently
    - Control for this using a *Fixed Effect* (ISIN uniquely identifies companies)

```{r}
forecast3.1 <-
  lm(revt_lead ~ revt + act + che + lct + dp + ebit + factor(isin),
     data=df_clean[df_clean$fic=="SGP",])
# n=15 to prevent outputting every fixed effect
library(broom)
print(tidy(forecast3.1), n=15)
```

## Removing noise: Singapore model {data-background="../Backgrounds/case.jpg"}

```{r}
glance(forecast3.1)
anova(forecast3, forecast3.1, test="Chisq")
```

> This only a bit better.  Why?  There was another source of noise within Singapore real estate companies

## Another way to do fixed effects {data-background="../Backgrounds/aside.jpg"}

- The library `{fixest}` has `feols()`: **f**ixed **e**ffects **OLS**
    - Better for complex models
    - Extremely efficient computationally

```{r, message=F, error=F, warning=F}
library(fixest)
forecast3.2 <-
  feols(revt_lead ~ revt + act + che + lct + dp + ebit | isin,
       data=df_clean[df_clean$fic=="SGP",])
summary(forecast3.2)
```

## Why exactly would we use fixed effects? {data-background="../Backgrounds/reference.jpg"}

:::: {.columns}

::: {.column width="50%"}
- Fixed effects are used when the average of $\hat{y}$ varies by some group in our data
    - In our problem, the average revenue of each firm is different
- Fixed effects absorb this difference

```{r, echo=F, message=F, fig.height=4}
library(tidyverse)
df_clean %>%
  filter(fic=="SGP") %>%
  group_by(isin) %>%
  mutate(mean_revt_lead=mean(revt_lead, na.rm=T)) %>%
  slice(1) %>%
  ungroup() %>%
  ggplot(aes(x=mean_revt_lead)) + geom_histogram(aes(y = ..density..)) +  geom_density(alpha=.4, fill="#FF6666")
```
:::

::: {.column width="50%"}
- Further reading:
    - Introductory Econometrics by Jeffrey M. Wooldridge

![](../Figures/Econometrics_text.jpg){height="400px" fig-align="center"}
:::

::::

# Data preparation {data-background-video="../Backgrounds/section.mp4"}

## The data {data-background="../Backgrounds/default.jpg"}

- From quarterly reports
- Two sets of firms:
    - US "Hypermarkets & Super Centers" (GICS: 30101040)
    - US "Multiline Retail" (GICS: 255030)
- Data from Compustat - Capital IQ > North America - Daily > Fundamentals Quarterly

![](../Figures/walmart.jpg){height="300px" fig-align="center"}

## Formalization {data-background="../Backgrounds/default.jpg"}

1. Question
    - How can we predict quarterly revenue for large retail companies?
2. Hypothesis (just the alternative ones)
    1. Current quarter revenue helps predict next quarter revenue
    2. 3 quarters ago revenue helps predict next quarter revenue (Year-over-year)
    3. Different quarters exhibit different patterns (seasonality)
    4. A long-run autoregressive model helps predict next quarter revenue
3. Prediction
    - Use OLS for all the above, $t$-tests for coefficients
    - Hold out sample: 2018-2023

## Variable generation {data-background="../Backgrounds/coding.jpg"}

```{r, message=F, warning=F}
library(tidyverse)  # As always
library(plotly)  # interactive graphs
library(lubridate)  # import some sensible date functions

# Generate quarter over quarter growth "revtq_gr"
df <- df %>% group_by(gvkey) %>% mutate(revtq_gr=revtq / lag(revtq) - 1) %>% ungroup()

# Generate year-over-year growth "revtq_yoy"
df <- df %>% group_by(gvkey) %>% mutate(revtq_yoy=revtq / lag(revtq, 4) - 1) %>% ungroup()

# Generate first difference "revtq_d"
df <- df %>% group_by(gvkey) %>% mutate(revtq_d=revtq - lag(revtq)) %>% ungroup()

# Generate a proper date
# Date was YYYY-MM-DD, can be converted from text to date easily
df$date <- ymd(df$datadate)  # From lubridate
df$qtr <- quarter(df$date)   # From lubridate
```

- Use mutate for variables using lags
- `ymd()` from `{lubridate}` is handy to convert any date listing year, then month, then day.
    - It also has `ydm()`, `mdy()`, `myd()`, `dmy()` and `dym()`
    - It can handle quarters, times, and date-times as well
    - [Cheat sheet](https://rawgit.com/rstudio/cheatsheets/master/lubridate.pdf)


## Base R Date manipulation {data-background="../Backgrounds/coding.jpg"}

- `as.Date()` can take a date formatted as "YYYY/MM/DD" and convert to a proper date value
    - You can convert other date types using the `format=` argument
        - i.e., "DD.MM.YYYY" is format code "%d.%m.%Y"
        - [Full list of date codes](https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/strptime)
        
## Example output {data-background="../Backgrounds/coding.jpg"}

```{r, echo=F}
html_df(head(df[,c("conm","date","revtq","revtq_gr", "revtq_yoy", "revtq_d")]))

head(df[,c("conm","date", "datadate")])
```

## Creating 8 quarters (2 years) of lags {data-background="../Backgrounds/coding.jpg"}

- Eight is easy to do by copying and pasting code, but...
    - What if we wanted 100 lags?
        - That would be a lot of copying
            - The code for this is in the appendix
- Three solutions, though none is straightforward
    1. Use sapply and a function to reassemble output into a data frame
        - We will ignore this -- option 2 is always better
    2. Use `{purrr}` to loop into a data frame
    3. Use advanced dplyr programming with quosures
- To compare the methods, we'll make a copy of the data frame

```{r}
# Make a copy of our data frame to compare later
df2 <- df
```

## Method 2: Map with Purrr {data-background="../Backgrounds/code.jpg"}

```{r, message=F}
# Approach #2: Mixing purrr with dplyr across()
library(purrr)

multi_lag_map <- function(df, lags, var, postfix="") {
  new_columns <- map(lags,
    function(x) df %>% 
      group_by(gvkey) %>%
      transmute(across(all_of(var),
                       .fns = list(~ lag(., x)),
                       .names = paste0('{col}', postfix, x) ) ) %>%
      ungroup() %>%
      select(-gvkey)
      )
  cbind(df, list_cbind(new_columns))
}
```

- `map()` is the key function. It takes a vector and iterates over it applying a specified function, returning a list of data frames
- `transmute()`: like mutate except it extracts the column it makes automatically
- `cbind_list()`: takes a list of data frames and combines them into 1 data frame
- `cbind()` is a base R function to combine 

## Method 3: Advanced programming {data-background="../Backgrounds/code.jpg"}

```{r, message=F}
# Approach #3: Advanced programming using quosures...
library(rlang)
multi_lag <- function(df, lags, var, postfix="") {
  var <- enquo(var)
  quosures <- map(lags, ~quo(lag(!!var, !!.x))) %>%
    set_names(paste0(quo_text(var), postfix, lags))
  return(ungroup(mutate(group_by(df, gvkey), !!!quosures)))
}
```

- `setNames()`: allows for storing a value and name simultaneously

::: {.callout-warning appearance="simple"}

Quosures are beyond the scope of the course. However, they are very useful if:

1) You need to do non-standard calculations at scale
2) You are programming a library.

:::

## Comparing the methods {data-background="../Backgrounds/default.jpg"}

2. Using `{purrr}` + `{dplyr}`

```{r}
df <- multi_lag_map(df, 1:8, 'revtq', "_l")  # Generate lags "revtq_l#"
df <- multi_lag_map(df, 1:8, 'revtq_gr')     # Generate changes "revtq_gr#"
df <- multi_lag_map(df, 1:8, 'revtq_yoy')    # Generate year-over-year changes "revtq_yoy#"
df <- multi_lag_map(df, 1:8, 'revtq_d')      # Generate first differences "revtq_d#"
```

3. Using advanced programming techniques

```{r}
df2 <- multi_lag(df2, 1:8, revtq, "_l")  # Generate lags "revtq_l#"
df2 <- multi_lag(df2, 1:8, revtq_gr)     # Generate changes "revtq_gr#"
df2 <- multi_lag(df2, 1:8, revtq_yoy)    # Generate year-over-year changes "revtq_yoy#"
df2 <- multi_lag(df2, 1:8, revtq_d)      # Generate first differences "revtq_d#"
```

- Verify that the output is the same

```{r}
all(df==df2, na.rm=T)
```

## Example output {data-background="../Backgrounds/coding.jpg"}

```{r, echo=F}
html_df(head(df[,c("conm","date","revtq","revtq_l1", "revtq_l2", "revtq_l3", "revtq_l4")]))
```

## Clean and split into training and testing {data-background="../Backgrounds/coding.jpg"}

```{r}
# Clean the data: Replace NaN, Inf, and -Inf with NA
df <- df %>%
  mutate(across(where(is.numeric), ~replace(., !is.finite(.), NA)))

# Split into training and testing data
# Training data: We'll use data released before 2018
train <- filter(df, year(date) < 2018)

# Testing data: We'll use data released 2018 through 2023
test <- filter(df, year(date) >= 2018)
```

- Same cleaning function as last week:
    - Replaces all `NaN`, `Inf`, and `-Inf` with `NA`
- `year()` comes from `{lubridate}`

# Univariate stats {data-background-video="../Backgrounds/section.mp4"}

## Univariate stats {data-background="../Backgrounds/default.jpg"}

- To get a better grasp on the problem, looking at univariate stats can help
    - Summary stats (using `summary()`)
    - Correlations using `cor()`
    - Plots using your preferred package such as `{ggplot2}`

```{r, message=F, warning=F}
summary(df[,c("revtq","revtq_gr","revtq_yoy", "revtq_d","qtr")])
```

## ggplot2 for visualization {data-background="../Backgrounds/reference.jpg"}

```{r, include=F}
# These functions are a bit ugly, but can construct many charts quickly
# eval(parse(text=var)) is just a way to convert the string name to a variable reference
# Density plot for 1st to 99th percentile of data
plt_dist <- function(df,var) {
  df %>%
    filter(eval(parse(text=var)) < quantile(eval(parse(text=var)),0.99, na.rm=TRUE),
           eval(parse(text=var)) > quantile(eval(parse(text=var)),0.01, na.rm=TRUE)) %>%
    ggplot(aes(x=eval(parse(text=var)))) + 
    geom_density() + xlab(var)
}
```

```{r, include=F}
# Density plot for 1st to 99th percentile of both columns
plt_bar <- function(df,var) {
  df %>%
    filter(eval(parse(text=var)) < quantile(eval(parse(text=var)),0.99, na.rm=TRUE),
           eval(parse(text=var)) > quantile(eval(parse(text=var)),0.01, na.rm=TRUE)) %>%
    ggplot(aes(y=eval(parse(text=var)), x=qtr)) + 
    geom_bar(stat = "summary", fun.y = "mean") + ylab(var) + xlab("Calendar quarter")
}
```

```{r, include=F}
# Scatter plot with lag for 1st to 99th percentile of data
plt_sct <- function(df,var1, var2) {
  df %>%
    filter(eval(parse(text=var1)) < quantile(eval(parse(text=var1)),0.99, na.rm=TRUE),
           eval(parse(text=var2)) < quantile(eval(parse(text=var2)),0.99, na.rm=TRUE),
           eval(parse(text=var1)) > quantile(eval(parse(text=var1)),0.01, na.rm=TRUE),
           eval(parse(text=var2)) > quantile(eval(parse(text=var2)),0.01, na.rm=TRUE)) %>%
    ggplot(aes(y=eval(parse(text=var1)), x=eval(parse(text=var2)), color=factor(qtr))) + 
    geom_point() + geom_smooth(method = "lm") + ylab(var1) + xlab(var2)
}
```

- The next slides will use some custom functions using `{ggplot2}`
- `{ggplot2}` has an odd syntax:
    - It doesn't use pipes (`%>%`), but instead adds everything together (`+`)

```{r, eval=F}
library(ggplot2)  # or tidyverse -- it's part of tidyverse
df %>%
  ggplot(aes(y=var_for_y_axis, x=var_for_x_axis)) +
  geom_point()  # scatterplot
```

- `aes()` is for aesthetics -- how the chart is set up
- Other useful aesthetics:
    - `group=` to set groups to list in the legend.  Not needed if using the below though
    - `color=` to set color by some grouping variable.  Put `factor()` around the variable if you want discrete groups, otherwise it will do a color scale (light to dark)
    - `shape=` to set shapes for points -- <a href="https://cran.r-project.org/web/packages/ggplot2/vignettes/ggplot2-specs.html">see here for a list</a>

## ggplot2 for visualization {data-background="../Backgrounds/reference.jpg"}

```{r, eval=F}
library(ggplot2)  # or tidyverse -- it's part of tidyverse
df %>%
  ggplot(aes(y=var_for_y_axis, x=var_for_x_axis)) +
  geom_point()  # scatterplot
```

- `geom` stands for geometry -- any shapes, lines, etc. start with `geom`
- Other useful geoms:
    - `geom_line()`: makes a line chart
    - `geom_bar()`: makes a bar chart -- y is the height, x is the category
    - `geom_smooth(method="lm")`: Adds a linear regression into the chart
    - `geom_abline(slope=1)`: Adds a 45° line
- Add `xlab("Label text here")` to change the x-axis label
- Add `ylab("Label text here")` to change the y-axis label
- Add `ggtitle("Title text here")` to add a title
- Plenty more details in the [Data Visualization Cheat Sheet](https://github.com/rstudio/cheatsheets/blob/main/data-visualization.pdf)

## Plotting: Distribution of revenue {data-background="../Backgrounds/default.jpg"}

:::: {.columns}

::: {.column width="50%"}
1. Revenue
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_dist(train, "revtq")
```

2. Quarterly growth
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_dist(train, "revtq_gr")
```
:::

::: {.column width="50%"}
3. Year-over-year growth
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_dist(train, "revtq_yoy")
```

4. First difference
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_dist(train, "revtq_d")
```
:::

::::


## What do we learn from these graphs? {data-background="../Backgrounds/question.jpg"}

<div class="alignleft">

1. Revenue
    - $~$
    - $~$
2. Quarterly growth
    - $~$
    - $~$
3. Year-over-year growth
    - $~$
    - $~$
4. First difference
    - $~$
    - $~$
</div>

## What do we learn from these graphs? {data-background="../Backgrounds/answer.jpg" class="default exclude"}

1. Revenue
    - This is really skewed data -- a lot of small revenue quarters, but a significant amount of large revenue quarters in the tail
        - Potential fix: use `log(revtq)`?
2. Quarterly growth
    - Quarterly growth is reasonably close to normally distributed
        - Good for OLS
3. Year-over-year growth
    - Year over year growth is reasonably close to normally distributed
        - Good for OLS
4. First difference
    - Reasonably close to normally distributed, with really long tails
        - Good enough for OLS

## Plotting: Mean revenue (y) by quarter (x) {data-background="../Backgrounds/default.jpg"}

:::: {.columns}

::: {.column width="50%"}
1. Revenue
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_bar(train, "revtq")
```

2. Quarterly growth
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_bar(train, "revtq_gr")
```
:::

::: {.column width="50%"}
3. Year-over-year growth
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_bar(train, "revtq_yoy")
```

4. First difference
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_bar(train, "revtq_d")
```
:::

::::


## What do we learn from these graphs? {data-background="../Backgrounds/question.jpg"}

<div class="alignleft">

1. Revenue
    - $~$
    - $~$
2. Quarterly growth
    - $~$
    - $~$
3. Year-over-year growth
    - $~$
    - $~$
4. First difference
    - $~$
    - $~$
</div>

## What do we learn from these graphs? {data-background="../Backgrounds/answer.jpg" class="default exclude"}

1. Revenue
    - Revenue seems cyclical!
2. Quarterly growth
    - Definitely cyclical!
3. Year-over-year growth
    - Year over year difference is less cyclical -- more constant
4. First difference
    - Definitely cyclical!

## Plotting: Revenue (y) vs lag (x) by quarter {data-background="../Backgrounds/default.jpg"}

:::: {.columns}

::: {.column width="50%"}
1. Revenue
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_sct(train, "revtq", "revtq_l1")
```

2. Quarterly growth
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_sct(train, "revtq_gr", "revtq_gr1")
```
:::

::: {.column width="50%"}
3. Year-over-year growth
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_sct(train, "revtq_yoy", "revtq_yoy1")
```

4. First difference
```{r, echo=F, message=F, warning=F, fig.height=3}
plt_sct(train, "revtq_d", "revtq_d1")
```
:::

::::


## What do we learn from these graphs? {data-background="../Backgrounds/answer.jpg"}

1. Revenue
    - Revenue is really linear!  But each quarter has a distinct linear relation.
2. Quarterly growth
    - All over the place.  Each quarter appears to have a different pattern though.  Quarters will matter.
3. Year-over-year growth
    - Linear but noisy.
4. First difference
    - Again, all over the place.  Each quarter appears to have a different pattern though.  Quarters will matter.

## Correlation matrices {data-background="../Backgrounds/default.jpg"}

```{r}
cor(train[,c("revtq","revtq_l1","revtq_l2","revtq_l3", "revtq_l4")],
    use="complete.obs")
```

```{r}
cor(train[,c("revtq_gr","revtq_gr1","revtq_gr2","revtq_gr3", "revtq_gr4")],
    use="complete.obs")
```

> Retail revenue has **high** autocorrelation!  There is a concern for multicolinearity.  Revenue growth is less autocorrelated and oscillates.

## Correlation matrices {data-background="../Backgrounds/default.jpg"}

```{r}
cor(train[,c("revtq_yoy","revtq_yoy1","revtq_yoy2","revtq_yoy3", "revtq_yoy4")],
    use="complete.obs")
```

```{r}
cor(train[,c("revtq_d","revtq_d1","revtq_d2","revtq_d3", "revtq_d4")],
    use="complete.obs")
```

> Year over year change fixes the multicollinearity issue. First difference oscillates like quarter over quarter growth.

## R Practice {data-background="../Backgrounds/coding.jpg"}

- This practice will look at predicting Walmart's quarterly revenue using:
    - 1 lag
    - Cyclicality
- Practice using:
    - `lm()`
    - `{ggplot2}`
- Do the exercises in today's practice file
    - [R Practice](https://rmc.link/420r3)
    - Short link: [rmc.link/420r3](https://rmc.link/420r3)

> In this practice, we'll work through an easier variant of today's problem

# Forecasting {data-background-video="../Backgrounds/section.mp4"}

## The models {data-background="../Backgrounds/default.jpg"}

<div class="alignleft">

1. 1 Quarter lag: We saw a very strong linear pattern here earlier (*AR(1) model*)
```{r}
mod1 <- lm(revtq ~ revtq_l1, data=train)
```

2. Quarter and year lag: Year-over-year seemed pretty constant
```{r}
mod2 <- lm(revtq ~ revtq_l1 + revtq_l4, data=train)
```

3. 2 years of lags: Other lags could also help us predict (*AR(8) model*)
```{r}
mod3 <- lm(revtq ~ revtq_l1 + revtq_l2 + revtq_l3 + revtq_l4 + revtq_l5 +
           revtq_l6 + revtq_l7 + revtq_l8, data=train)
```

</div>

## Quarter lag {data-background="../Backgrounds/case.jpg"}

```{r}
summary(mod1)
```

## Quarter and year lag {data-background="../Backgrounds/case.jpg"}

```{r}
summary(mod2)
```

## 2 years of lags {data-background="../Backgrounds/case.jpg"}

```{r}
summary(mod3)
```

## Testing out of sample {data-background="../Backgrounds/reference.jpg"}

- RMSE: Root mean square Error
- RMSE is very affected by outliers, and a bad choice for noisy data where you are OK with missing a few outliers here and there
    - Doubling error *quadruples* the penalty

```{r}
rmse <- function(v1, v2) {
  sqrt(mean((v1 - v2)^2, na.rm=T))
}
```

- MAE: Mean absolute error
- MAE is measures average accuracy with no weighting
    - Doubling error *doubles* the penalty
```{r}
mae <- function(v1, v2) {
  mean(abs(v1-v2), na.rm=T)
}
```

> Both are commonly used for evaluating OLS out of sample

## Testing out of sample: Lags of revenue {data-background="../Backgrounds/case.jpg" class="small"}

```{r, echo=F, warning=F, message=F}
models <- list(mod1, mod2, mod3)
model_names <- c("1 quarter", "1 and 4 quarters", "8 quarter")

# sapply can be replaced by map here, but it would be a good bit more verbose
df_test <- data.frame(adj_r_sq=sapply(models, function(x)summary(x)[["adj.r.squared"]]),
                      rmse_in=sapply(models, function(x)rmse(train$revtq, predict(x,train))),
                      mae_in=sapply(models, function(x)mae(train$revtq, predict(x,train))),
                      rmse_out=sapply(models, function(x)rmse(test$revtq, predict(x,test))),
                      mae_out=sapply(models, function(x)mae(test$revtq, predict(x,test))))
rownames(df_test) <- model_names
html_df(df_test)  # Custom function using knitr and kableExtra
```


:::: {.columns .fragment}

::: {.column width="50%"}
1 quarter model
```{r, echo=F, warning=F, fig.height=4.5}
test %>%
  ggplot(aes(y=revtq,x=predict(mod1,test), color=factor(qtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 1 quarter model")
```
:::

::: {.column width="50%"}
8 quarter model
```{r, echo=F, warning=F, fig.height=4.5}
test %>%
  ggplot(aes(y=revtq,x=predict(mod3,test), color=factor(qtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 8 quarter model")
```
:::

::::

## What about for revenue growth? {data-background="../Backgrounds/case.jpg" class="small"}

$$
revt_t=(1+growth_t)\times revt_{t-1}
$$

```{r, echo=F, warning=F, message=F}
# models
mod1g <- lm(revtq_gr ~ revtq_gr1, data=train)
mod2g <- lm(revtq_gr ~ revtq_gr1 + revtq_gr4, data=train)
mod3g <- lm(revtq_gr ~ revtq_gr1 + revtq_gr2 + revtq_gr3 + revtq_gr4 + revtq_gr5 + revtq_gr6 + revtq_gr7 + revtq_gr8, data=train)

models <- list(mod1g, mod2g, mod3g)
model_names <- c("1 quarter", "1 and 4 quarters", "8 quarter")

df_test <- data.frame(adj_r_sq=sapply(models, function(x)summary(x)[["adj.r.squared"]]),
                      rmse_in=sapply(models, function(x)rmse(train$revtq, (1+predict(x,train))*train$revtq_l1)),
                      mae_in=sapply(models, function(x)mae(train$revtq, (1+predict(x,train))*train$revtq_l1)),
                      rmse_out=sapply(models, function(x)rmse(test$revtq, (1+predict(x,test))*test$revtq_l1)),
                      mae_out=sapply(models, function(x)mae(test$revtq, (1+predict(x,test))*test$revtq_l1)))
rownames(df_test) <- model_names
html_df(df_test)
```


:::: {.columns .fragment}

::: {.column width="50%"}
1 quarter model
```{r, echo=F, warning=F, fig.height=5}
test %>%
  ggplot(aes(y=revtq,x=(1+predict(mod1g,test))*test$revtq_l1, color=factor(qtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 1 quarter model")
```
:::

::: {.column width="50%"}
8 quarter model
```{r, echo=F, warning=F, fig.height=5}
test %>%
  ggplot(aes(y=revtq,x=(1+predict(mod3g,test))*test$revtq_l1, color=factor(qtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 8 quarter model")
```
:::

::::

## What about for YoY revenue growth? {data-background="../Backgrounds/case.jpg" class="small"}

$$
revt_t=(1+yoy\_growth_t)\times revt_{t-4}
$$

```{r, echo=F, warning=F, message=F}
# models
mod1y <- lm(revtq_yoy ~ revtq_yoy1, data=train)
mod2y <- lm(revtq_yoy ~ revtq_yoy1 + revtq_yoy4, data=train)
mod3y <- lm(revtq_yoy ~ revtq_yoy1 + revtq_yoy2 + revtq_yoy3 + revtq_yoy4 + revtq_yoy5 + revtq_yoy6 + revtq_yoy7 + revtq_yoy8, data=train)

models <- list(mod1y, mod2y, mod3y)
model_names <- c("1 quarter", "1 and 4 quarters", "8 quarter")

df_test <- data.frame(adj_r_sq=sapply(models, function(x)summary(x)[["adj.r.squared"]]),
                      rmse_in=sapply(models, function(x)rmse(train$revtq, (1+predict(x,train))*train$revtq_l4)),
                      mae_in=sapply(models, function(x)mae(train$revtq, (1+predict(x,train))*train$revtq_l4)),
                      rmse_out=sapply(models, function(x)rmse(test$revtq, (1+predict(x,test))*test$revtq_l4)),
                      mae_out=sapply(models, function(x)mae(test$revtq, (1+predict(x,test))*test$revtq_l4)))
rownames(df_test) <- model_names
html_df(df_test)
```

:::: {.columns .fragment}

::: {.column width="50%"}
1 quarter model
```{r, echo=F, warning=F, fig.height=5}
test %>%
  ggplot(aes(y=revtq,x=(1+predict(mod1y,test))*test$revtq_l4, color=factor(qtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 1 quarter model")
```
:::

::: {.column width="50%"}
8 quarter model
```{r, echo=F, warning=F, fig.height=5}
test %>%
  ggplot(aes(y=revtq,x=(1+predict(mod3y,test))*test$revtq_l4, color=factor(qtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 8 quarter model")
```
:::

::::

## What about for first difference? {data-background="../Backgrounds/case.jpg" class="small"}

$$
revt_t= change_t+ revt_{t-1}
$$

```{r, echo=F, warning=F, message=F}
# models
mod1d <- lm(revtq_d ~ revtq_d1, data=train)
mod2d <- lm(revtq_d ~ revtq_d1 + revtq_d4, data=train)
mod3d <- lm(revtq_d ~ revtq_d1 + revtq_d2 + revtq_d3 + revtq_d4 + revtq_d5 + revtq_d6 + revtq_d7 + revtq_d8, data=train)

models <- list(mod1d, mod2d, mod3d)
model_names <- c("1 quarter", "1 and 4 quarters", "8 quarter")

df_test <- data.frame(adj_r_sq=sapply(models, function(x)summary(x)[["adj.r.squared"]]),
                      rmse_in=sapply(models, function(x)rmse(train$revtq, predict(x,train)+train$revtq_l1)),
                      mae_in=sapply(models, function(x)mae(train$revtq, predict(x,train)+train$revtq_l1)),
                      rmse_out=sapply(models, function(x)rmse(test$revtq, predict(x,test)+test$revtq_l1)),
                      mae_out=sapply(models, function(x)mae(test$revtq, predict(x,test)+test$revtq_l1)))
rownames(df_test) <- model_names
html_df(df_test)
```

:::: {.columns .fragment}

::: {.column width="50%"}
1 quarter model
```{r, echo=F, warning=F, fig.height=5}
test %>%
  ggplot(aes(y=revtq,x=predict(mod1d,test)+test$revtq_l1, color=factor(qtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 1 quarter model")
```
:::

::: {.column width="50%"}
8 quarter model
```{r, echo=F, warning=F, fig.height=5}
test %>%
  ggplot(aes(y=revtq,x=predict(mod3d,test)+test$revtq_l1, color=factor(qtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 8 quarter model")
```
:::

::::

## Takeaways {data-background="../Backgrounds/reference.jpg"}

1. The first difference model works a bit better than the revenue model at predicting next quarter revenue
    - It also doesn't suffer (as much) from multicollinearity
        - This is why time series analysis is often done on first differences
            - Or second differences (difference in differences)
2. The other models perform pretty well as well
3. Extra lags generally seems helpful when accounting for cyclicality
4. Regressing by quarter helps a bit, particularly with revenue growth



# Case: Advanced revenue prediction {data-background="../Backgrounds/case.jpg"}

## RS Metrics' approach {data-background="../Backgrounds/group.jpg"}

How they do it: <a href="https://rsmetrics.com/satellite-technology-tells-us-exactly-how-well-a-business-performs-rs-metrics/">rmc.link/420class3</a>

- How does RS Metrics approach revenue prediction?
- What other creative ways might there be?

> Come up with ~3 creative approaches

- Don't worry whether the data exists or is easy to get






# End Matter {data-background-video="../Backgrounds/section.mp4"}

## Wrap up {data-background="../Backgrounds/ending.jpg"}

- For next week:
    - First individual assignment is due!
        - Finish by next class 
        - Submit on eLearn
    - Datacamp
        - Practice a bit more to keep up to date
            - Using R more will make it more natural
- Survey on the class session at this QR code:

{{< qrcode rmc.link/420survey3 qrsurvey width=200 height=200 >}}


## Packages used for these slides {data-background="../Backgrounds/code.jpg" visibility="uncounted"}

- `{broom}`
- `{DT}`
- `{downlit}`
- `{fixest}`
- `{kableExtra}`
- `{knitr}`
- `{lubridate}`
- `{plotly}`
- `{purrr}`
- `{quarto}`
- `{revealjs}`
- `{rlang}`
- `{tidyverse}`

## Custom code {data-background="../Backgrounds/code.jpg" visibility="uncounted" class="small"}

```{r, eval=F}
# Brute force code for variable generation of quarterly data lags
df <- df %>%
  group_by(gvkey) %>%
  mutate(revtq_lag1=lag(revtq), revtq_lag2=lag(revtq, 2),
         revtq_lag3=lag(revtq, 3), revtq_lag4=lag(revtq, 4),
         revtq_lag5=lag(revtq, 5), revtq_lag6=lag(revtq, 6),
         revtq_lag7=lag(revtq, 7), revtq_lag8=lag(revtq, 8),
         revtq_lag9=lag(revtq, 9), revtq_gr=revtq / revtq_lag1 - 1,
         revtq_gr1=lag(revtq_gr), revtq_gr2=lag(revtq_gr, 2),
         revtq_gr3=lag(revtq_gr, 3), revtq_gr4=lag(revtq_gr, 4),
         revtq_gr5=lag(revtq_gr, 5), revtq_gr6=lag(revtq_gr, 6),
         revtq_gr7=lag(revtq_gr, 7), revtq_gr8=lag(revtq_gr, 8),
         revtq_yoy=revtq / revtq_lag4 - 1, revtq_yoy1=lag(revtq_yoy),
         revtq_yoy2=lag(revtq_yoy, 2), revtq_yoy3=lag(revtq_yoy, 3),
         revtq_yoy4=lag(revtq_yoy, 4), revtq_yoy5=lag(revtq_yoy, 5),
         revtq_yoy6=lag(revtq_yoy, 6), revtq_yoy7=lag(revtq_yoy, 7),
         revtq_yoy8=lag(revtq_yoy, 8), revtq_d=revtq - revtq_l1,
         revtq_d1=lag(revtq_d), revtq_d2=lag(revtq_d, 2),
         revtq_d3=lag(revtq_d, 3), revtq_d4=lag(revtq_d, 4),
         revtq_d5=lag(revtq_d, 5), revtq_d6=lag(revtq_d, 6),
         revtq_d7=lag(revtq_d, 7), revtq_d8=lag(revtq_d, 8)) %>%
  ungroup()
```

## Custom code {data-background="../Backgrounds/code.jpg" visibility="uncounted" class="small"}

```{r, eval=F}
# Custom html table for small data frames
library(knitr)
library(kableExtra)
html_df <- function(text, cols=NULL, col1=FALSE, full=F) {
  if(!length(cols)) {
    cols=colnames(text)
  }
  if(!col1) {
    kable(text,"html", col.names = cols, align = c("l",rep('c',length(cols)-1))) %>%
      kable_styling(bootstrap_options = c("striped","hover"), full_width=full)
  } else {
    kable(text,"html", col.names = cols, align = c("l",rep('c',length(cols)-1))) %>%
      kable_styling(bootstrap_options = c("striped","hover"), full_width=full) %>%
      column_spec(1,bold=T)
  }
}
```


## Custom code {data-background="../Backgrounds/code.jpg" visibility="uncounted" class="small"}

```{r, eval=F}
# Density plot for 1st to 99th percentile of both columns
plt_bar <- function(df,var) {
  df %>%
    filter(eval(parse(text=var)) < quantile(eval(parse(text=var)),0.99, na.rm=TRUE),
           eval(parse(text=var)) > quantile(eval(parse(text=var)),0.01, na.rm=TRUE)) %>%
    ggplot(aes(y=eval(parse(text=var)), x=qtr)) + 
    geom_bar(stat = "summary", fun.y = "mean") + xlab(var)
}
```

```{r, eval=F}
# Scatter plot with lag for 1st to 99th percentile of data
plt_sct <- function(df,var1, var2) {
  df %>%
    filter(eval(parse(text=var1)) < quantile(eval(parse(text=var1)),0.99, na.rm=TRUE),
           eval(parse(text=var2)) < quantile(eval(parse(text=var2)),0.99, na.rm=TRUE),
           eval(parse(text=var1)) > quantile(eval(parse(text=var1)),0.01, na.rm=TRUE),
           eval(parse(text=var2)) > quantile(eval(parse(text=var2)),0.01, na.rm=TRUE)) %>%
    ggplot(aes(y=eval(parse(text=var1)), x=eval(parse(text=var2)), color=factor(qtr))) + 
    geom_point() + geom_smooth(method = "lm") + ylab(var1) + xlab(var2)
}
```

```{r, eval=F}
# Calculating various in and out of sample statistics
models <- list(mod1,mod2,mod3)
model_names <- c("1 quarter", "1 and 4 quarters", "8 quarter")

df_test <- data.frame(adj_r_sq=sapply(models, function(x)summary(x)[["adj.r.squared"]]),
                      rmse_in=sapply(models, function(x)rmse(train$revtq, predict(x,train))),
                      mae_in=sapply(models, function(x)mae(train$revtq, predict(x,train))),
                      rmse_out=sapply(models, function(x)rmse(test$revtq, predict(x,test))),
                      mae_out=sapply(models, function(x)mae(test$revtq, predict(x,test))))
rownames(df_test) <- model_names
html_df(df_test)  # Custom function using knitr and kableExtra
```

## NOTES {class="never"}

Add an illustration of fixed effects
