---
title: "Code for Session 2"
author: "Dr. Richard M. Crowley"
date: ""
output:
  html_notebook
---

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.

```{r}
library(tidyverse)
df <- read.csv("../../Data/Session_2-1.csv", stringsAsFactors=FALSE)
df_full <- df
uol <- filter(df, isin == "SG1S83002349")
#clean_df <- subset(df,fyear==2017 & !is.na(revt) & !is.na(ni) & revt > 1)
```

```{r}
# revt: Revenue, at: Assets
summary(uol[,c("revt", "at")])
```

```{r}
mod1 <- lm(revt ~ at, data = uol)
summary(mod1)
```

```{r, message=F, warning=F, fig.width=8, fig.height=1.5}
# Graph showing squared error
uolg <- uol[,c("at","revt")]
uolg$resid <- mod1$residuals
uolg$xleft <- ifelse(uolg$resid < 0,uolg$at,uolg$at - uolg$resid)
uolg$xright <- ifelse(uolg$resid < 0,uolg$at - uolg$resid, uol$at)
uolg$ytop <- ifelse(uolg$resid < 0,uolg$revt - uolg$resid,uol$revt)
uolg$ybottom <- ifelse(uolg$resid < 0,uolg$revt, uolg$revt - uolg$resid)
uolg$point <- TRUE

uolg2 <- uolg
uolg2$point <- FALSE
uolg2$at <- ifelse(uolg$resid < 0,uolg2$xright,uolg2$xleft)
uolg2$revt <- ifelse(uolg$resid < 0,uolg2$ytop,uolg2$ybottom)

uolg <- rbind(uolg, uolg2)


uolg %>% ggplot(aes(y=revt, x=at)) + 
         geom_point(aes(shape=point, group=point)) + 
         scale_shape_manual(values=c(NA,18)) + 
         geom_smooth(method="lm", se=FALSE) +
         geom_errorbarh(aes(xmax=xright, xmin = xleft)) + 
         geom_errorbar(aes(ymax=ytop, ymin = ybottom)) + 
         theme(legend.position="none") + xlim(0, 20000) + ylim(0, 2500)
```

```{r, message=F}
# tidyverse
uol <- uol %>%
  mutate(revt_growth1 = revt / lag(revt) - 1)

# R way
uol$revt_growth2 = uol$revt / c(NA, uol$revt[-length(uol$revt)]) - 1

identical(uol$revt_growth1, uol$revt_growth2)
```

```{r}
# Make the other needed change
uol <- uol %>%
  mutate(at_growth = at / lag(at) - 1) %>%  # Calculate asset growth
  rename(revt_growth = revt_growth1)        # Rename for readability
# Run the OLS model
mod2 <- lm(revt_growth ~ at_growth, data = uol)
summary(mod2)
```

```{r}
# lct: short term liabilities, che: cash and equivalents, ebit: EBIT
uol <- uol %>%
  mutate_at(vars(lct, che, ebit), list(growth = ~. / lag(.) - 1))  # Calculate 3 growths
mod3 <- lm(revt_growth ~ lct_growth + che_growth + ebit_growth, data=uol)
summary(mod3)
```

```{r}
detector <- function() {
  dice <- sample(1:6, size=2, replace=TRUE)
  if (sum(dice) == 12) {
    "exploded"
  } else {
    "still there"
  }
}

experiment <- replicate(1000,detector())
# p value
p <- sum(experiment == "still there") / 1000
if (p < 0.05) {
  paste("p-value: ", p, "-- Fail to reject H_A, sun appears to have exploded")
} else {
  paste("p-value: ", p, "-- Reject H_A that sun exploded")
}
```

```{r, warning=F, message=F}
library(tidyverse)
read_csv('../../Data/Session_3-1.csv') %>%
  ggplot(aes(y=revtq, x=atq)) + 
  geom_point() + 
  geom_smooth(method="lm") + 
  xlab("Assets") + 
  ylab("Revenue")
```

```{r, fig.height=3, fig.width=4}
plot_norm <- function(bound, outer) {
  x <- seq(-outer,outer,length=100)
  hx <- dnorm(x)
  plot(x, hx, type="n", xlab="z values", ylab="Normal PDF", main="Normal Distribution", axes=FALSE)
  lines(x, hx)
  i <- x < -bound
  polygon(c(-outer,x[i],-bound,-bound), c(0,hx[i],max(hx[i]),0), col="red")
  i <- x > bound
  polygon(c(bound, bound, x[i],outer), c(0,max(hx[i]), hx[i],0), col="red")
  axis(1, at=-outer:outer, pos=0)
}
plot_norm(1.96, 4)
```

```{r}
anova(mod2, mod3, test="Chisq")
```

```{r}
# Ensure firms have at least $1M (local currency), and have revenue
# df contains all real estate companies excluding North America
df_clean <- df %>%
  filter(at>1, revt>0)

# We cleaned out 578 observations!
print(c(nrow(df), nrow(df_clean)))

# Another useful cleaning funtion:
# Replaces NaN, Inf, and -Inf with NA for all numeric variables in the data!
df_clean <- df_clean %>%
  mutate_if(is.numeric, list(~replace(., !is.finite(.), NA)))
```

```{r, message=F, warning=F}
uol <- uol %>% mutate(revt_lead = lead(revt))  # From dplyr
forecast1 <- lm(revt_lead ~ lct + che + ebit, data=uol)
library(broom)  # Lets us view bigger regression outputs in a tidy fashion
tidy(forecast1)  # Present regression output
glance(forecast1)  # Present regression statistics
```

```{r}
forecast2 <- 
  lm(revt_lead ~ revt + act + che + lct + dp + ebit , data=uol)
tidy(forecast2)
```

```{r}
glance(forecast2)

anova(forecast1, forecast2, test="Chisq")
```

```{r}
# Note the group_by -- without it, lead() will pull from the subsequent firm!
# ungroup() tells R that we finished grouping
df_clean <- df_clean %>% 
  group_by(isin) %>% 
  mutate(revt_lead = lead(revt)) %>%
  ungroup()
```

```{r}
forecast3 <-
  lm(revt_lead ~ revt + act + che + lct + dp + ebit,
     data=df_clean[df_clean$fic=="SGP",])
tidy(forecast3)
```

```{r}
glance(forecast3)
```

```{r}
forecast4 <-
  lm(revt_lead ~ revt + act + che + lct + dp + ebit , data=df_clean)
tidy(forecast4)
```

```{r}
glance(forecast4)
```

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

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

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

```{r, message=F, fig.height=4}
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")
```

```{r}
# Exports for the following week
save(df_clean, forecast2, uol, forecast4, file = "../../Data/Session_2_export.RData")
```

