/* TESTING FOR CRS IN A NONLINEAR C0BB-DOUGLAS PRODUCTION FN */

/* 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);

/* 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);

/* *********************************** */
/* Unrestricted Gauss-Newton Algorithm */
/* *********************************** */

loops = 1;
g = zeros(nobs,rows(b));
print "UNRESTRICTED 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-6;
	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;

/* ********************************* */
/* Restricted Gauss-Newton Algorithm */
/* ********************************* */

rb = {1, 0.5};
loops = 1;
rg = zeros(nobs,rows(b)-1);
print;
print "RESTRICTED GAUSS-NEWTON ALGORITHM";
print "       loops            b1              b2 ";
print loops rb[1,1] rb[2,1];
cc = 1;  
 
do while cc>1e-6;
	h = (L^rb[2,1]).*(K^(1-rb[2,1]));
	rg[.,1] = h;
	rg[.,2] = (ln(L)-ln(K)).*h*rb[1,1];
	rresid = Q - rb[1,1]*h; 
	oldrb = rb;
	rb = rb + inv(rg'*rg)*(rg'*rresid);
	cc = (rb - oldrb)'*(rb - oldrb);	
	loops = loops + 1;
	print loops rb[1,1] rb[2,1];
endo;

/* ****************** */
/* Hypothesis Testing */ 
/* ****************** */

/* Asymptotically Valid F Test*/

sb = resid'*resid;
srb = rresid'*rresid;
F = (srb - sb)/(sb/(nobs-3));
print;
print "Asymptotically Valid F Statistic = " F;
print;

/* Wald Test */

R = {0 1 1};
s2 = sb/nobs;
V = s2*inv(g'*g);
W = (R*b - 1)'*inv(R*V*R')*(R*b - 1);
print "Wald Statistic = " W;
print;

/* Likelihood Ratio Test */

LR = nobs*(ln(srb)-ln(sb));
print "Likelihood Ratio Statistic = " LR;
print;

/* Lagrange Multiplier Test */

h = (L^rb[2,1]).*(K^(1-rb[2,1]));
g[.,1] = h;
g[.,2] = ln(L).*h*rb[1,1];
g[.,3] = ln(K).*h*rb[1,1];
LM = (rresid'*g*inv(g'*g)*g'*rresid)/(srb/nobs);
print "Lagrange Multiplier Test = " LM;