Question

Como fazer uma simulação Monte Carlo para explorar computacionalmente o fato que o estimador de mínimos quadrados é consistente?

Answer 1

Em python, eu fiz assim (se você ainda não programa em Python, você pode dar uma olhada aqui):

import random
import matplotlib.pyplot as plt
from scipy import stats
import numpy as np

def geracaoDados(n):
    alpha=0.5
    beta=1.1
    x=np.random.normal(0.5,0.25,n)
    u=np.random.normal(0,0.01,n)
    y=alpha*np.ones(n)+beta*x+u
    return y,x

def estimacaoOLS(x,y):
    beta, alpha, r_value, p_value, std_err = stats.linregress(x,y)
    return alpha,beta

if __name__ == '__main__':
    numeroAmostras=100
    n=30 # numero de observacoes
    vetor=np.empty([numeroAmostras])
    for i in range(0,numeroAmostras):
        [y,x]=geracaoDados(n)
        [alpha,beta]=estimacaoOLS(x,y)
        vetor[i]=beta
    axes=plt.subplot(221)    
    plt.setp(axes, xticks=[1.08,1.09, 1.1,1.11,1.12],yticks=[0, 5, 10,15])   
    plt.plot([1.1,1.1],[0,20],'r-') 
    plt.axis([1.08, 1.12, 0, 15]) 
    matplotlib.pyplot.hist(vetor,bins=30)
    plt.title('MQ: Amostra de tamanho 30')

    numeroAmostras=100
    n=100 # numero de observacoes
    vetor=np.empty([numeroAmostras])
    for i in range(0,numeroAmostras):
        [y,x]=geracaoDados(n)
        [alpha,beta]=estimacaoOLS(x,y)
        vetor[i]=beta
    axes=plt.subplot(222)    
    plt.setp(axes, xticks=[1.08,1.09, 1.1,1.11,1.12],yticks=[0, 5, 10,15])   
    plt.plot([1.1,1.1],[0,20],'r-') 
    plt.axis([1.08, 1.12, 0, 15]) 
    matplotlib.pyplot.hist(vetor,bins=30)
    plt.title('MQ: Amostra de tamanho 100')

    numeroAmostras=100
    n=250 # numero de observacoes
    vetor=np.empty([numeroAmostras])
    for i in range(0,numeroAmostras):
        [y,x]=geracaoDados(n)
        [alpha,beta]=estimacaoOLS(x,y)
        vetor[i]=beta
    axes=plt.subplot(223)    
    plt.setp(axes, xticks=[1.08,1.09, 1.1,1.11,1.12],yticks=[0, 5, 10,15])   
    plt.plot([1.1,1.1],[0,20],'r-') 
    plt.axis([1.08, 1.12, 0, 15]) 
    matplotlib.pyplot.hist(vetor,bins=30)
    plt.title('MQ: Amostra de tamanho 250')

    numeroAmostras=100
    n=1000 # numero de observacoes
    vetor=np.empty([numeroAmostras])
    for i in range(0,numeroAmostras):
        [y,x]=geracaoDados(n)
        [alpha,beta]=estimacaoOLS(x,y)
        vetor[i]=beta
    axes=plt.subplot(224)    
    plt.setp(axes, xticks=[1.08,1.09, 1.1,1.11,1.12],yticks=[0, 5, 10,15])   
    plt.plot([1.1,1.1],[0,20],'r-') 
    plt.axis([1.08, 1.12, 0, 15]) 
    matplotlib.pyplot.hist(vetor,bins=30)
    plt.title('MQ: Amostra de tamanho 1000')

Esse código gerará os seguintes painéis.