In [1]:
# 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
3
3
In [2]:
#==============================================================================
# Estimativa de minimos quadrados
#==============================================================================
print(betahat1)

            [,1]
[1,] -0,50063897
[2,] -0,01719844
[3,]  0,65372331
In [3]:
#==============================================================================
# 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
In [4]:
#==============================================================================
# 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!
In [5]:
#==============================================================================
# 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
5
5
In [6]:
#==============================================================================
# Estimativa de minimos quadrados
#==============================================================================
print(betahat2)

              [,1]
[1,] -0,5090707909
[2,] -0,0165803945
[3,]  0,6703834376
[4,] -0,0023259283
[5,] -0,0000940107
In [7]:
#==============================================================================
# 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
In [8]:
#==============================================================================
# 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
In [9]:
#==============================================================================
# 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
In [10]:
#============================================================================
# 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
In [11]:
#============================================================================
# Correlacao parcial
#============================================================================
cor.par <- rep(NA,4)
In [12]:
#============================================================================
# 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)
In [13]:
#============================================================================
# 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)
In [14]:
#============================================================================
# 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)
In [15]:
#============================================================================
# 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)
In [16]:
#============================================================================
# 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