setwd("/Users/hedibert/Dropbox/Econometria16-2/aula5-2408-suposicoes")

data = read.table("http://hedibert.org/wp-content/uploads/2016/08/100metrosrasos.txt",header=TRUE)

attach(data)

n = nrow(data)


par(mfrow=c(1,1))
plot(ano,homen,ylim=range(homen,mulher),pch=16,ylab="segundos")
points(ano,mulher,col=2,pch=16)
legend("topright",legend=c("homen","mulher"),col=1:2,pch=16)
text(1987.5,10.4,"Florence Delorez Griffith Joyner",col=2)
text(2004,9.6,"Usain Bolt")

# Homens
pdf(file="100metrosrasos.pdf",width=10,height=8)


par(mfrow=c(1,1))
plot(ano,homen,pch=16,ylab="segundos",main="100 metros rasos para homens")

m.homen = mean(homen)
m.ano = mean(ano)

b1 = sum((homen-m.homen)*(ano-m.ano))/sum((ano-m.ano)^2)
b0 = m.homen-b1*m.ano
c(b0,b1)

SSR = sum((homen-b0-b1*ano)^2)

sig2hat = SSR/(n-2)

SST = sum((ano-m.ano)^2)
varb0 = sig2hat*(1/n+m.ano^2/SST)

varb1 = sig2hat/SST


c(b0-2*sqrt(varb0),b0+2*sqrt(varb0))

c(b1-2*sqrt(varb1),b1+2*sqrt(varb1))


covb0b1 = -m.ano*varb1

corrb0b1 = covb0b1/sqrt(varb0*varb1)

corrb0b1

Lb0 = b0+qt(0.025,n-2)*sqrt(varb0)
Lb1 = b1+qt(0.025,n-2)*sqrt(varb1)
Ub0 = b0+qt(0.975,n-2)*sqrt(varb0)
Ub1 = b1+qt(0.975,n-2)*sqrt(varb1)

c(Lb0,Ub0)
c(Lb1,Ub1)

par(mfrow=c(1,1))
plot(ano,homen,pch=16,ylab="segundos",main="100 metros rasos para homens")
abline(b0,b1,col=2,lwd=2)
text(2000,10.1,paste("b0 = ",round(b0,3),"  -  s.e.(b0) = ",round(sqrt(varb0),3),sep=""))
text(2000,10.05,paste("b1 = ",round(b1,7),"  -  s.e.(b1) = ",round(sqrt(varb1),7),sep=""))
text(2000,10,paste("cor(b0,b1) = ",round(corrb0b1,7),sep=""))
text(1985,9.65,paste("IC 95% b0 : (",round(Lb0,3),",",round(Ub0,3),")",sep=""))
text(1985,9.6,paste("IC 95% b1 : (",round(Lb1,7),",",round(Ub1,7),")",sep=""))

b1s = rnorm(1000,b1,sqrt(varb1))
b0s = m.homen-b1s*m.ano
for (i in 1:100)
  lines(ano,b0s[i]+b1s[i]*ano,col=grey(0.9))
abline(b0,b1,col=2,lwd=2)
points(ano,homen,pch=16)


par(mfrow=c(1,1))
plot(ano,homen,pch=16,ylab="segundos",main="100 metros rasos para homens")
abline(b0,b1,col=2,lwd=2)
text(1985,9.65,paste("IC 95% b0 : (",round(Lb0,3),",",round(Ub0,3),")",sep=""))
text(1985,9.6,paste("IC 95% b1 : (",round(Lb1,6),",",round(Ub1,6),")",sep=""))
text(2000,10.15,expression(paste(H[0]," : ",beta[1]," = ",-0.009,sep="")))
text(2000,10.1,expression(paste(H[1]," : ",H[0]," e' falsa",sep="")))
text(2000,10.05,paste("t=(",round(b1,7),"-(-0.009))/",round(sqrt(varb1),7),")",sep=""))
text(2000,10.0,paste("P-valor=",round(1-pt((b1+0.009)/sqrt(varb1),n-2),3),sep=""))


dev.off()