open "greene11_3.gdt"
setobs 1 1 --cross-section
# Set the actual values of the parameters
scalar A = 180
scalar B = .25   
scalar G = 1

# Starting values
ols C 0 Y
genr alpha0 = $coeff(0)
genr beta0 = $coeff(Y)
genr gamma0 = 1


# Set the number of Simulated Samples
scalar NMC = 1000

# Create an empty matrix to store results
matrix coeffs = zeros(NMC, 3)
matrix vcvs   = zeros(NMC, 6)

# Create systematic portion of the model and log(y)
series C0 = A + B*Y^G
series lY = log(Y)

# The loop
set warnings off
loop i=1..NMC --quiet
    # generate new sample
    genr C = C0 + normal(0,10)
    # Initialize parameters
    alpha = alpha0
    beta  = beta0
    gamma = gamma0

    # Estimate parameters via NLS
    nls C = alpha + beta * Y^gamma 
        deriv alpha = 1
        deriv beta = Y^gamma
        deriv gamma = beta * Y^gamma * lY
    end nls --quiet

    # Collect the coefficients into a vector and matrix
    matrix coeffs[i,] = {alpha, beta, gamma}
    matrix vcvs[i,] = vech($vcv)'
endloop

# open a new, empty dataset
nulldata NMC --preserve

# Convert the columns of matrix to data
series a = coeffs[,1]
series b = coeffs[,2]
series c = coeffs[,3] 

# Print the summary Statistics
summary a b c

printf "Monte Carlo vcv vs average estimated vcv\n"

MCV = mcov(coeffs)
AEV = unvech(meanc(vcvs)')

print MCV AEV