#####################################################################################
#
# Variaveis:
# Analfabetos - pessoas 15 anos e mais (2000)	
# Renda per capita (2000)	
# Razao entre a renda dos 10% mais ricos e 40% mais pobres (2000)	
#
#####################################################################################
rm(list=ls())
nomes    = c("Analfabetos - 15 anos ou mais","Renda per capita",
             "Razao renda - top 10% e bottom 40%")
nomesl = c("Analfabetos - 15 anos ou  mais (LOG)",
             "Renda per capita (LOG)",
             "Razao renda - top 10% e bottom 40% (LOG)")
data  = read.csv("dadosmunicipais.csv",header=TRUE)
n     = nrow(data)
data1 = data[,3:5]
data2 = log(data[,3:5])


par(mfrow=c(1,3))
hist(data1[,1],prob=TRUE,main=nomes[1],xlab="");box()
hist(data1[,2],prob=TRUE,main=nomes[2],xlab="");box()
hist(data1[,3],prob=TRUE,main=nomes[3],xlab="");box()

par(mfrow=c(1,3))
hist(data2[,1],prob=TRUE,main=nomesl[1],xlab="");box()
hist(data2[,2],prob=TRUE,main=nomesl[2],xlab="");box()
hist(data2[,3],prob=TRUE,main=nomesl[3],xlab="");box()


par(mfrow=c(1,2))
plot(data1[,1],data1[,3],xlab=nomes[1],ylab=nomes[3])
plot(data1[,2],data1[,3],xlab=nomes[2],ylab=nomes[3])

par(mfrow=c(1,2))
plot(data1[,1],data1[,3],xlab=nomes[1],ylab=nomes[3],ylim=c(0,40))
text(50,0,paste(round(100*mean(data1[,3]<=40),2),"%",sep=""),col=2)
plot(data1[,2],data1[,3],xlab=nomes[2],ylab=nomes[3],ylim=c(0,40))

par(mfrow=c(1,2))
plot(data2[,1],data2[,3],xlab=nomesl[1],ylab=nomesl[3])
abline(lm(data2[,3]~data2[,1]),col=2,lwd=3)
plot(data2[,2],data2[,3],xlab=nomesl[2],ylab=nomesl[3])
abline(lm(data2[,3]~data2[,2]),col=2,lwd=3)

# Regression 1
# ------------
y   = data2[,3]
x1  = data2[,1]
x2  = data2[,2]
reg = lm(y~x1+x2)
summary(reg)

# RESET test
yfit2 = reg$fit^2
yfit3 = reg$fit^3
reg1  = lm(y~x1+x2+yfit2+yfit3)
summary(reg1)

e      = reg1$res
SSR    = t(e)%*%e
bhat   = reg1$coef
X      = cbind(1,x1,x2,yfit2,yfit3)
p      = ncol(X)
iXtX   = solve(t(X)%*%X)


g      = 2
R      = matrix(c(0,0,0,1,0,0,0,0,0,1),g,p,byrow=TRUE)
r      = matrix(0,2,1)
Rbr    = R%*%bhat-r
num   = t(Rbr)%*%solve(R%*%iXtX%*%t(R))%*%(Rbr)/g
den   = SSR/(n-p)
Ftest = num/den
criticalvalue = qf(0.95,g,n-p)
pvalue = 1-pf(Ftest,g,n-p)
c(Ftest,criticalvalue,pvalue)