In [1]:
# rm(list=ls())
options(OutDec = ",")
#===================================================================
# Exemplo do investimento nos EUA de jan/1950 a dez/2000.
#===================================================================
# Parte I - Ajustes e testes.
#===================================================================
# Análise: Equação de investimento
#===================================================================
dataxy <- read.table("../dados/exemplo_10.txt",header=T)
print(head(dataxy))
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 6 1951 2 1803,7 1091,4 249,2 480,1 1225,8 77,6 116,19 1,53 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 6 3,1 152,917 1,5494 -0,0161
In [2]:
#===================================================================
# Análise do modelo 1: irrestrito.
#===================================================================
yfit1 <- lm(lninvest ~ txjuros + inflacao + lnPIB + tempo)
R2.yfit1 <- 1 - t(yfit1$res)%*%yfit1$res / SST
print(summary(yfit1))
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
In [3]:
#===================================================================
# Tabela anova
#===================================================================
X <- cbind(txjuros,inflacao,lnPIB,tempo)
print("===========================================================")
print(anova(yfit1)) # versão 1
yfit1x <- lm(lninvest ~ X)
print("===========================================================")
print(anova(yfit1x)) # versão 2
[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
In [4]:
#===================================================================
# Matriz de covariâncias.
#===================================================================
options(digits=4)
vcovfit1 <- vcov(yfit1)
print(vcovfit1)
(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
In [5]:
#===================================================================
# Intervalos de confiança.
#===================================================================
print(confint(yfit1))
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
In [6]:
#===================================================================
# Analise do modelo 2: restrito
#===================================================================
txjurosx <- txjuros - inflacao
yfit2 <- lm(lninvest ~ txjurosx + lnPIB + tempo)
R2.yfit2 <- 1 - t(yfit2$res)%*%yfit2$res / SST
print(summary(yfit2))
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
In [7]:
#===================================================================
# Tabela anova
#===================================================================
X <- cbind(txjurosx,lnPIB,tempo)
print("===========================================================")
print(anova(yfit2)) # versão 1
yfit2x <- lm(lninvest ~ X)
print("===========================================================")
print(anova(yfit2x)) # versão 2
[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
In [8]:
#===================================================================
# Matriz de covariâncias.
#===================================================================
options(digits=4)
print(vcov(yfit2))
(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
In [9]:
#===================================================================
# Intervalos de confiança.
#===================================================================
print(confint(yfit2))
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
In [10]:
#===================================================================
# Testando a restrição através do teste t.
#===================================================================
qchapeu <- sum(yfit1$coef[2:3])
# -0.00860 + 0.00331
EP.q <- sqrt(vcovfit1[2,2]+vcovfit1[3,3]+2*vcovfit1[2,3])
# 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
In [11]:
#===================================================================
# 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(pf(c(valorF),3,(n-1)-5,lower.tail=F)) # p-valor
[1] TRUE [1] 6,593e-42
In [12]:
#===================================================================
# Análise do modelo 3: restrito.
# Calculando a estatística 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 )
print(summary(yfit3))
R2.yfit3 <- 1 - t(yfit3$res)%*%yfit3$res / SST
valorFx <- (R2.yfit1-R2.yfit3)/3 / ( (1- R2.yfit1)/yfit1$df.residual )
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
In [13]:
#===================================================================
# Verificando os valores da estatística F.
#===================================================================
print("===========================================================")
print(as.numeric(valorF))
print("===========================================================")
print(as.numeric(valorFx))
[1] "===========================================================" [1] 109,8 [1] "===========================================================" [1] 109,8
In [14]:
#===================================================================
# Verificando a diminuição no R^2 (irrestrito) para R*^2 (restrito).
#===================================================================
print("===========================================================")
print(R2.yfit1) # R^2 (irrestrito)
print("===========================================================")
print(R2.yfit3) # R*^2 (restrito)
[1] "==========================================================="
[,1]
[1,] 0,9798
[1] "==========================================================="
[,1]
[1,] 0,9461
In [15]:
#===================================================================
# Parte 2 - Predição.
#===================================================================
# c(cst,txjuros,inflacao,lnPIB,tempo)
x0 <- c(1,4.48,log(528.0/521.1)*4*100,log(9316.8),204)
# print(x0)
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)
# Estimativa pontual e erro padrão.
print(c(ypred, ep.ypred))
[1] 7,3312 0,0877
In [16]:
#===================================================================
# Intervalo de predição.
#===================================================================
print(c(ypred-talpha*ep.ypred, ypred+talpha*ep.ypred))
[1] 7,159 7,503
In [17]:
#===================================================================
# Fim
#===================================================================