# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Equacao de investimento:
# Modelo: y = beta0 + beta1 T + beta2 G + epsilon
#==============================================================================
library("Matrix")
dados <- matrix(scan("../dados/d02_investimento.txt"),15,5,byrow=T)
y <- dados[,1]
n <- length(y)
X <- cbind(rep(1,n),dados[,2:3])
rankMatrix(X)
XtX <- t(X)%*%X
Xty <- t(X)%*%y
rankMatrix(XtX)
XtXInv <- solve(XtX)
betahat1 <- XtXInv%*%Xty
#==============================================================================
# Estimativa de minimos quadrados
#==============================================================================
print(betahat1)
[,1] [1,] -0,50063897 [2,] -0,01719844 [3,] 0,65372331
#==============================================================================
# Correlacao do b_i's
#==============================================================================
print(cov2cor(XtXInv))
[,1] [,2] [,3] [1,] 1,0000000 0,9644496 -0,9977924 [2,] 0,9644496 1,0000000 -0,9786266 [3,] -0,9977924 -0,9786266 1,0000000
#==============================================================================
# Analisando a tendencia e PIB
#==============================================================================
plot(X[,2],X[,3],pch=15,xlab="Tendência",ylab="PIB",col="red")
# Dados ordenados no tempo. Possiveis problemas!
#==============================================================================
# Agora, no exemplo temos o modelo (mais as covariaveis R e P):
# y = beta0 + beta1 T + beta2 G + beta3 R + beta4 P + epsilon
#==============================================================================
X <- cbind(rep(1,n),dados[,2:5])
rankMatrix(X)
XtX <- t(X)%*%X
Xty <- t(X)%*%y
rankMatrix(XtX)
XtXInv <- solve(XtX)
betahat2 <- XtXInv%*%Xty
#==============================================================================
# Estimativa de minimos quadrados
#==============================================================================
print(betahat2)
[,1] [1,] -0,5090707909 [2,] -0,0165803945 [3,] 0,6703834376 [4,] -0,0023259283 [5,] -0,0000940107
#==============================================================================
# Correlacao do b_i's
#==============================================================================
print(cov2cor(XtXInv))
[,1] [,2] [,3] [,4] [,5] [1,] 1,00000000 0,9415070 -0,99286402 0,08332984 -0,04318183 [2,] 0,94150703 1,0000000 -0,93846805 -0,12002463 -0,01602380 [3,] -0,99286402 -0,9384680 1,00000000 -0,10858380 -0,03075184 [4,] 0,08332984 -0,1200246 -0,10858380 1,00000000 -0,45708148 [5,] -0,04318183 -0,0160238 -0,03075184 -0,45708148 1,00000000
#==============================================================================
# Comparando os coeficiente da primeira e segunda regressoes
#==============================================================================
print("==== Coeficientes da primeira regressão ====")
print(betahat1)
print("==== Coeficientes da segunda regressão ====")
print(betahat2)
[1] "==== Coeficientes da primeira regressão ===="
[,1]
[1,] -0,50063897
[2,] -0,01719844
[3,] 0,65372331
[1] "==== Coeficientes da segunda regressão ===="
[,1]
[1,] -0,5090707909
[2,] -0,0165803945
[3,] 0,6703834376
[4,] -0,0023259283
[5,] -0,0000940107
#==============================================================================
# Note que os coeficientes diferem de uma regressao para outra
# Isto se deve ao fato das regressoras nao serem ortogonais
#==============================================================================
print(cor(X[2:5,2:5]))
[,1] [,2] [,3] [,4] [1,] 1,0000000 0,9125219 -0,9262413 -0,8183521 [2,] 0,9125219 1,0000000 -0,9293711 -0,9818378 [3,] -0,9262413 -0,9293711 1,0000000 0,8755363 [4,] -0,8183521 -0,9818378 0,8755363 1,0000000
#============================================================================
# Calcularemos os coeficientes de correlacao parcial
# y = beta0 + beta1 T + beta2 G + beta3 R + beta4 P + epsilon
#============================================================================
# Correlacao simples entre investimento e as demais variaveis
#============================================================================
cor.inv <- cor(dados)[1,2:5]
print(cor.inv)
[1] 0,7495873 0,8632077 0,5871350 0,4777117
#============================================================================
# Correlacao parcial
#============================================================================
cor.par <- rep(NA,4)
#============================================================================
# Correlacao parcial entre investimento e G, R, e P
#============================================================================
X <- cbind(rep(1,n),dados[,3:5])
XtXInv <- solve(t(X)%*%X)
M <- diag(rep(1,n))- X%*%XtXInv%*%t(X)
res1 <- M%*%y
res2 <- M%*%dados[,2]
cor.par[1] <- cor(res1,res2)
#============================================================================
# Correlacao parcial entre investimento e T, R, e P
#============================================================================
X <- cbind(rep(1,n),dados[,c(2,4,5)])
XtXInv <- solve(t(X)%*%X)
M <- diag(rep(1,n))- X%*%XtXInv%*%t(X)
res1 <- M%*%y
res2 <- M%*%dados[,3]
cor.par[2] <- cor(res1,res2)
#============================================================================
# Correlacao parcial entre investimento e T, G, e P
#============================================================================
X <- cbind(rep(1,n),dados[,c(2,3,5)])
XtXInv <- solve(t(X)%*%X)
M <- diag(rep(1,n))- X%*%XtXInv%*%t(X)
res1 <- M%*%y
res2 <- M%*%dados[,4]
cor.par[3] <- cor(res1,res2)
#============================================================================
# Correlacao parcial entre investimento e T, G, e R
#============================================================================
X <- cbind(rep(1,n),dados[,2:4])
XtXInv <- solve(t(X)%*%X)
M <- diag(rep(1,n))- X%*%XtXInv%*%t(X)
res1 <- M%*%y
res2 <- M%*%dados[,5]
cor.par[4] <- cor(res1,res2)
#============================================================================
# Combinando as correlações
#============================================================================
print(cbind(cor.inv,cor.par))
cor.inv cor.par [1,] 0,7495873 -0,93600197 [2,] 0,8632077 0,96795713 [3,] 0,5871350 -0,51666472 [4,] 0,4777117 -0,02205714
#============================================================================
# Calcularemos as somas de quadrados e o R^2
# y = beta0 + beta1 T + beta2 G + beta3 R + beta4 P + epsilon
# Analise para os dados completos
#============================================================================
n <- length(y)
X <- cbind(rep(1,n),dados[,2:5])
print("==== Verificando o posto da Matriz X ====")
print(rankMatrix(X)[1])
XtX <- t(X)%*%X
Xty <- t(X)%*%y
print("==== Verificando o posto da Matriz XtX ====")
print(rankMatrix(XtX)[1])
XtXInv <- solve(XtX)
M0 <- matrix(-1/n,n,n)
for(i in 1:n) M0[i,i] <- M0[i,i] + 1
XtM0X <- t(X)%*%M0%*%X
betahat <- XtXInv%*%Xty
print("==== Estimativas dos Coeficientes ====")
print(betahat)
print("==== Estrutura de Correlação das Estimativas ====")
print(cov2cor(XtXInv))
[1] "==== Verificando o posto da Matriz X ===="
[1] 5
[1] "==== Verificando o posto da Matriz XtX ===="
[1] 5
[1] "==== Estimativas dos Coeficientes ===="
[,1]
[1,] -0,5090707909
[2,] -0,0165803945
[3,] 0,6703834376
[4,] -0,0023259283
[5,] -0,0000940107
[1] "==== Estrutura de Correlação das Estimativas ===="
[,1] [,2] [,3] [,4] [,5]
[1,] 1,00000000 0,9415070 -0,99286402 0,08332984 -0,04318183
[2,] 0,94150703 1,0000000 -0,93846805 -0,12002463 -0,01602380
[3,] -0,99286402 -0,9384680 1,00000000 -0,10858380 -0,03075184
[4,] 0,08332984 -0,1200246 -0,10858380 1,00000000 -0,45708148
[5,] -0,04318183 -0,0160238 -0,03075184 -0,45708148 1,00000000
#============================================================================
# Somas de Quadrados
#============================================================================
SQReg <- t(betahat)%*%XtM0X%*%betahat
resid <- y - X%*%betahat
SQRes <- t(resid)%*%resid # sum(resid^2)
SQTot <- var(y)*(n-1)
print("==== Soma de Quadrados da Regressão ====")
print(c(SQReg)) # 4 graus de liberdade
print("==== Soma de Quadrados dos Resíduos ====")
print(c(SQRes)) # 10 graus de liberdade
print("==== Soma de Quadrados Totais ====")
print(c(SQTot)) # 14 graus de liberdade
[1] "==== Soma de Quadrados da Regressão ====" [1] 0,01590252 [1] "==== Soma de Quadrados dos Resíduos ====" [1] 0,0004508118 [1] "==== Soma de Quadrados Totais ====" [1] 0,01635333
#============================================================================
# Calcularemos as somas de quadrados e o R^2
# y = beta0 + beta1 T + beta2 G + beta3 R + beta4 P + epsilon
# Analise para os dados completos
# Utilizamos como alternativa as funcoes "lm" e "anova"
#============================================================================
# Observe a diferenca entre as duas formas abaixo:
options(decimals=10,digits=10)
#============================================================================
# A funcao 'lsfit' faz ajuste de minimos quadrados.
# Entretanto, recomenda-se o uso da funcao 'lm' por ser mais completa.
# A funcao 'lm' sera usada a partir daqui para ajuste de minimos quadrados.
#============================================================================
yfit1 <- lm(y ~ X[,2:5])
print("==== Resumo do Ajuste ====")
print(summary(yfit1))
print("==== ANOVA ====")
print(anova(yfit1))
[1] "==== Resumo do Ajuste ===="
Call:
lm(formula = y ~ X[, 2:5])
Residuals:
Min 1Q Median 3Q Max
-0,0102778515 -0,0022946023 0,0004118672 0,0029377204 0,0080417956
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0,5090707909 0,0551276899 -9,23439 3,2821e-06 ***
X[, 2:5]1 -0,0165803945 0,0019717611 -8,40893 7,5891e-06 ***
X[, 2:5]2 0,6703834376 0,0549972149 12,18941 2,5210e-07 ***
X[, 2:5]3 -0,0023259283 0,0012188677 -1,90827 0,085449 .
X[, 2:5]4 -0,0000940107 0,0013474804 -0,06977 0,945754
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Residual standard error: 0,006714252 on 10 degrees of freedom
Multiple R-squared: 0,972433, Adjusted R-squared: 0,9614062
F-statistic: 88,18825 on 4 and 10 DF, p-value: 9,332607e-08
[1] "==== ANOVA ===="
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
X[, 2:5] 4 0,0159025215 0,0039756304 88,18825 9,3326e-08 ***
Residuals 10 0,0004508118 0,0000450812
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
#============================================================================
# Outra forma
#============================================================================
yfit2 <- lm(y ~ X[,2]+X[,3]+X[,4]+X[,5])
print("==== Resumo do Ajuste ====")
print(summary(yfit2))
print("==== ANOVA ====")
print(anova(yfit2))
[1] "==== Resumo do Ajuste ===="
Call:
lm(formula = y ~ X[, 2] + X[, 3] + X[, 4] + X[, 5])
Residuals:
Min 1Q Median 3Q Max
-0,0102778515 -0,0022946023 0,0004118672 0,0029377204 0,0080417956
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0,5090707909 0,0551276899 -9,23439 3,2821e-06 ***
X[, 2] -0,0165803945 0,0019717611 -8,40893 7,5891e-06 ***
X[, 3] 0,6703834376 0,0549972149 12,18941 2,5210e-07 ***
X[, 4] -0,0023259283 0,0012188677 -1,90827 0,085449 .
X[, 5] -0,0000940107 0,0013474804 -0,06977 0,945754
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Residual standard error: 0,006714252 on 10 degrees of freedom
Multiple R-squared: 0,972433, Adjusted R-squared: 0,9614062
F-statistic: 88,18825 on 4 and 10 DF, p-value: 9,332607e-08
[1] "==== ANOVA ===="
Analysis of Variance Table
Response: y
Df Sum Sq Mean Sq F value Pr(>F)
X[, 2] 1 0,0091886286 0,0091886286 203,82405 5,6164e-08 ***
X[, 3] 1 0,0064991615 0,0064991615 144,16574 2,9056e-07 ***
X[, 4] 1 0,0002145120 0,0002145120 4,75835 0,054121 .
X[, 5] 1 0,0000002194 0,0000002194 0,00487 0,945754
Residuals 10 0,0004508118 0,0000450812
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1