import numpy as np

np.random.seed(1234)

# DGP
n =100
x = np.random.uniform(size=(n,1))
eps = np.random.standard_normal(size=(n,1))
beta  = [2, 2, 2]
y = beta[0]+beta[1]*x**beta[2]+eps

# DGP ENDS


# ESTIMATION

# Starting values

beta0=np.array([1, 1, 1])

tol=1
iter=1

while tol>0.0001:
    y0=y-beta0[0]-beta0[1]*x**beta0[2]
    # First derivatives
    xd1 = np.ones((n,1))
    xd2 = x**beta[2]
    xd3 = (beta0[1]*x**beta0[2])*np.log(x)
    xd=np.matrix(np.hstack([xd1,xd2,xd3]))
    #GNR 
    b=(xd.T*xd).I*xd.T*y0
    tol=b.T*b
    beta0=np.asarray(b.T)+beta0
    
    beta0=beta0.reshape(3)
    iter=iter+1

bnls=beta0

# Var Cov Matrix of the NLS estimator

def vcov(bnls):

    k=len(bnls)

    e=np.matrix(y-bnls[0]-bnls[2]*x**bnls[2])
    sig2hat=(e.T*e)/(n-k)
    
    #First derivatives
    xd1 = np.ones((n,1))
    xd2 = x**bnls[2]
    xd3 = (bnls[1]*x**bnls[2])*np.log(x)
    xdm=np.matrix(np.hstack([xd1,xd2,xd3]))
    return (np.multiply(sig2hat,(xdm.T*xdm).I))

#Standard Errors
serr=np.sqrt(np.diag(vcov(bnls)))


# OUTPUT
print '''
==========================
GNR Estimation Results
model : y = b1+b2+x^b3+eps
==========================
'''
bnlst = ['beta1','beta2','beta3']
print 'Param.','Coef.','  SERR' 
print bnlst[0],"%.4f"% bnls[0],"(%.4f)"% serr[0]
print bnlst[1],"%.4f"% bnls[1],"(%.4f)"% serr[1]
print bnlst[2],"%.4f"% bnls[2],"(%.4f)"% serr[2]
print '''
===========================
'''

print 'Number of iterations:', iter