# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo:
# Investimento nos EUA de jan/1950 a dez/2000.
#==============================================================================
# Parte I - Ajustes e testes
#==============================================================================
# Analise: Equacao de Investimento
#==============================================================================
dataxy <- read.table("../dados/d06_investimento_eua.txt",header=T)
print(head(dataxy[1:5,]))
n <- dim(dataxy)[1]
k <- dim(dataxy)[2]
lninvest <- log(dataxy[2:n,5])
txjuros <- dataxy[2:n,10]
inflacao <- diff(log(dataxy[,8]))*4*100 # Dados trimestrais, 100%
lnPIB <- log(dataxy[2:n,3])
tempo <- 1:(n-1)
SST <- t(lninvest-mean(lninvest))%*%(lninvest-mean(lninvest))
Year qtr realgdp realcons realinvs realgovt realdpi cpi_u M1 tbilrate 1 1950 1 1610,5 1058,9 198,1 361,0 1186,1 70,6 110,20 1,12 2 1950 2 1658,8 1075,9 220,4 366,4 1178,1 71,4 111,75 1,17 3 1950 3 1723,0 1131,0 239,7 359,6 1196,5 73,2 112,95 1,23 4 1950 4 1753,9 1097,6 271,8 382,5 1210,0 74,9 113,93 1,35 5 1951 1 1773,5 1122,8 242,9 421,9 1207,9 77,3 115,08 1,40 unemp pop infl realint 1 6,4 149,461 0,0000 0,0000 2 5,6 150,260 4,5071 -3,3404 3 4,6 151,064 9,9590 -8,7290 4 4,2 151,871 9,1834 -7,8301 5 3,5 152,393 12,6160 -11,2160
#==============================================================================
# Analise do modelo 1: irrestrito
#==============================================================================
yfit1 <- lm(lninvest ~ txjuros + inflacao + lnPIB + tempo)
print(summary(yfit1))
print("==== ==== ==== ==== ==== ====")
X <- cbind(txjuros,inflacao,lnPIB,tempo)
print(anova(yfit1))
print("==== ==== ==== ==== ==== ====")
yfit1x <- lm(lninvest ~ X)
print(anova(yfit1x)) # Tabela anova
print("==== ==== ==== ==== ==== ====")
options(digits=4)
print(vcov(yfit1)) # Matriz de covariancia
print("==== ==== ==== ==== ==== ====")
print(confint(yfit1)) # Intervalo de confianca
R2.yfit1 <- 1 - t(yfit1$res)%*%yfit1$res / SST
Call:
lm(formula = lninvest ~ txjuros + inflacao + lnPIB + tempo)
Residuals:
Min 1Q Median 3Q Max
-0,22313 -0,05540 -0,00312 0,04246 0,31989
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -9,134092 1,366459 -6,684 2,3e-10 ***
txjuros -0,008598 0,003196 -2,691 0,00774 **
inflacao 0,003306 0,002337 1,415 0,15872
lnPIB 1,930156 0,183272 10,532 < 2e-16 ***
tempo -0,005659 0,001488 -3,803 0,00019 ***
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Residual standard error: 0,08618 on 198 degrees of freedom
Multiple R-squared: 0,9798, Adjusted R-squared: 0,9793
F-statistic: 2395 on 4 and 198 DF, p-value: < 2,2e-16
[1] "==== ==== ==== ==== ==== ===="
Analysis of Variance Table
Response: lninvest
Df Sum Sq Mean Sq F value Pr(>F)
txjuros 1 19,631 19,631 2643,118 < 2,2e-16 ***
inflacao 1 1,575 1,575 212,079 < 2,2e-16 ***
lnPIB 1 49,845 49,845 6711,270 < 2,2e-16 ***
tempo 1 0,107 0,107 14,463 0,0001904 ***
Residuals 198 1,471 0,007
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
[1] "==== ==== ==== ==== ==== ===="
Analysis of Variance Table
Response: lninvest
Df Sum Sq Mean Sq F value Pr(>F)
X 4 71,158 17,7896 2395,2 < 2,2e-16 ***
Residuals 198 1,471 0,0074
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
[1] "==== ==== ==== ==== ==== ===="
(Intercept) txjuros inflacao lnPIB tempo
(Intercept) 1,8672103 1,079e-03 8,312e-04 -0,2504208 2,026e-03
txjuros 0,0010789 1,021e-05 -3,717e-06 -0,0001465 9,861e-07
inflacao 0,0008312 -3,717e-06 5,462e-06 -0,0001121 9,697e-07
lnPIB -0,2504208 -1,465e-04 -1,121e-04 0,0335887 -2,718e-04
tempo 0,0020256 9,861e-07 9,697e-07 -0,0002718 2,214e-06
[1] "==== ==== ==== ==== ==== ===="
2,5 % 97,5 %
(Intercept) -11,828773 -6,439411
txjuros -0,014899 -0,002296
inflacao -0,001302 0,007915
lnPIB 1,568740 2,291572
tempo -0,008593 -0,002724
#==============================================================================
# Analise do modelo 2: restrito
#==============================================================================
txjurosx <- txjuros - inflacao
yfit2 <- lm(lninvest ~ txjurosx + lnPIB + tempo)
print(summary(yfit2))
print("==== ==== ==== ==== ==== ====")
X <- cbind(txjurosx,lnPIB,tempo)
print(anova(yfit2))
print("==== ==== ==== ==== ==== ====")
yfit2x <- lm(lninvest ~ X)
print(anova(yfit2x)) # Tabela anova
print("==== ==== ==== ==== ==== ====")
options(digits=4)
print(vcov(yfit2)) # Matriz de covariancia
print("==== ==== ==== ==== ==== ====")
print(confint(yfit2)) # Intervalo de confianca
R2.yfit2 <- 1 - t(yfit2$res)%*%yfit2$res / SST
Call:
lm(formula = lninvest ~ txjurosx + lnPIB + tempo)
Residuals:
Min 1Q Median 3Q Max
-0,22790 -0,05454 -0,00243 0,03999 0,31393
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -7,90716 1,20063 -6,59 3,9e-10 ***
txjurosx -0,00443 0,00227 -1,95 0,0526 .
lnPIB 1,76406 0,16056 10,99 < 2e-16 ***
tempo -0,00440 0,00133 -3,31 0,0011 **
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Residual standard error: 0,0867 on 199 degrees of freedom
Multiple R-squared: 0,979, Adjusted R-squared: 0,979
F-statistic: 3,15e+03 on 3 and 199 DF, p-value: <2e-16
[1] "==== ==== ==== ==== ==== ===="
Analysis of Variance Table
Response: lninvest
Df Sum Sq Mean Sq F value Pr(>F)
txjurosx 1 6,3 6,3 837,7 <2e-16 ***
lnPIB 1 64,8 64,8 8614,8 <2e-16 ***
tempo 1 0,1 0,1 10,9 0,0011 **
Residuals 199 1,5 0,0
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
[1] "==== ==== ==== ==== ==== ===="
Analysis of Variance Table
Response: lninvest
Df Sum Sq Mean Sq F value Pr(>F)
X 3 71,1 23,71 3154 <2e-16 ***
Residuals 199 1,5 0,01
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
[1] "==== ==== ==== ==== ==== ===="
(Intercept) txjurosx lnPIB tempo
(Intercept) 1,4415149 -4,319e-04 -1,928e-01 1,591e-03
txjurosx -0,0004319 5,154e-06 5,802e-05 -5,623e-07
lnPIB -0,1927647 5,802e-05 2,578e-02 -2,129e-04
tempo 0,0015911 -5,623e-07 -2,129e-04 1,771e-06
[1] "==== ==== ==== ==== ==== ===="
2,5 % 97,5 %
(Intercept) -10,274751 -5,540e+00
txjurosx -0,008903 5,015e-05
lnPIB 1,447442 2,081e+00
tempo -0,007027 -1,778e-03
#==============================================================================
# Testando a restricao atraves do teste t
#==============================================================================
qchapeu <- -0.00860 + 0.00331
EP.q <- sqrt( 0.00320^2 + 0.00234^2 + 2*(-3.717e-06) )
valort <- qchapeu / EP.q
print(abs(valort) > qt(0.975,(n-1)-5))
[1] FALSE
#==============================================================================
# Testando H_0: beta_1+beta_2=0, beta_4=1, e beta_5=0
#==============================================================================
R <- matrix(0,3,5)
R[1,2:3] <- c(1,1)
R[2,4] <- 1
R[3,5] <- 1
q <- c(0,1,0)
m <- R%*%yfit1$coef-q
estvarm <- R%*%vcov(yfit1)%*%t(R) # variancia estimada de m
valorF <- t(m)%*%solve(estvarm)%*%m/3
print(c(valorF > qf(0.95,3,(n-1)-5)))
print(c(pf(valorF,3,(n-1)-5,lower.tail=F))) # p-valor
[1] TRUE [1] 6,593e-42
#==============================================================================
# Analise do modelo 3: restrito
# Calculando a estatistica F com R^2 (irrestrito) e R*^2 (restrito)
# Testando H_0: beta_1+beta_2=0, beta_4=1, e beta_5=0
#==============================================================================
lninvestx <- lninvest - lnPIB # porque beta_4=1
yfit3 <- lm(lninvestx ~ txjurosx )
summary(yfit3)
print("==== ==== ==== ==== ==== ====")
R2.yfit3 <- 1 - t(yfit3$res)%*%yfit3$res / SST
valorFx <- (R2.yfit1-R2.yfit3)/3 / ( (1- R2.yfit1)/yfit1$df.residual )
# Verificando os valores da estatistica F
print(c(valorF))
print(c(valorFx))
print("==== ==== ==== ==== ==== ====")
# Verificando a diminuicao no R^2 (irrestrito) para R*^2 (restrito)
print(c(R2.yfit1)) # R^2 (irrestrito)
print(c(R2.yfit3)) # R*^2 (restrito)
Call:
lm(formula = lninvestx ~ txjurosx)
Residuals:
Min 1Q Median 3Q Max
-0,2915 -0,0977 -0,0070 0,0885 0,3627
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2,0165 0,0108 -187,32 <2e-16 ***
txjurosx 0,0069 0,0034 2,03 0,044 *
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Residual standard error: 0,14 on 201 degrees of freedom
Multiple R-squared: 0,0201, Adjusted R-squared: 0,0152
F-statistic: 4,12 on 1 and 201 DF, p-value: 0,0438
[1] "==== ==== ==== ==== ==== ====" [1] 109,8 [1] 109,8 [1] "==== ==== ==== ==== ==== ====" [1] 0,9798 [1] 0,9461
#==============================================================================
# Parte 2 - Predicao
#==============================================================================
# c(cst, txjuros, inflacao, lnPIB, tempo)
x0 <- c( 1, 4.48, log(528.0/521.1)*4*100, log(9316.8), 204 )
ypred <- t(x0)%*%yfit1$coef
s2 <- t(yfit1$res)%*%yfit1$res/yfit1$df.residual
ep.ypred <- sqrt( s2 + t(x0)%*%vcov(yfit1)%*%x0 )
#talpha <- qt(0.975,yfit1$df.residual)
talpha <- qnorm(0.975)
print(c(x0))
print(c(ypred, ep.ypred))
print(c(ypred-talpha*ep.ypred, ypred+talpha*ep.ypred)) # Intervalo de confianca
[1] 1,000 4,480 5,262 9,140 204,000 [1] 7,3312 0,0877 [1] 7,159 7,503