3a Lista de Econometria Computacional

Desempenho dos estimadores \(\hat{\beta}_0\) e \(\hat{\beta}_1\)

# Vetor de Probabilidade, Gradiente e Hessiana
Matriz<-function(b0,b1,X,y,P,W){
    
             beta<-c(b0,b1)
          
             P <- 1 / ( 1+exp(-X %*% beta) )
             
             Grad<-t(X) %*% (y-P)
  
             Hess<- -t(X) %*% W %*% X
  
 return(list(prob=P, grad=Grad, hess=Hess))
}

# Parametros fornecidos 
b0<-0.2 ; b1<-0.5 ; beta<-c(b0,b1) 

# Método Monte Carlo com Newton-Rhapson
N=1000 ; Est.mc_b0<-numeric(N) ; Est.mc_b1<-numeric(N)
    
f.est<-function(n) {   
       
       for (j in 1:N){ 
   
          # Matriz de dados
          X<-matrix( c(rep(1,n),rnorm(n,0,1)), ncol=2 )

          # Vetor de Probabilidades       
          P <- 1 / ( 1+exp(-X %*% beta) ) 

          # Vetor Resposta
          y<-rbinom(n,1,P) 

          # Matriz de pesos
          W<-matrix(0,nrow=n,ncol=n) ; diag(W)<-P
      
          # equacao de atualizacao 
          tol=1e-5; check<-1 ; theta<-beta
          Mat<-Matriz(theta[1],theta[2],X,y,P,W)
      
          while( check > tol){
                
                Mat<-Matriz(theta[1],theta[2],X,y,P,W)
    
                theta <- theta - solve( Mat$hess ) %*%  Mat$grad 
                
                check <- sqrt( t(Mat$grad) %*% Mat$grad )
                         
                theta<-theta
           }
       Est.mc_b0[j]<-theta[1] ; Est.mc_b1[j]<-theta[2]
    }
return(c(mean(Est.mc_b0),mean(Est.mc_b0)-b0,mean(Est.mc_b1),mean(Est.mc_b1)-b1))
}

Resultados

## Instalando e Carregando os pacotes (bibliotecas)
lbs<-c('DT','kableExtra')
req <- substitute(require(x, character.only = TRUE))
sapply(lbs, function(x) eval(req) || {install.packages(x,dependencies = TRUE); eval(req)})

##         DT kableExtra 
##       TRUE       TRUE

df<-data.frame(f.est(25),f.est(60),f.est(120))
colnames(df)<-c("Estimativa (n=25)","Estimativa (n=60)", "Estimativa (n=120)")
row.names(df)<-c("beta0 est.","Vies0 est.","beta1 est.","Vies1 est.")

kable(df) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

	Estimativa (n=25)	Estimativa (n=60)	Estimativa (n=120)
beta0 est.	0.2018072	0.2115743	0.2048631
Vies0 est.	0.0018072	0.0115743	0.0048631
beta1 est.	0.6097041	0.5349989	0.5034653
Vies1 est.	0.1097041	0.0349989	0.0034653