rm(list=ls())
data = read.table("temco.txt",header=TRUE)
n = nrow(data)
attach(data)

########################################################################
# PARTE 1
#
# Modelo irrestrito: 
#
# SALARIO = beta0+beta1*EDUC+beta2*ANOSEMP+beta3*EXPPREV
#
# Modelo restrito:
#
# H0: beta1+beta2=2700
#
# H0: SALARIO-2700*ANOSEMP = beta0+beta1*(EDUC-ANOSEMP)+beta3*EXPPREV
#
######################################################################### 

# Modelo irrestrito
# -----------------
y     = SALARIO
X     = cbind(1,EDUC,ANOSEMP,EXPPREV)
iXtX  = solve(t(X)%*%X)
bhat  = iXtX%*%t(X)%*%y
e     = y-X%*%bhat
SSR   = t(e)%*%e

# Modelo restrito
educmanosemp = EDUC-ANOSEMP
y1    = SALARIO-2700*ANOSEMP
X1    = cbind(1,educmanosemp,EXPPREV)
bhat1 = solve(t(X1)%*%X1)%*%t(X1)%*%y1
e1    = y1-X1%*%bhat1
SSR1  = t(e1)%*%e1

# F test
num   = (SSR1-SSR)/1
den   = SSR/(n-4)
Ftest = num/den

# Alternativamente
# ----------------
R      = matrix(c(0,1,1,0),1,4)
r      = 2700
g      = nrow(R)
Rbr    = R%*%bhat-r
num1   = t(Rbr)%*%solve(R%*%iXtX%*%t(R))%*%(Rbr)/g
den1   = SSR/(n-4)
Ftest1 = num1/den1

Ftest
Ftest1

########################################################################
# PARTE 2
#
# Modelo irrestrito: 
#
# SALARIO = beta0 + beta1*EDUC + beta2*D1 + beta3*D2 + beta4*D3 + 
#           beta5*EDUC*D1 + beta6*EDUC*D2 + beta7*EDUC*D3
# onde
#
# D1 dummy para DEPT=1
# D2 dummy para DEPT=2
# D3 dummy para DEPT=3
#
# Modelo restrito:
#
# H0: beta5=beta6=beta7=0
#
# H0: SALARIO = beta0 + beta1*EDUC + beta2*D1 + beta3*D2 + beta4*D3
#
######################################################################### 
D1 = rep(0,n)
D2 = rep(0,n)
D3 = rep(0,n)
D1[DEPT==1]=1
D2[DEPT==2]=1
D3[DEPT==3]=1
ED1 = EDUC*D1
ED2 = EDUC*D2
ED3 = EDUC*D3

# Modelo irrestrito
# -----------------
y     = SALARIO
X     = cbind(1,EDUC,D1,D2,D3,ED1,ED2,ED3)
k     = ncol(X)-1
iXtX  = solve(t(X)%*%X)
bhat  = iXtX%*%t(X)%*%y
e     = y-X%*%bhat
SSR   = t(e)%*%e

# Modelo restrito
# ---------------
R        = matrix(0,3,8)
R[1,6]   = 1
R[2,7]   = 1
R[3,8]   = 1
r        = matrix(0,3,1)
g        = nrow(R)
Rbr      = R%*%bhat-r
num      = t(Rbr)%*%solve(R%*%iXtX%*%t(R))%*%(Rbr)/g
den      = SSR/(n-4)
Ftest    = num/den
Fcritico = qf(0.95,g,n-k-1)
pvalue   = 1-pf(Ftest,g,n-k-1)
c(Ftest,Fcritico,pvalue)