# ********************************************************
# Função 'multlm': calcula os coeficientes e matriz de 
# Variância-Covariância do modelo de regressão multivariada
# Model: 
#	Y = B_0 + X * B + U , U ~ N( 0 , Sigma )	
# Argumentos: 
#	Y : Matriz de dados de dimensão n x p
#	X : Matriz de covariáveis de dimensão n x q
# Saida:
#	B.hat : Matriz de coeficientes de dimensão q x p
#	Sigma.hat : Variância-covariância de dimensão p x p 
# *********************************************************

multlm <- function(Y,X){
	# model: Y = B_0 + X * B + U, U ~ N(0,Sigma)	
        p     <- ncol(Y)
        q     <- ncol(X)
        n     <- nrow(Y)
        meanX    <- c(1:q)
        for(i in 1:q) meanX[i] <- mean(X[,i])
        mX    <- t(matrix(meanX,q,n)) 
        meanY    <- c(1:p)
        for(i in 1:p) meanY[i] <- mean(Y[,i])
        mY    <- t(matrix(meanY,p,n)) 
        Xc    <- X - mX        
        Yc    <- Y - mY
        VarB  <- solve(t(Xc)%*%Xc)
	B.hat     <- (t(Yc)%*%Xc)%*%VarB
	Sigma.hat <- (t(Yc)%*%Yc)-(t(Yc)%*%Xc)%*%VarB%*%(t(Xc)%*%Yc)
        B_0.hat   <- meanY - (t(Yc)%*%Xc)%*%VarB%*%meanX
	res <- Y - ( t(matrix(B_0.hat,p,n)) + X%*%t(B.hat) )
        #testes F para cada regressão
        Freg      <- 1:p
        for(j in 1:p)  Freg[j] <- (n-(q+1))*(B.hat[j,]%*%(t(Xc)%*%Xc)%*%(B.hat[j,]))/(q*Sigma.hat[j,j])
        #teste F global para todas as p regressões
        globalF <- (n-(q+1))*sum(B.hat%*%(t(Xc)%*%Xc)%*%t(B.hat))/(q*sum(Sigma.hat))
        #testes individuais para cada coeficiente de regressao
        tcoef   <- matrix(0,p,q)
        for(i in 1:p){for (j in 1:q){tcoef[i,j] <- B.hat[i,j]/sqrt(Sigma.hat[i,i]*VarB[j,j]/(n-(q+1)))  }}
	output <- list(coef=B.hat,Sigma=Sigma.hat,resid=res,Ftest=Freg,Fglobal=globalF,ttest=tcoef)
	output
} 


# ********************************************************
# Exemplo: Dados macroeconómicos - J. W. Longley (1967) 
# Variáveis: 
# 	GNP.deflator: GNP implicit price deflator (1954=100)
#	GNP: Gross National Product.
# 	Unemployed: number of unemployed.
# 	Armed.Forces: number of people in the armed forces.
# 	Population: ‘noninstitutionalized’ population >= 14 years of age.
# 	Year: the year (time).
# 	Employed: number of people employed. 
# Variáveis dependentes: Y = (GNP deflator, GNP), X = (Unemployed, Armed.Forces, Population, Year, Employed)
# ********************************************************

# Dados:
data(longley)
longley
pairs(longley, main = "longley data")

# Y, X:
longley.y <- data.matrix(longley[, 1:2])
longley.x <- data.matrix(longley[, 3:7])

# Ajuste de regressão multivariada
reg.longley <- multlm(longley.y,longley.x)

reg.longley$coef
reg.longley$Sigma
reg.longley$resid

# Gráfico dos residuos:
par(mfrow=c(2,1))
plot(reg.longley$resid[,"GNP"],ylab="",main="Residuos: GNP.deflator")
abline(h=0,lty=2,col=2,lwd=2)
plot(reg.longley$resid[,"GNP"],ylab="",main="Residuos: GNP")
abline(h=0,lty=2,col=2,lwd=2)


# ************************************************
# Análise multivariada no R:
# 	- Análise de componentes principais (princomp)
# 	- Análise Fatorial (factanal)
# ************************************************