/* NLS ESTIMATION OF COBB-DOUGLAS PRODUCTION FUNCTION */

/* Load Data and Define Variables */

load path = c:\gauss35\classes\econ5340\data\;
load x[30,3]=table12.3;
L= x[.,1];
K = x[.,2];
Q = x[.,3];
nobs = rows(x);

/* *********************** */
/* Gauss-Newton Estimation */
/* *********************** */

/* Initial Values */

constant = ones(nobs,1);
lnL = ln(L);
lnK = ln(K);
lnXmat = constant~lnL~lnK;
lnQ = ln(Q);
b = inv(lnXmat'*lnXmat)*(lnXmat'*lnQ);

/* Algorithm */

loops = 1;
g = zeros(nobs,rows(b));  
print "GAUSS-NEWTON ALGORITHM";
print "       loops            b1              b2               b3 ";
print loops b[1,1] b[2,1] b[3,1];
cc = 1;  

do while cc>1e-8;
	h = (L^b[2,1]).*(K^b[3,1]);
	g[.,1] = h;
	g[.,2] = ln(L).*h*b[1,1];
	g[.,3] = ln(K).*h*b[1,1];
	resid = Q - b[1,1]*h; 
	oldb = b;
	b = b + inv(g'*g)*(g'*resid);
	cc = (b - oldb)'*(b -oldb);	
	loops = loops +1;
	print loops b[1,1] b[2,1] b[3,1];
endo;

/* ************************** */
/* Newton's Method Estimation */
/* ************************** */

/* Initial Values */

b = {1,0.5,0.5};

/* Algorithm */

loops = 1;
g = zeros(rows(b),1);
Hs = zeros(rows(b),rows(b));  
print;
print "NEWTON-RAPHSON ALGORITHM";
print "       loops            b1              b2               b3 ";
print loops b[1,1] b[2,1] b[3,1];
cc = 1;  

do while cc>1e-6;
	h = (L^b[2,1]).*(K^b[3,1]);
	resid = Q - b[1,1]*h; 
	g[1,1] = -2*sumc(resid.*h);
	g[2,1] = -2*sumc(resid.*(ln(L).*b[1,1].*h));
	g[3,1] = -2*sumc(resid.*(ln(K).*b[1,1].*h));
	Hs[1,1] = 2*sumc(h^2);
	Hs[1,2] = -2*sumc(resid.*ln(L).*h - b[1,1]*(h^2).*ln(L));
	Hs[1,3] = -2*sumc(resid.*ln(K).*h - b[1,1]*(h^2).*ln(K));
	Hs[2,1] = Hs[1,2];
	Hs[2,2] = -2*sumc(resid.*(ln(L)^2).*h*b[1,1] - (b[1,1]^2)*(h^2).*(ln(L)^2));
	Hs[2,3] = -2*sumc(resid.*ln(L).*ln(K).*h*b[1,1] - (b[1,1]^2)*(h^2).*ln(L).*ln(K));
	Hs[3,1] = Hs[1,3];
	Hs[3,2] = Hs[2,3];
	Hs[3,3] = -2*sumc(resid.*(ln(K)^2).*h*b[1,1] - (b[1,1]^2)*(h^2).*(ln(K)^2));
	
	oldb = b;
	b = b - inv(Hs)*g;
	cc = (b - oldb)'*(b -oldb);	
	loops = loops + 1;
	print loops b[1,1] b[2,1] b[3,1];
endo;

/* ************* */
/* Graph Surface */
/* ************* */

b2 = seqa(0.5,0.02,30);
b3 = seqa(0.5,0.02,30);
ssresid = zeros(30,30);

for ii(1,30,1);
for jj(1,30,1); 
	resid = Q - b[1,1]*(L^b2[ii,1]).*(K^b3[jj,1]);
	ssresid[ii,jj] = resid'*resid;
endfor;
endfor;

library pgraph;
graphset;
title("Nonlinear Minimization of the Sum of Squared Residual");
xlabel("Beta2");
ylabel("Beta3");
zlabel("Sum of Squared Residuals");
surface(b2',b3,ssresid);