In [1]:
# rm(list=ls())
options(OutDec = ",")
#===================================================================
# Construindo exemplo com dados artificiais.
#===================================================================
set.seed(1234)
n <- 1000
beta0 <- 2
beta1 <- 0.3
x <- rnorm(n)
lambda <- exp(beta0 + beta1*x)
y <- rpois(n,lambda)
In [2]:
#===================================================================
# Nas análises abaixo, você deve ignorar o conhecimento do processo
# gerador dos dados.
#===================================================================
# Teste de correlação
#===================================================================
print(cor.test(x,y))
Pearson's product-moment correlation
data: x and y
t = 26,009, df = 998, p-value < 2,2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0,5971316 0,6711445
sample estimates:
cor
0,6355962
In [3]:
#===================================================================
# Histograma dos dados y e dispersão de x contra y.
#===================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(4.5,4.5,1,1),cex.main=2.0,bty="n")
hist(y,main="",prob=T,lwd=2,col="darkorange",xlab=expression(y),
ylab="Densidade",xlim=c(0,30))
plot(x,y,pch=15,lwd=2,col="red",xlim=c(-4,4),ylim=c(0,30),
xlab=expression(x),ylab=expression(y))
In [4]:
#===================================================================
# Ajuste do modelo de regressão linear:
# y = beta_0 + beta_1*x + erro.
#===================================================================
ajuste <- lm(y ~ x)
print(summary(ajuste))
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-9,1162 -1,9212 -0,2016 1,7298 14,3707
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7,80615 0,08815 88,56 <2e-16 ***
x 2,29912 0,08840 26,01 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0,001 ‘**’ 0,01 ‘*’ 0,05 ‘.’ 0,1 ‘ ’ 1
Residual standard error: 2,787 on 998 degrees of freedom
Multiple R-squared: 0,404, Adjusted R-squared: 0,4034
F-statistic: 676,4 on 1 and 998 DF, p-value: < 2,2e-16
In [5]:
#===================================================================
# Comparar os coeficientes com exp(beta0) e exp(beta1+sigma2/2).
#===================================================================
print(c(exp(beta0),exp(beta1+1/2)))
[1] 7,389056 2,225541
In [6]:
#===================================================================
# Análises dos resíduos.
#===================================================================
res <- ajuste$residual
yajustado <- ajuste$fitted.values
print(mean(res))
[1] -6,534351e-17
In [7]:
#===================================================================
# Histograma dos resíduos.
#===================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(4.5,4.5,1,1),cex.main=2.0,bty="n")
hist(res,main="",prob=T,xlab="Resíduos",ylab="Densidade",lwd=2,
col="darkorange",xlim=c(-15,20),ylim=c(0,0.16))
In [8]:
#===================================================================
# Boxplot dos resíduos.
#===================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(4.5,4.5,1,1),cex.main=2.0,bty="n")
boxplot(res,ylim=c(-15,20),col="darkorange")
In [9]:
#===================================================================
# qq-norm dos resíduos.
#===================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(4.5,4.5,1,1),cex.main=2.0,bty="n")
qqnorm(res,col="red",main="",xlab="Percentis teóricos",xlim=c(-4,4),
ylim=c(-15,20),ylab="Percentis amostrais",pch=16,xpd=F)
qqline(res,lwd=2,bty="n",lty=2)
In [10]:
#===================================================================
# Resíduos ao longo dos índices (tempo).
#===================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(6,8,5),
mar=c(4.5,4.5,1,1),cex.main=2.0,bty="n")
ts.plot(res,lwd=2,col="blue",ylim=c(-15,20),xlab=expression(i),
ylab=expression(e[i]))
points(res,pch=16,lwd=3,col=1)
In [11]:
#===================================================================
# Dispersão entre valores ajustados e os resíduos.
#===================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(4.5,4.5,1,1),cex.main=2.0,bty="n")
plot(yajustado,res,lwd=2,col="blue",xlim=c(0,16),ylim=c(-15,20),
pch=16,xlab=expression(hat(y)[i]),ylab=expression(e[i]))
print(cor(yajustado,res))
[1] -1,440739e-17
In [12]:
#===================================================================
# Dispersão entre valores da regressora e os resíduos.
#===================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(4.5,4.5,1,1),cex.main=2.0,bty="n")
plot(x,res,lwd=2,col="blue",ylim=c(-15,20),
pch=16,xlab=expression(x[i]),ylab=expression(e[i]))
print(cor(x,res))
[1] -2,336724e-17
In [13]:
#===================================================================
# Qual o problema encontrado nas analises graficas?
#===================================================================
In [14]:
#===================================================================
# Verificando se os resíduos são correlacionados utilizando a função
# de autocorrelação amostral
#===================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(4.5,4.5,1,1),cex.main=2.0,bty="n")
acf(res,lwd=2,col="black",main="",xlab="Defasagem",ylab="ACF",
xpd=F)
In [15]:
#===================================================================
# Verificando se os resíduos são correlacionados utilizando a função
# de autocorrelação amostral (quadrado dos resíduos).
# Homocedástico ou heterocedástico (ao longo das observações)?
#===================================================================
par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,8,5),
mar=c(4.5,4.5,1,1),cex.main=2.0,bty="n")
acf(res^2,lwd=2,col="black",main="",xlab="Defasagem",ylab="ACF",
xpd=F)
In [16]:
#===================================================================
# Teste formal de normalidade
#===================================================================
# Teste parametrico
silence <- suppressPackageStartupMessages
silence(library(tseries))
jarque.bera.test(res)
Jarque Bera Test data: res X-squared = 68,45, df = 2, p-value = 1,332e-15
In [17]:
#===================================================================
# Fim
#===================================================================