Capítulo 14 Estimador de Maxima Verossimilhanca

14.1 Distribuição não condicional

**Distribuição Normal

\(l = -0.5 n ln(2\pi) - 0.5nln(\sigma^2) - \frac{1}{2\sigma^2} \sum_i (y_i-\mu)^2\)

1º passo: definir a função de log-verossimilhança

normal = function(theta,y){
  mu = theta[1]
  sigma2= theta[2]
  n = nrow(y)
  logl = -.5*n*log(2*pi) - .5*n*log(sigma2) - (1/(2*sigma2))*sum((y-mu)^2)
  return(-logl)
}

2º passo: extrair uma amostra aleatória da distribuição normal

n = 1000
y = data.frame(rnorm(n,10,5))

3º passo: encontrar os valores dos parâmetros que maximizam a função

par = optim(c(0,1),normal,y=y,method="BFGS")

## Warning in log(sigma2): NaNs produzidos

## Warning in log(sigma2): NaNs produzidos

## Warning in log(sigma2): NaNs produzidos

## Warning in log(sigma2): NaNs produzidos

## Warning in log(sigma2): NaNs produzidos

## Warning in log(sigma2): NaNs produzidos

par$par

## [1] 10.03299 26.03589

14.2 Distribuição Condicional

**Modelo linear

\(f(y|x) = \frac{1}{(2\pi\sigma^2)^{1/2}} exp (-\frac{1}{2\sigma^2}(y-x'\beta)^2)\)

ols.lf<-function(theta,y,X){
  n<-nrow(X)
  k<-ncol(X)
  beta<-theta[1:k]
  sigma2<-theta[k+1]
  e<- y-X%*%beta
  logl<- -(n/2)*log(2*pi*sigma2) - ((t(e)%*%e)/(2*sigma2))
  return(-logl)
}

n = 1000
X<-cbind(1,runif(n))
theta.true<-c(2,3,1)
y<-X%*%theta.true[1:2] + rnorm(n,0,theta.true[3])

p<-optim(c(1,1,1),ols.lf,method="BFGS",hessian=T,y=y,X=X)

## Warning in log(2 * pi * sigma2): NaNs produzidos

## Warning in log(2 * pi * sigma2): NaNs produzidos

## Warning in log(2 * pi * sigma2): NaNs produzidos

## Warning in log(2 * pi * sigma2): NaNs produzidos

## Warning in log(2 * pi * sigma2): NaNs produzidos

## Warning in log(2 * pi * sigma2): NaNs produzidos

## Warning in log(2 * pi * sigma2): NaNs produzidos

p$par

## [1] 2.0027331 3.0270868 0.9925956

OI<-solve(p$hessian)
OI

##               [,1]          [,2]          [,3]
## [1,]  3.917773e-03 -5.897796e-03  1.446704e-11
## [2,] -5.897796e-03  1.189125e-02 -1.975390e-11
## [3,]  1.446704e-11 -1.975390e-11  1.970483e-03

se<-sqrt(diag(OI))
se

## [1] 0.06259211 0.10904699 0.04439012

t = p$par/se
pval = 2*(1-pt(abs(t),nrow(X)-ncol(X)))
results = cbind(p$par,se,t,pval)
colnames(results) = c("b","se","t","p")
rownames(results) = c("Const","X1","Sigma2")
print(results,digits=3)

##            b     se    t p
## Const  2.003 0.0626 32.0 0
## X1     3.027 0.1090 27.8 0
## Sigma2 0.993 0.0444 22.4 0

B = solve(t(X)%*%X)%*%t(X)%*%y
B

##          [,1]
## [1,] 2.002733
## [2,] 3.027087

e_hat = y - X%*%B
k = ncol(X)
sigma_hat = (1/(n-k)) * t(e_hat)%*%e_hat   # n-k correção de viés
sigma_hat

##           [,1]
## [1,] 0.9945841

e_hat = y - X%*%B
k = ncol(X)
sigma_hat = (1/(n)) * t(e_hat)%*%e_hat   # n-k correção de viés
sigma_hat

##          [,1]
## [1,] 0.992595

V = sigma_hat[1,1] * solve(t(X)%*%X)
V 

##              [,1]         [,2]
## [1,]  0.003917770 -0.005897792
## [2,] -0.005897792  0.011891237

ep = sqrt(diag(V))
ep

## [1] 0.06259209 0.10904695

t = B / ep
t

##          [,1]
## [1,] 31.99658
## [2,] 27.75948

2 * (1 - pt(abs(t), df = n-k))

##      [,1]
## [1,]    0
## [2,]    0