# rm(list=ls()) 
options(OutDec = ",")
#==============================================================================
# Exemplo:
#     Companhias aereas. 
#==============================================================================
dataxy     <- read.table("../dados/d09_companhias_aereas.txt",header=T)
print(head(dataxy))
n          <- dim(dataxy)[1]
k          <- dim(dataxy)[2]
lncit      <- log(dataxy[,3])
lnqit      <- log(dataxy[,4])
lnqitsq    <- lnqit^2
lnpit      <- log(dataxy[,5])
lfit       <- dataxy[,6]
SST        <- t(lncit-mean(lncit))%*%(lncit-mean(lncit))

  I T       C        Q     PF       LF
1 1 1 1140640 0,952757 106650 0,534487
2 1 2 1215690 0,986757 110307 0,532328
3 1 3 1309570 1,091980 110574 0,547736
4 1 4 1511530 1,175780 121974 0,540846
5 1 5 1676730 1,160170 196606 0,591167
6 1 6 1823740 1,173760 265609 0,575417

#==============================================================================
# Criando variaveis binarias
#==============================================================================
DB         <- matrix(0,n,14)
i          <- 1
j          <- 1
while(i <= n){
    DB[i,j] <- 1
    i       <- i + 1
    j       <- j + 1
    if( j==15){
        i <- i+1
        j <- 1 
    }
}
FB         <- matrix(0,n,5)
i          <- 1
j          <- 1
while(i <= (n-15)){
    FB[i,j] <- 1
    if( !(i %% 15) ){
        j <- j + 1 
    }
    i <- i + 1
}

#==============================================================================
# Analise do modelo
#==============================================================================
# Modelo sem efeitos # 14+5 = 19 restricoes
X          <- cbind(lnqit,lnqitsq,lnpit,lfit)
yfit1      <- lm(lncit ~ X)
SSE1       <- sum(yfit1$res^2)
print(SSE1)

# Modelo com efeitos da firma somente # 14 restricoes
X          <- cbind(lnqit,lnqitsq,lnpit,lfit,FB)
yfit2      <- lm(lncit ~ X)
SSE2       <- sum(yfit2$res^2)
print(SSE2)

# Modelo com efeitos temporais somente # 5 restricoes
X          <- cbind(lnqit,lnqitsq,lnpit,lfit,DB)
yfit3      <- lm(lncit ~ X)
SSE3       <- sum(yfit3$res^2)
print(SSE3)

# Modelo completo # 0 restricao
X          <- cbind(lnqit,lnqitsq,lnpit,lfit,DB,FB)
yfit4      <- lm(lncit ~ X)
SSE4       <- sum(yfit4$res^2)
print(SSE4)

# Os coeficientes sao conjuntamente iguais a zero
F1         <- ((SSE1 - SSE4) / 19 )/ ( SSE4/yfit4$df  )
print(F1 > qf(0.95,19,yfit4$df)) # O modelo completo eh melhor
# Efeitos da firma e temporais sao importantes

# Os coeficientes sao conjuntamente iguais a zero
F2         <- ((SSE2 - SSE4) / 14 )/ ( SSE4/yfit4$df  )
print(F2 > qf(0.95,14,yfit4$df)) # O modelo completo eh melhor
# Efeitos temporais sao importantes

# Os coeficientes sao conjuntamente iguais a zero
F3         <- ((SSE3 - SSE4) / 5 )/ ( SSE4/yfit4$df  )
print(F3 > qf(0.95,5,yfit4$df)) # O modelo completo eh melhor
# Efeitos da firma sao importantes

[1] 1,274921
[1] 0,2681542
[1] 1,034704
[1] 0,1725702
[1] TRUE
[1] TRUE
[1] TRUE