```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

#< load libraries
  library(foreign)
  library(haven)
  library(dplyr)
#>


```




# Dummies

A binary variable capturing a qualitative state in the data; e.g. female where we set the variable to 1 for women

```{r}
  wage1 <-  read.csv("https://www.dropbox.com/s/9agc2vmamfztlel/WAGE1.csv?dl=1")


  r1 <- lm(wage ~ female, wage1)
  r1 %>% summary()
  
```
How can we interpret the coefficients?

$$Wage=\beta_0+\beta_1\times FEMALE+\epsilon$$

As before: by how much we can expect the dependent variable to change if we move one unit in the explanatory variable. Here this means if we go from a man to a women (or from a women to a man). In this case: women earn on average $2.51 less.

Note that

```{r}

wage1 %>% group_by(female) %>% summarize(mean(wage))

```

Notice, that the intercept ($\beta_0$) is the average wage for men. To get the average wage for women we need to add the coefficient: 
$$ E\{wage| Women\}=\beta_0+\beta_1=4.587...$$
where $E\{wage| Women\}$ is a mathematical way of saying the average (or expected) wage for women.



```{r}
  summary(r1)
  
  summary(lm(wage ~ 0+educ, wage1))

  rows=nrow(wage1)
  wage1["one"]<-rep(1,rows)
  head(wage1) 
  
  summary(lm(wage ~ 0+educ+one, wage1))
  summary(lm(wage ~ educ+one, wage1))
  summary(lm(wage ~one, wage1))
  summary(wage1$wage)
  summary(lm(wage ~ female, wage1))
  
  aggregate(wage1$wage, list(wage1$female), mean)
  aggregate(wage1$wage, list(wage1$female), mean)
   
  wage1["male"]<-1-wage1$female
  
  summary(lm(wage ~ male +female, wage1))
  
  summary(lm(wage ~ female+male , wage1))  
  summary(lm(wage ~ 0+male +female, wage1))  
  
```

# Dummies and categorical variables


```{r}
wage1= wage1 %>% mutate(educats=0,
                        educats=ifelse(educ==12,1,educats),
                        educats=ifelse(educ>12,2,educats))

wage1 %>% group_by(educats) %>% summarize(n())

```

```{r}

summary(lm(wage~educats,wage1))

```

Making a set of dummy variables

```{r}
wage1=wage1 %>% mutate(edu_low=educats==0,
                       edu_normal=educats==1,
                       edu_high=educats==2)

reg=lm(wage~edu_low+edu_normal+edu_high,wage1)

reg %>% summary()

```



## Testing for a linear relationship
```{r,message=FALSE}
  library(car)
  
  summary(lm(wage ~ edu_low+edu_normal+edu_high, wage1))
  reg=lm(wage ~ edu_low+edu_high, wage1)
  
  linearHypothesis(reg,c("edu_lowTRUE=-edu_highTRUE"))

  
  
```

## Factor variables


```{r}
wage1 =wage1 %>% mutate(educatsf=factor(educats,label=c("low","normal","high")))
```

Compare:
```{r}
wage1 %>% group_by(educatsf) %>% summarise(n())
wage1 %>% group_by(educats) %>% summarise(n())

```

Running the same regression again:
```{r}
reg2=lm(wage ~ educatsf, wage1)
reg2 %>% summary()
```

# Dummies as dependent variable


```{r}
  bm=read.csv("https://www.dropbox.com/s/e4w6pjm6b4wz9do/bm.csv?dl=1")
  

  summary(lm(call~ofjobs+computerskills+education,bm))
  
  summary(lm(call~education,bm))
  
  
```


```{r,message=FALSE}
summary(lm(call~black,bm))
library(ggplot2)
agg=bm %>% group_by(black) %>% summarise(call=mean(call))
ggplot(agg,aes(x=black,y=call))+geom_bar(stat="identity")+theme_minimal()

```


What about continous explanatory variables?

```{r}
(linear0 <- lm(call ~ yearsexp        , data = bm)) %>% summary()
(linear1 <- lm(call ~ factor(yearsexp), data = bm)) %>% summary()

```
Student Question: A stepwise linear model?


```{r}

bm= bm %>% mutate(moret10=yearsexp>10,
                  
                  moret10Xyearsexp=moret10*yearsexp,
                  
                  
                  dev=ifelse(moret10==T, yearsexp-10,yearsexp))


lm(call ~ yearsexp        , data = bm) %>% summary()


lm(call ~ yearsexp+moret10        , data = bm) %>% summary()


res=lm(call ~ yearsexp+moret10+moret10Xyearsexp        , data = bm) 
res %>% summary()


lm(call ~ yearsexp+moret10Xyearsexp        , data = bm) %>% summary()
```


$$ y=\beta_1 + \beta_2  exper + \beta_3 MT10 + \beta_4 MT10 \times exper$$


say we are at 15 years


$$ y=\beta_1+ \beta_3 + (\beta_2 + \beta_4 )exper $$


```{r}


bm=bm %>% mutate(hat=res$coefficients[[1]] + 
                     res$coefficients[[2]] * yearsexp+
                     res$coefficients[[3]] * moret10 +
                     res$coefficients[[4]] * moret10Xyearsexp )

ggplot(bm, aes(x=yearsexp,y=hat))+geom_point()

```


```{r,message=FALSE}
agg2=agg=bm %>% group_by(yearsexp) %>% summarise(call=mean(call))
ggplot(bm,aes(x=yearsexp,y=call) ) + geom_jitter()+
           geom_smooth(method="lm")+
           geom_bar(data=agg2,aes(x=yearsexp,y=call),
                    stat="identity",color="orange",
                    fill="orange",alpha=0.3)+
          theme_minimal()

```


```{r}

  summary(lm(call~yearsexp,bm))

  summary(linear2 <- lm(call ~ black, data = bm))
  summary(linear1 <- lm(call ~ yearsexp, data = bm))
  
  library(aod)
  library(margins)
  summary(logit1 <- glm(call ~ yearsexp, data = bm, family = "binomial"))
  
  
  summary(linear1)
  margins(logit1)
  log(100)
```
Note that the linear model is fairly close to the average marginal effect


```{r}
  
  summary(logit2 <- glm(call ~ black, data = bm, family = "binomial"))
  
  
  #summary(linear2)
  margins(logit2)

  
  summary(logit3 <- glm(call ~ black+yearsexp, data = bm, family = "binomial"))
  summary(linear3<- lm(call ~ black+ yearsexp, data = bm))
  margins(logit3)
  margins(linear3)
  
  
  #test=lm(call~ yearsexp*education,bm)
  
  #mar=margins(test,type="response")
  
  plot(plogis,-5 ,5, )
  plogis(2)
  
  bm["callstar1"] = plogis(-2.7+bm$yearsexp * .3)
  bm["callstar2"] = plogis( -2.7+bm$yearsexp * 0.06)
  bm["callstar3"] = plogis( -5+bm$yearsexp * 0.06)
  
  bm["callstar4"] = plogis(-1 +bm$yearsexp*.8 )
  
  summary(bm$callstar1)
  
  bm=bm[order(bm$yearsexp),]

  bm["yearsexp_rand"]=bm$yearsexp+runif( nrow(bm),min = -0.1, max = .1)
  bm["call_rand"]=bm$call+runif( nrow(bm),min = -0.02, max = .02)
    
  plot(bm$yearsexp_rand, bm$call_rand,  main="", 
       xlab="Experience", ylab="Call", pch=1)
  
  lines(bm$yearsexp, bm$callstar1  ,col=2)
  lines(bm$yearsexp, bm$callstar2  ,col=3)
  lines(bm$yearsexp, bm$callstar3  ,col=4)
  lines(bm$yearsexp, bm$callstar4  ,col=5)
  
  
```  
  

# Nonlinear
```{r}
  wage1=wage1 %>% mutate(educ2 =educ^2) 
  summary(lm(wage ~ educ+educ2, wage1))
  
  
  library(ggplot2)
  
  ggplot(wage1, aes(x=educ, y=wage))+geom_point()+theme_minimal()+
    geom_smooth(method=lm,se=FALSE,formula=y ~ poly(x, 2))
  
```  
`

# log linear
```{r}
  library(ggplot2)
  z=seq(-.1,.1,.001)
  ln1Pz=log(1+z)
  df=data.frame(z,ln1Pz)
  ggplot(df,aes(x=z))+geom_line(aes(y=ln1Pz),color="green")+
            geom_line(aes(y=z),color="red") + theme_minimal() + ylab("ln(1+z)")

```


# Many dummies as control
```{r}
library(dplyr)
prod=read.csv("https://www.dropbox.com/s/v1j4xzkaado8zlz/prod_balanced.csv?dl=1")
prod=prod %>% mutate(lnl=log(l),lnva=log(va))
prod=prod %>% arrange(id,year) %>% select(id,year,va,l,k,lnk,lnl,lnva)


reg1<-lm(lnva~lnl+lnk,prod) 
reg1 %>% summary()

reg2<-lm(lnva~lnl+lnk+factor(year),prod) 
reg2 %>% summary()

reg2<-lm(lnva~lnl+lnk+factor(year)+factor(id),prod) 
reg2 %>% summary()


```


# Heteroskedasticity
```{r,eval=FALSE}
  # load necessary packages for importing the function
  library(RCurl)
  
  # import the function from repository
  url_robust <- "https://raw.githubusercontent.com/IsidoreBeautrelet/economictheoryblog/master/robust_summary.R"
  eval(parse(text = getURL(url_robust, ssl.verifypeer = FALSE)),
       envir=.GlobalEnv)
  
  r_hetero<-lm(wage ~ educatsf*exper, wage1)
  summary(r_hetero)  
  r_hetero<-lm(wage ~ educatsf*exper, wage1)
  summary(r_hetero, robust=T)  
  
  
```

