# rm(list=ls()) 
options(OutDec = ",")
#==============================================================================
# Exemplo:
#     Propensao marginal a consumir (PMC) de longo prazo.
#
#     Embora os dados sejam os mesmo do exemplo_13.r, a selecao das 
#     variaveis pertinentes aqui sao diferentes.
#==============================================================================
dataxy     <- read.table("../dados/d06_investimento_eua.txt",header=T)
print(head(dataxy))
n          <- dim(dataxy)[1]
k          <- dim(dataxy)[2]
lncons     <- log(dataxy[2:n,4])
lnconstm1  <- log(dataxy[1:(n-1),4])
lnyt       <- log(dataxy[2:n,7])
SST        <- t(lncons-mean(lncons))%*%(lncons-mean(lncons))

  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

#==============================================================================
# Ajuste do modelo e teste de hipoteses
#==============================================================================
options(digits=7,decimals=7)
yfit1      <- lm(lncons ~ lnyt + lnconstm1)
print("==== ==== ==== ==== ==== ====")
print(summary(yfit1))
covb       <- vcov(yfit1)          # Matriz de covariancia 
delta      <- yfit1$coef[2] / (1 - yfit1$coef[3]) # PMC longo prazo
gbeta      <- 1 / (1 - yfit1$coef[3])
ggamma     <- delta / (1 - yfit1$coef[3])
g          <- c(gbeta,ggamma)
ep.delta   <- sqrt( t(g)%*%covb[2:3,2:3]%*%g )
valorz     <- (delta-1) / ep.delta
print("==== ==== ==== ==== ==== ====")
print(cbind(delta, ep.delta, abs(valorz), qnorm(0.975)))
print("==== ==== ==== ==== ==== ====")
print(c(abs(valorz) > qnorm(0.975)))
# A hipotese nula nao eh rejeitada.
# Assim, a hipotese de delta = 1 nao eh rejeitada.

[1] "==== ==== ==== ==== ==== ===="

Call:
lm(formula = lncons ~ lnyt + lnconstm1)

Residuals:
      Min        1Q    Median        3Q       Max 
-0,035243 -0,004606  0,000496  0,005147  0,041754 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 0,003142   0,010553   0,298  0,76624    
lnyt        0,074958   0,028727   2,609  0,00976 ** 
lnconstm1   0,924625   0,028594  32,337  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1

Residual standard error: 0,008742 on 200 degrees of freedom
Multiple R-squared:  0,9997,	Adjusted R-squared:  0,9997 
F-statistic: 3,476e+05 on 2 and 200 DF,  p-value: < 2,2e-16

[1] "==== ==== ==== ==== ==== ===="
         delta                              
lnyt 0,9944606 0,01635798 0,3386337 1,959964
[1] "==== ==== ==== ==== ==== ===="
[1] FALSE