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.
