Code for Session 3

Note that the directories used to store data are likely different on your computer, and such references will need to be changed before using any such code.

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)
  }
}

library(tidyverse)
df <- read.csv("../../Data/Session_3-1.csv", stringsAsFactors=FALSE)
wmt <- filter(df, tic == "WMT")
# load in relevant data from Session 2
load("../../Data/Session_2_export.RData")

expectations <- read_csv("../../Data/general-business-expectations-by-detailed-services-industry-quarterly.csv") %>%
  mutate(year = as.numeric(substr(quarter, 1, 4))) %>%    # split out year
  mutate(quarter = as.numeric(substr(quarter, 7, 7))) %>% # split out quarter
  mutate(value = as.numeric(value))                       # Ensue value is numeric

Parsed with column specification:
cols(
  quarter = col_character(),
  level_1 = col_character(),
  level_2 = col_character(),
  level_3 = col_character(),
  value = col_character()
)
NAs introduced by coercion

# extract out Q1, finance only
expectations_avg <- expectations %>%
  filter(quarter == 1,                               # Keep only the first quarter
         level_2 == "Financial & Insurance") %>%     # Keep only financial responses
  group_by(year) %>%                                 # Group data by year
  mutate(fin_sentiment=mean(value, na.rm=TRUE)) %>%  # Calculate average
  slice(1)                                           # Take only 1 row per group

library(DT)

expectations %>%
  arrange(level_2, level_3, desc(year)) %>%  # sort the data
  select(year, quarter, level_2, level_3, value) %>%  # keep only these variables
  datatable(options = list(pageLength = 5), rownames=FALSE)  # display using DT

# subset out our Singaporean data, since our macro data is Singapore-specific
df_SG <- df_clean %>% filter(fic == "SGP")
# Create year in df_SG (date is given by datadate as YYYYMMDD)
df_SG$year = round(df_SG$datadate / 10000, digits=0)
# Combine datasets
# Notice how it automatically figures out to join by "year"
df_SG_macro <- left_join(df_SG, expectations_avg[,c("year","fin_sentiment")])

Joining, by = "year"

macro1 <- lm(revt_lead ~ revt + act + che + lct + dp + ebit + fin_sentiment,
             data=df_SG_macro)
library(broom)
tidy(macro1)

df_SG_macro %>%
  ggplot(aes(y=revt_lead,
             x=fin_sentiment)) + 
  geom_point()

df_SG_macro %>%
  ggplot(aes(y=revt_lead,
    x=scale(fin_sentiment) * revt)) + 
  geom_point()

# Scale creates z-scores, but returns a matrix by default.  [,1] gives a vector
df_SG_macro$fin_sent_scaled <- scale(df_SG_macro$fin_sentiment)[,1]
macro3 <-
  lm(revt_lead ~ revt + act + che + lct + dp + ebit + fin_sent_scaled:revt,
     data=df_SG_macro)
tidy(macro3)

glance(macro3)

baseline <-
  lm(revt_lead ~ revt + act + che + lct + dp + ebit,
     data=df_SG_macro[!is.na(df_SG_macro$fin_sentiment),])
glance(baseline)

glance(macro3)

anova(baseline, macro3, test="Chisq")

Analysis of Variance Table

Model 1: revt_lead ~ revt + act + che + lct + dp + ebit
Model 2: revt_lead ~ revt + act + che + lct + dp + ebit + fin_sent_scaled:revt
  Res.Df      RSS Df Sum of Sq Pr(>Chi)   
1    304 14285622                         
2    303 13949301  1    336321 0.006875 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

r_sd <- round(sd(df_SG_macro$fin_sentiment, na.rm=T),1)
r_min <- min(df_SG_macro$fin_sentiment, na.rm=T)
r_max <- max(df_SG_macro$fin_sentiment, na.rm=T)
rev <- macro3$coefficients[["revt:fin_sent_scaled"]]
r_rev = round(100 * rev,0)
rev_min <- round((r_min / r_sd * rev + rev)*100,1)
rev_max <- round((r_max / r_sd * rev + rev)*100,1)

p_uol <- predict(forecast2, uol[uol$fyear==2017,])
p_base <- predict(baseline,
  df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear==2017,])
p_macro <- predict(macro3,
  df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear==2017,])
p_world <- predict(forecast4,
  df_clean[df_clean$isin=="SG1S83002349" & df_clean$fyear==2017,])
preds <- c(p_uol, p_base, p_macro, p_world)
names(preds) <- c("UOL 2018 UOL", "UOL 2018 Base", "UOL 2018 Macro",
                  "UOL 2018 World")
preds

  UOL 2018 UOL  UOL 2018 Base UOL 2018 Macro UOL 2018 World 
      3177.073       2086.437       2024.842       2589.636 

library(plotly)
df_SG_macro$pred_base <- predict(baseline, df_SG_macro)
df_SG_macro$pred_macro <- predict(macro3, df_SG_macro)
df_clean$pred_world <- predict(forecast4, df_clean)
uol$pred_uol <- predict(forecast2, uol)
df_preds <- data.frame(preds=preds, fyear=c(2018,2018,2018,2018), model=c("UOL only", "Base", "Macro", "World"))
plot <- ggplot() + 
  geom_point(data=df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear < 2017,], aes(y=revt_lead,x=fyear, color="Actual")) +
  geom_line(data=df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear < 2017,], aes(y=revt_lead,x=fyear, color="Actual")) + 
  geom_point(data=uol[uol$fyear < 2017,], aes(y=pred_uol,x=fyear, color="UOL only")) +
  geom_line(data=uol[uol$fyear < 2017,], aes(y=pred_uol,x=fyear, color="UOL only")) +
  geom_point(data=df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear < 2017,], aes(y=pred_base,x=fyear, color="Base")) +
  geom_line(data=df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear < 2017,], aes(y=pred_base,x=fyear, color="Base")) +
  geom_point(data=df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear < 2017,], aes(y=pred_macro,x=fyear, color="Macro")) +
  geom_line(data=df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear < 2017,], aes(y=pred_macro,x=fyear, color="Macro")) + 
  geom_point(data=df_clean[df_clean$isin=="SG1S83002349" & df_clean$fyear < 2017,], aes(y=pred_world,x=fyear, color="World")) +
  geom_line(data=df_clean[df_clean$isin=="SG1S83002349" & df_clean$fyear < 2017,], aes(y=pred_world,x=fyear, color="World")) + 
  geom_point(data=df_preds, aes(y=preds, x=fyear, color=model), size=1.5, shape=18)
ggplotly(plot)

actual_series <- df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear < 2017,]$revt_lead
uol_series <- uol[uol$fyear < 2017,]$pred_uol
base_series <- df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear < 2017,]$pred_base
macro_series <- df_SG_macro[df_SG_macro$isin=="SG1S83002349" & df_SG_macro$fyear < 2017,]$pred_macro
world_series <- df_clean[df_clean$isin=="SG1S83002349" & df_clean$fyear < 2017,]$pred_world

# series vectors calculated here -- See appendix
rmse <- function(v1, v2) {
  sqrt(mean((v1 - v2)^2, na.rm=T))
}
rmse <- c(rmse(actual_series, uol_series), rmse(actual_series, base_series),
          rmse(actual_series, macro_series), rmse(actual_series, world_series))
names(rmse) <- c("UOL 2018 UOL", "UOL 2018 Base", "UOL 2018 Macro", "UOL 2018 World")
rmse

  UOL 2018 UOL  UOL 2018 Base UOL 2018 Macro UOL 2018 World 
      175.5609       301.3161       344.9681       332.8101 

preds

  UOL 2018 UOL  UOL 2018 Base UOL 2018 Macro UOL 2018 World 
      3177.073       2086.437       2024.842       2589.636 

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 YYMMDDs10: YYYY/MM/DD, can be converted from text to date easily
df$date <- as.Date(df$datadate)  # Built in to R

html_df(head(df[,c("conm","date","revtq","revtq_gr", "revtq_yoy", "revtq_d")]))

conm	date	revtq	revtq_gr	revtq_yoy	revtq_d
ALLIED STORES	1962-04-30	156.5	NA	NA	NA
ALLIED STORES	1962-07-31	161.9	0.0345048	NA	5.4
ALLIED STORES	1962-10-31	176.9	0.0926498	NA	15.0
ALLIED STORES	1963-01-31	275.5	0.5573770	NA	98.6
ALLIED STORES	1963-04-30	171.1	-0.3789474	0.0932907	-104.4
ALLIED STORES	1963-07-31	182.2	0.0648743	0.1253860	11.1

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

# Custom Function to generate a series of lags
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(mutate(group_by(df, gvkey), !!!quosures))
}
# Generate lags "revtq_l#"
df <- multi_lag(df, 1:8, revtq, "_l")
# Generate changes "revtq_gr#"
df <- multi_lag(df, 1:8, revtq_gr)
# Generate year-over-year changes "revtq_yoy#"
df <- multi_lag(df, 1:8, revtq_yoy)
# Generate first differences "revtq_d#"
df <- multi_lag(df, 1:8, revtq_d)
# Equivalent brute force code for this is in the appendix

html_df(head(df[,c("conm","date","revtq","revtq_l1", "revtq_l2", "revtq_l3", "revtq_l4")]))

conm	date	revtq	revtq_l1	revtq_l2	revtq_l3	revtq_l4
ALLIED STORES	1962-04-30	156.5	NA	NA	NA	NA
ALLIED STORES	1962-07-31	161.9	156.5	NA	NA	NA
ALLIED STORES	1962-10-31	176.9	161.9	156.5	NA	NA
ALLIED STORES	1963-01-31	275.5	176.9	161.9	156.5	NA
ALLIED STORES	1963-04-30	171.1	275.5	176.9	161.9	156.5
ALLIED STORES	1963-07-31	182.2	171.1	275.5	176.9	161.9

# Clean the data: Replace NaN, Inf, and -Inf with NA
df <- df %>%
  mutate_if(is.numeric, list(~replace(., !is.finite(.), NA)))

`mutate_if()` ignored the following grouping variables:
Column `gvkey`

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

summary(df[,c("revtq","revtq_gr","revtq_yoy", "revtq_d","fqtr")])

     revtq              revtq_gr         revtq_yoy          revtq_d               fqtr      
 Min.   :     0.00   Min.   :-1.0000   Min.   :-1.0000   Min.   :-24325.21   Min.   :1.000  
 1st Qu.:    64.46   1st Qu.:-0.1112   1st Qu.: 0.0077   1st Qu.:   -19.33   1st Qu.:1.000  
 Median :   273.95   Median : 0.0505   Median : 0.0740   Median :     4.30   Median :2.000  
 Mean   :  2439.38   Mean   : 0.0650   Mean   : 0.1273   Mean   :    22.66   Mean   :2.478  
 3rd Qu.:  1254.21   3rd Qu.: 0.2054   3rd Qu.: 0.1534   3rd Qu.:    55.02   3rd Qu.:3.000  
 Max.   :136267.00   Max.   :14.3333   Max.   :47.6600   Max.   : 15495.00   Max.   :4.000  
 NA's   :367         NA's   :690       NA's   :940       NA's   :663                        

# 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)
}

# 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=fqtr)) + 
    geom_bar(stat = "summary", fun.y = "mean") + xlab(var)
}

# 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(fqtr))) + 
    geom_point() + geom_smooth(method = "lm") + ylab(var1) + xlab(var2)
}

plt_dist(train, "revtq")

plt_dist(train, "revtq_gr")

plt_dist(train, "revtq_yoy")

plt_dist(train, "revtq_d")

plt_bar(train, "revtq")

plt_bar(train, "revtq_gr")

plt_bar(train, "revtq_yoy")

plt_bar(train, "revtq_d")

plt_sct(train, "revtq", "revtq_l1")

plt_sct(train, "revtq_gr", "revtq_gr1")

plt_sct(train, "revtq_yoy", "revtq_yoy1")

plt_sct(train, "revtq_d", "revtq_d1")

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

             revtq  revtq_l1  revtq_l2  revtq_l3  revtq_l4
revtq    1.0000000 0.9916167 0.9938489 0.9905522 0.9972735
revtq_l1 0.9916167 1.0000000 0.9914767 0.9936977 0.9898184
revtq_l2 0.9938489 0.9914767 1.0000000 0.9913489 0.9930152
revtq_l3 0.9905522 0.9936977 0.9913489 1.0000000 0.9906006
revtq_l4 0.9972735 0.9898184 0.9930152 0.9906006 1.0000000

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

             revtq_gr   revtq_gr1   revtq_gr2   revtq_gr3   revtq_gr4
revtq_gr   1.00000000 -0.32291329  0.06299605 -0.22769442  0.64570015
revtq_gr1 -0.32291329  1.00000000 -0.31885020  0.06146805 -0.21923630
revtq_gr2  0.06299605 -0.31885020  1.00000000 -0.32795121  0.06775742
revtq_gr3 -0.22769442  0.06146805 -0.32795121  1.00000000 -0.31831023
revtq_gr4  0.64570015 -0.21923630  0.06775742 -0.31831023  1.00000000

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

           revtq_yoy revtq_yoy1 revtq_yoy2 revtq_yoy3 revtq_yoy4
revtq_yoy  1.0000000  0.6554179  0.4127263  0.4196003  0.1760055
revtq_yoy1 0.6554179  1.0000000  0.5751128  0.3665961  0.3515105
revtq_yoy2 0.4127263  0.5751128  1.0000000  0.5875643  0.3683539
revtq_yoy3 0.4196003  0.3665961  0.5875643  1.0000000  0.5668211
revtq_yoy4 0.1760055  0.3515105  0.3683539  0.5668211  1.0000000

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

            revtq_d   revtq_d1   revtq_d2   revtq_d3   revtq_d4
revtq_d   1.0000000 -0.6181516  0.3309349 -0.6046998  0.9119911
revtq_d1 -0.6181516  1.0000000 -0.6155259  0.3343317 -0.5849841
revtq_d2  0.3309349 -0.6155259  1.0000000 -0.6191366  0.3165450
revtq_d3 -0.6046998  0.3343317 -0.6191366  1.0000000 -0.5864285
revtq_d4  0.9119911 -0.5849841  0.3165450 -0.5864285  1.0000000

mod1 <- lm(revtq ~ revtq_l1, data=train)

mod2 <- lm(revtq ~ revtq_l1 + revtq_l4, data=train)

mod3 <- lm(revtq ~ revtq_l1 + revtq_l2 + revtq_l3 + revtq_l4 + 
             revtq_l5 + revtq_l6 + revtq_l7 + revtq_l8, data=train)

mod4 <- lm(revtq ~ (revtq_l1 + revtq_l2 + revtq_l3 + revtq_l4 +
             revtq_l5 + revtq_l6 + revtq_l7 + revtq_l8):factor(fqtr),
           data=train)

summary(mod1)

Call:
lm(formula = revtq ~ revtq_l1, data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-24438.7    -34.0    -11.7     34.6  15200.5 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 15.639975  13.514877   1.157    0.247    
revtq_l1     1.003038   0.001556 644.462   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1152 on 7676 degrees of freedom
  (662 observations deleted due to missingness)
Multiple R-squared:  0.9819,    Adjusted R-squared:  0.9819 
F-statistic: 4.153e+05 on 1 and 7676 DF,  p-value: < 2.2e-16

summary(mod2)

Call:
lm(formula = revtq ~ revtq_l1 + revtq_l4, data = train)

Residuals:
     Min       1Q   Median       3Q      Max 
-20245.7    -18.4     -4.4     19.1   9120.8 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 5.444986   7.145633   0.762    0.446    
revtq_l1    0.231759   0.005610  41.312   <2e-16 ***
revtq_l4    0.815570   0.005858 139.227   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 592.1 on 7274 degrees of freedom
  (1063 observations deleted due to missingness)
Multiple R-squared:  0.9954,    Adjusted R-squared:  0.9954 
F-statistic: 7.94e+05 on 2 and 7274 DF,  p-value: < 2.2e-16

summary(mod3)

Call:
lm(formula = revtq ~ revtq_l1 + revtq_l2 + revtq_l3 + revtq_l4 + 
    revtq_l5 + revtq_l6 + revtq_l7 + revtq_l8, data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-5005.6   -12.9    -3.7     9.3  5876.3 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.02478    4.37003   0.921   0.3571    
revtq_l1     0.77379    0.01229  62.972  < 2e-16 ***
revtq_l2     0.10497    0.01565   6.707 2.16e-11 ***
revtq_l3    -0.03091    0.01538  -2.010   0.0445 *  
revtq_l4     0.91982    0.01213  75.800  < 2e-16 ***
revtq_l5    -0.76459    0.01324 -57.749  < 2e-16 ***
revtq_l6    -0.08080    0.01634  -4.945 7.80e-07 ***
revtq_l7     0.01146    0.01594   0.719   0.4721    
revtq_l8     0.07924    0.01209   6.554 6.03e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 346 on 6666 degrees of freedom
  (1665 observations deleted due to missingness)
Multiple R-squared:  0.9986,    Adjusted R-squared:  0.9986 
F-statistic: 5.802e+05 on 8 and 6666 DF,  p-value: < 2.2e-16

summary(mod4)

Call:
lm(formula = revtq ~ (revtq_l1 + revtq_l2 + revtq_l3 + revtq_l4 + 
    revtq_l5 + revtq_l6 + revtq_l7 + revtq_l8):factor(fqtr), 
    data = train)

Residuals:
    Min      1Q  Median      3Q     Max 
-6066.6   -13.9     0.1    15.1  4941.1 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -0.201107   4.004046  -0.050 0.959944    
revtq_l1:factor(fqtr)1  0.488584   0.021734  22.480  < 2e-16 ***
revtq_l1:factor(fqtr)2  1.130563   0.023017  49.120  < 2e-16 ***
revtq_l1:factor(fqtr)3  0.774983   0.028727  26.977  < 2e-16 ***
revtq_l1:factor(fqtr)4  0.977353   0.026888  36.349  < 2e-16 ***
revtq_l2:factor(fqtr)1  0.258024   0.035136   7.344 2.33e-13 ***
revtq_l2:factor(fqtr)2 -0.100284   0.024664  -4.066 4.84e-05 ***
revtq_l2:factor(fqtr)3  0.212954   0.039698   5.364 8.40e-08 ***
revtq_l2:factor(fqtr)4  0.266761   0.035226   7.573 4.14e-14 ***
revtq_l3:factor(fqtr)1  0.124187   0.036695   3.384 0.000718 ***
revtq_l3:factor(fqtr)2 -0.042214   0.035787  -1.180 0.238197    
revtq_l3:factor(fqtr)3 -0.005758   0.024367  -0.236 0.813194    
revtq_l3:factor(fqtr)4 -0.308661   0.038974  -7.920 2.77e-15 ***
revtq_l4:factor(fqtr)1  0.459768   0.038266  12.015  < 2e-16 ***
revtq_l4:factor(fqtr)2  0.684943   0.033366  20.528  < 2e-16 ***
revtq_l4:factor(fqtr)3  0.252169   0.035708   7.062 1.81e-12 ***
revtq_l4:factor(fqtr)4  0.817136   0.017927  45.582  < 2e-16 ***
revtq_l5:factor(fqtr)1 -0.435406   0.023278 -18.704  < 2e-16 ***
revtq_l5:factor(fqtr)2 -0.725000   0.035497 -20.424  < 2e-16 ***
revtq_l5:factor(fqtr)3 -0.160408   0.036733  -4.367 1.28e-05 ***
revtq_l5:factor(fqtr)4 -0.473030   0.033349 -14.184  < 2e-16 ***
revtq_l6:factor(fqtr)1  0.059832   0.034672   1.726 0.084453 .  
revtq_l6:factor(fqtr)2  0.154990   0.025368   6.110 1.05e-09 ***
revtq_l6:factor(fqtr)3 -0.156840   0.041147  -3.812 0.000139 ***
revtq_l6:factor(fqtr)4 -0.106082   0.037368  -2.839 0.004541 ** 
revtq_l7:factor(fqtr)1  0.060031   0.038599   1.555 0.119936    
revtq_l7:factor(fqtr)2  0.061381   0.034510   1.779 0.075344 .  
revtq_l7:factor(fqtr)3  0.028149   0.025385   1.109 0.267524    
revtq_l7:factor(fqtr)4 -0.277337   0.039380  -7.043 2.08e-12 ***
revtq_l8:factor(fqtr)1 -0.016637   0.033568  -0.496 0.620177    
revtq_l8:factor(fqtr)2 -0.152379   0.028014  -5.439 5.54e-08 ***
revtq_l8:factor(fqtr)3  0.052208   0.027334   1.910 0.056179 .  
revtq_l8:factor(fqtr)4  0.103495   0.015777   6.560 5.78e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 299.7 on 6642 degrees of freedom
  (1665 observations deleted due to missingness)
Multiple R-squared:  0.9989,    Adjusted R-squared:  0.9989 
F-statistic: 1.935e+05 on 32 and 6642 DF,  p-value: < 2.2e-16

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

mae <- function(v1, v2) {
  mean(abs(v1-v2), na.rm=T)
}

models <- list(mod1,mod2,mod3,mod4)
model_names <- c("1 period", "1 and 4 periods", "8 periods", "8 periods w/ quarters")
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

	adj_r_sq	rmse_in	mae_in	rmse_out	mae_out
1 period	0.9818514	1151.3535	322.73819	2947.3619	1252.5196
1 and 4 periods	0.9954393	591.9500	156.20811	1400.3841	643.9823
8 periods	0.9985643	345.8053	94.91083	677.6218	340.8236
8 periods w/ quarters	0.9989231	298.9557	91.28056	645.5415	324.9395

test %>%
  ggplot(aes(y=revtq,x=predict(mod1,test), color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 1 period model")

test %>%
  ggplot(aes(y=revtq,x=predict(mod4,test), color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 8 period X quarter model")

# 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)
mod4g <- lm(revtq_gr ~ (revtq_gr1 + revtq_gr2 + revtq_gr3 + revtq_gr4 + revtq_gr5 + revtq_gr6 + revtq_gr7 + revtq_gr8):factor(fqtr), data=train)
models <- list(mod1g, mod2g, mod3g, mod4g)
model_names <- c("1 period", "1 and 4 periods", "8 periods", "8 periods w/ quarters")
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)

	adj_r_sq	rmse_in	mae_in	rmse_out	mae_out
1 period	0.0910390	1106.3730	308.48331	3374.728	1397.6541
1 and 4 periods	0.4398456	530.6444	154.15086	1447.035	679.3536
8 periods	0.6761666	456.2551	123.34075	1254.201	584.9709
8 periods w/ quarters	0.7758834	378.4082	98.45751	1015.971	436.1522

test %>%
  ggplot(aes(y=revtq,x=(1+predict(mod1g,test))*test$revtq_l1, color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 1 period model")

test %>%
  ggplot(aes(y=revtq,x=(1+predict(mod4g,test))*test$revtq_l1, color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 8 period X quarter model")

# 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)
mod4y <- lm(revtq_gr ~ (revtq_yoy1 + revtq_yoy2 + revtq_yoy3 + revtq_yoy4 + revtq_yoy5 + revtq_yoy6 + revtq_yoy7 + revtq_yoy8):factor(fqtr), data=train)
models <- list(mod1y, mod2y, mod3y, mod4y)
model_names <- c("1 period", "1 and 4 periods", "8 periods", "8 periods w/ quarters")
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)

	adj_r_sq	rmse_in	mae_in	rmse_out	mae_out
1 period	0.4370372	513.3264	129.2309	1867.4957	798.0327
1 and 4 periods	0.5392281	487.6441	126.6012	1677.4003	731.2841
8 periods	0.5398870	384.2923	101.0104	822.0065	403.5445
8 periods w/ quarters	0.1563169	714.4285	195.3204	1231.8436	617.2989

test %>%
  ggplot(aes(y=revtq,x=(1+predict(mod1y,test))*test$revtq_l4, color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 1 period model")

test %>%
  ggplot(aes(y=revtq,x=(1+predict(mod3y,test))*test$revtq_l4, color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 8 period model")

# 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)
mod4d <- lm(revtq_d ~ (revtq_d1 + revtq_d2 + revtq_d3 + revtq_d4 + revtq_d5 + revtq_d6 + revtq_d7 + revtq_d8):factor(fqtr), data=train)
models <- list(mod1d, mod2d, mod3d, mod4d)
model_names <- c("1 period", "1 and 4 periods", "8 periods", "8 periods w/ quarters")
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)

	adj_r_sq	rmse_in	mae_in	rmse_out	mae_out
1 period	0.3532044	896.7969	287.77940	2252.7605	1022.0960
1 and 4 periods	0.8425348	454.8651	115.52694	734.8120	377.5281
8 periods	0.9220849	333.0054	95.95924	651.4967	320.0567
8 periods w/ quarters	0.9397434	292.3102	86.95563	659.4412	319.7305

test %>%
  ggplot(aes(y=revtq,x=predict(mod1d,test)+test$revtq_l1, color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 1 period model")

test %>%
  ggplot(aes(y=revtq,x=predict(mod4d,test)+test$revtq_l1, color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue") + 
  xlab("Prediction: 8 period X quarter model")

# 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)
mod4g <- lm(revtq_gr ~ (revtq_gr1 + revtq_gr2 + revtq_gr3 + revtq_gr4 + revtq_gr5 + revtq_gr6 + revtq_gr7 + revtq_gr8):factor(fqtr), data=train)
models <- list(mod1g, mod2g, mod3g, mod4g)
model_names <- c("1 period", "1 and 4 periods", "8 periods", "8 periods w/ quarters")
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_gr, predict(x,train))),
                      mae_in=sapply(models, function(x)mae(train$revtq_gr, predict(x,train))),
                      rmse_out=sapply(models, function(x)rmse(test$revtq_gr, predict(x,test))),
                      mae_out=sapply(models, function(x)mae(test$revtq_gr, predict(x,test))))
rownames(df_test) <- model_names
html_df(df_test)

	adj_r_sq	rmse_in	mae_in	rmse_out	mae_out
1 period	0.0910390	0.3509269	0.2105219	0.2257396	0.1750580
1 and 4 periods	0.4398456	0.2681899	0.1132003	0.1597771	0.0998087
8 periods	0.6761666	0.1761825	0.0867347	0.1545298	0.0845826
8 periods w/ quarters	0.7758834	0.1462979	0.0765792	0.1459460	0.0703554

test %>%
  ggplot(aes(y=revtq_gr,x=predict(mod1g,test), color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue growth") + 
  xlab("Prediction: 1 period model")

test %>%
  ggplot(aes(y=revtq_gr,x=predict(mod4g,test), color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue growth") + 
  xlab("Prediction: 8 period X quarter model")

# 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)
mod4y <- lm(revtq_gr ~ (revtq_yoy1 + revtq_yoy2 + revtq_yoy3 + revtq_yoy4 + revtq_yoy5 + revtq_yoy6 + revtq_yoy7 + revtq_yoy8):factor(fqtr), data=train)
models <- list(mod1y, mod2y, mod3y, mod4y)
model_names <- c("1 period", "1 and 4 periods", "8 periods", "8 periods w/ quarters")
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_yoy, predict(x,train))),
                      mae_in=sapply(models, function(x)mae(train$revtq_yoy, predict(x,train))),
                      rmse_out=sapply(models, function(x)rmse(test$revtq_yoy, predict(x,test))),
                      mae_out=sapply(models, function(x)mae(test$revtq_yoy, predict(x,test))))
rownames(df_test) <- model_names
html_df(df_test)

	adj_r_sq	rmse_in	mae_in	rmse_out	mae_out
1 period	0.4370372	0.3116645	0.1114610	0.1515638	0.0942544
1 and 4 periods	0.5392281	0.2451749	0.1015699	0.1498755	0.0896079
8 periods	0.5398870	0.1928940	0.0764447	0.1346238	0.0658011
8 periods w/ quarters	0.1563169	0.3006075	0.1402156	0.1841025	0.0963205

test %>%
  ggplot(aes(y=revtq_yoy,x=predict(mod1y,test), color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual year over year revenue growth") + 
  xlab("Prediction: 1 period model")

test %>%
  ggplot(aes(y=revtq_yoy,x=predict(mod3y,test), color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual year over year revenue growth") + 
  xlab("Prediction: 8 period model")

# 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)
mod4d <- lm(revtq_d ~ (revtq_d1 + revtq_d2 + revtq_d3 + revtq_d4 + revtq_d5 + revtq_d6 + revtq_d7 + revtq_d8):factor(fqtr), data=train)
models <- list(mod1d, mod2d, mod3d, mod4d)
model_names <- c("1 period", "1 and 4 periods", "8 periods", "8 periods w/ quarters")
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_d, predict(x,train))),
                      mae_in=sapply(models, function(x)mae(train$revtq_d, predict(x,train))),
                      rmse_out=sapply(models, function(x)rmse(test$revtq_d, predict(x,test))),
                      mae_out=sapply(models, function(x)mae(test$revtq_d, predict(x,test))))
rownames(df_test) <- model_names
html_df(df_test)

	adj_r_sq	rmse_in	mae_in	rmse_out	mae_out
1 period	0.3532044	896.7969	287.77940	2252.7605	1022.0960
1 and 4 periods	0.8425348	454.8651	115.52694	734.8120	377.5281
8 periods	0.9220849	333.0054	95.95924	651.4967	320.0567
8 periods w/ quarters	0.9397434	292.3102	86.95563	659.4412	319.7305

test %>%
  ggplot(aes(y=revtq_d,x=predict(mod1d,test), color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue first difference") + 
  xlab("Prediction: 1 period model")

test %>%
  ggplot(aes(y=revtq_d,x=predict(mod4d,test), color=factor(fqtr))) +
  geom_abline(slope=1) + geom_point() +
  ylab("Actual revenue first difference") + 
  xlab("Prediction: 8 period X quarter model")