# rm(list=ls())
options(OutDec = ",",digits=5)
#==============================================================================
# Exemplo:
# Funcao de producao translog
#==============================================================================
# Analise: funcao de producao
#==============================================================================
dataxy <- read.table("../dados/d07_funcao_producao.txt",header=T)
n <- dim(dataxy)[1]
lnprod <- log(dataxy[,1])
lntrab <- log(dataxy[,2])
lncapi <- log(dataxy[,3])
lntrabx <- 0.5*(lntrab^2)
lncapix <- 0.5*(lncapi^2)
SST <- t(lnprod-mean(lnprod))%*%(lnprod-mean(lnprod))
#==============================================================================
# Analise do modelo 1: irrestrito (maior)
#==============================================================================
yfit1 <- lm(lnprod ~ lntrab + lncapi + lntrabx + lncapix + lntrab:lncapi)
print("==== ==== ==== ==== ==== ====")
print(summary(yfit1))
X <- cbind(lntrab,lncapi,lntrabx,lncapix,lntrab*lncapi)
print("==== ==== ==== ==== ==== ====")
print(anova(yfit1))
yfit1x <- lm(lnprod ~ X)
print("==== ==== ==== ==== ==== ====")
print(anova(yfit1x)) # Tabela anova
print("==== ==== ==== ==== ==== ====")
print(vcov(yfit1)) # Matriz de covariancia
print("==== ==== ==== ==== ==== ====")
print(confint(yfit1)) # Intervalo de confianca
R2.yfit1 <- 1 - t(yfit1$res)%*%yfit1$res / SST
[1] "==== ==== ==== ==== ==== ===="
Call:
lm(formula = lnprod ~ lntrab + lncapi + lntrabx + lncapix + lntrab:lncapi)
Residuals:
Min 1Q Median 3Q Max
-0,3399 -0,1011 -0,0124 0,0461 0,3928
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0,9442 2,9108 0,32 0,749
lntrab 3,6136 1,5481 2,33 0,030 *
lncapi -1,8931 1,0163 -1,86 0,077 .
lntrabx -0,9641 0,7074 -1,36 0,187
lncapix 0,0853 0,2926 0,29 0,774
lntrab:lncapi 0,3124 0,4389 0,71 0,484
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Residual standard error: 0,18 on 21 degrees of freedom
Multiple R-squared: 0,955, Adjusted R-squared: 0,944
F-statistic: 88,8 on 5 and 21 DF, p-value: 2,12e-13
[1] "==== ==== ==== ==== ==== ===="
Analysis of Variance Table
Response: lnprod
Df Sum Sq Mean Sq F value Pr(>F)
lntrab 1 13,52 13,52 417,69 2,4e-15 ***
lncapi 1 0,69 0,69 21,24 0,00015 ***
lntrabx 1 0,00 0,00 0,04 0,83928
lncapix 1 0,15 0,15 4,75 0,04074 *
lntrab:lncapi 1 0,02 0,02 0,51 0,48448
Residuals 21 0,68 0,03
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
[1] "==== ==== ==== ==== ==== ===="
Analysis of Variance Table
Response: lnprod
Df Sum Sq Mean Sq F value Pr(>F)
X 5 14,38 2,877 88,8 2,1e-13 ***
Residuals 21 0,68 0,032
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
[1] "==== ==== ==== ==== ==== ===="
(Intercept) lntrab lncapi lntrabx lncapix lntrab:lncapi
(Intercept) 8,47249 -2,387903 -0,331293 -0,08760 -0,233173 0,36354
lntrab -2,38790 2,396529 -1,231016 -0,66580 0,034767 0,18311
lncapi -0,33129 -1,231016 1,032787 0,52305 0,026369 -0,22554
lntrabx -0,08760 -0,665804 0,523052 0,50039 0,146743 -0,28803
lncapix -0,23317 0,034767 0,026369 0,14674 0,085620 -0,11604
lntrab:lncapi 0,36354 0,183113 -0,225542 -0,28803 -0,116040 0,19266
[1] "==== ==== ==== ==== ==== ===="
2,5 % 97,5 %
(Intercept) -5,10905 6,99744
lntrab 0,39425 6,83303
lncapi -4,00654 0,22032
lntrabx -2,43514 0,50704
lncapix -0,52322 0,69381
lntrab:lncapi -0,60041 1,22519
#==============================================================================
# Analise do modelo 2: restrito (menor)
#==============================================================================
yfit2 <- lm(lnprod ~ lntrab + lncapi)
print("==== ==== ==== ==== ==== ====")
print(summary(yfit2))
X <- cbind(lntrab,lncapi)
print("==== ==== ==== ==== ==== ====")
print(anova(yfit2))
yfit2x <- lm(lnprod ~ X)
print("==== ==== ==== ==== ==== ====")
print(anova(yfit2x)) # Tabela anova
print("==== ==== ==== ==== ==== ====")
print(vcov(yfit2)) # Matriz de covariancia
print("==== ==== ==== ==== ==== ====")
print(confint(yfit2)) # Intervalo de confianca
R2.yfit2 <- 1 - t(yfit2$res)%*%yfit2$res / SST
[1] "==== ==== ==== ==== ==== ===="
Call:
lm(formula = lnprod ~ lntrab + lncapi)
Residuals:
Min 1Q Median 3Q Max
-0,3038 -0,1012 -0,0182 0,0558 0,5056
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1,1706 0,3268 3,58 0,00150 **
lntrab 0,6030 0,1260 4,79 7,1e-05 ***
lncapi 0,3757 0,0853 4,40 0,00019 ***
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Residual standard error: 0,188 on 24 degrees of freedom
Multiple R-squared: 0,943, Adjusted R-squared: 0,939
F-statistic: 200 on 2 and 24 DF, p-value: 1,07e-15
[1] "==== ==== ==== ==== ==== ===="
Analysis of Variance Table
Response: lnprod
Df Sum Sq Mean Sq F value Pr(>F)
lntrab 1 13,52 13,52 381,1 3,1e-16 ***
lncapi 1 0,69 0,69 19,4 0,00019 ***
Residuals 24 0,85 0,04
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
[1] "==== ==== ==== ==== ==== ===="
Analysis of Variance Table
Response: lnprod
Df Sum Sq Mean Sq F value Pr(>F)
X 2 14,21 7,11 200 1,1e-15 ***
Residuals 24 0,85 0,04
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
[1] "==== ==== ==== ==== ==== ===="
(Intercept) lntrab lncapi
(Intercept) 0,1067865 -0,0198354 0,0011889
lntrab -0,0198354 0,0158644 -0,0096162
lncapi 0,0011889 -0,0096162 0,0072839
[1] "==== ==== ==== ==== ==== ===="
2,5 % 97,5 %
(Intercept) 0,49620 1,84509
lntrab 0,34304 0,86296
lncapi 0,19956 0,55186
#==============================================================================
# Calculando a estatistica F com R^2 (irrestrito) e R*^2 (restrito)
# Testando H_0: beta_4 = beta_5 = beta_6 = 0
#==============================================================================
valorF <- (R2.yfit1-R2.yfit2)/3 / ( (1- R2.yfit1)/yfit1$df.residual )
print(c(valorF > qf(0.95,3,yfit2$df.residual)))
# Hipotese nula nao eh rejeitada.
# Conclusao: a funcao Cobb-Douglas eh mais apropriada.
[1] FALSE
#==============================================================================
# Verificando a diminuicao no R^2 (irrestrito) para R*^2 (restrito)
#==============================================================================
print(c(R2.yfit1)) # R^2 (irrestrito)
print(c(R2.yfit2)) # R*^2 (restrito)
[1] 0,95486 [1] 0,94346
#==============================================================================
# Testando H_0: beta_2 + beta_3 = 1 utilizando o teste t e o F
#==============================================================================
m <- 0.6030 + 0.3757 - 1
ep.m <- sqrt(0.1260^2+0.0853^2-2*0.009616)
valort <- m / ep.m
print(c(abs(valort) > qt(0.975,yfit2$df.residual)))
# Hipotese nula nao eh rejeitada.
valorF <- valort^2
print(c(valorF > qf(0.95,1,yfit2$df.residual)))
# Hipotese nula nao eh rejeitada.
[1] FALSE [1] FALSE