## A tale of three students
## (I believe I invented it: but I am not sure. "Salomon saith,
## There is no new thing upon the earth. So that as Plato had an
## imagination, that all knowledge was but remembrance; so Salomon
## giveth his sentence, that all novelty is but oblivion."
## (Francis Bacon, in Borges))
## 
## S is smart, A who is average, D who is...

## y is grade, x is effort
## data generating process (DGP) for S:
x1 <- runif(10, min=1, max=5)
y1 <- 6 + 0.8*x1 + rnorm(10, sd=0.5)

plot(x1, y1, xlim=c(0,11), ylim=c(0,11), 
     xlab="Effort", ylab="Grade")

## A:
x2 <- runif(10, min=4, max=8)
y2 <- 3 + 0.8*x2 + rnorm(10, sd=0.5)

points(x2, y2, col="red")

## D:
x3 <- runif(10, min=6, max=10)
y3 <- -2 + 0.8*x3 + rnorm(10, sd=0.5)

points(x3, y3, col="blue")

## make pooled datasets of all students, all exams
x <- c(x1,x2,x3)
y <- c(y1,y2,y3)
students <- data.frame(y=y, x=x)

## OLS linear model of grade regressed on effort
olsmod <- lm(y ~ x, students)
summary(olsmod)

abline(olsmod, lwd=2)
## nonsense coefficient! The more you study, the worse
## the grade! What happens?

## by student:
abline(lm(y~x, students[1:10, ]))
abline(lm(y~x, students[11:20, ]), col="red")
abline(lm(y~x, students[21:30, ]), col="blue")

## a look at the model specifications:

## pooled: (what I imposed in estimation)
## y_it = a + b * x_it + u_it

## individual intercepts: (how I actually generated data)
## y_it = a_i + b * x_it + u_it

## what happens when I estimate imposing pooling:
## y_it = (a + mu_i) + b * x_it + u_it
## y_it = a + b * x_it + (mu_i + u_it)

## the problem here: mu_i is correlated with the effort
## (the smart student studies less etc.)

## make a regular panel
students$id <- rep(1:3, each=10)
students$ex <- rep(1:10, 3)
## reorder columns
students <- students[ , c(3,4,1,2)]

## LSDV (least squares dummy variables)
lsdvmod <- lm(y ~ x + as.factor(id) - 1, students)

summary(lsdvmod)
## sensible coefficient now!

## consider the different means of effort and grade
## for the different students
abline(h=mean(students[students$id == 1, "y"]), lty=2)
abline(v=mean(students[students$id == 1, "x"]), lty=2)

abline(h=mean(students[students$id == 2, "y"]), lty=2, col="red")
abline(v=mean(students[students$id == 2, "x"]), lty=2, col="red")

abline(h=mean(students[students$id == 3, "y"]), lty=2, col="blue")
abline(v=mean(students[students$id == 3, "x"]), lty=2, col="blue")

## one way to "kill" unobserved heterogeneity:
## demeaning:
## wy_it = y_it - y_i. , y_i. = mean_t(y_it) 

students$dx <- students$x-rep(tapply(students$x, students$id, mean),
                     each=10)
students$dy <- students$y-rep(tapply(students$y, students$id, mean),
                     each=10)

## graphically:
## student S
plot(students$dx[1:10], students$dy[1:10],
     xlim=c(-5,5), ylim=c(-5,5),
     xlab="Effort", ylab="Grade")
abline(h=0)
abline(v=0)
## students A and D
points(students$dx[11:20], students$dy[11:20], col="red")
points(students$dx[21:30], students$dy[21:30], col="blue")

## from original model:
##  y_it = a_i + b*x_it + u_it
## to demeaned ("within"):
##  y_it-y_i. = a_i-a_i + b*(x_it-x_i.) + (u_it-u_i.)

## within model
wmod <- lm(dy ~ dx - 1, students)
summary(wmod)
## notice same coefficient as LSDV
abline(wmod, lwd=2)

## is now consistent with individual models:
abline(lm(dy~dx, students[1:10, ]))
abline(lm(dy~dx, students[11:20, ]), col="red")
abline(lm(dy~dx, students[21:30, ]), col="blue")

save(students, file="students.rda")