new;
cls;
@ this implementation is due to Mr Uldis, I had a minor role in getting the @
@ intercept right in one of the equations @
/* Input the data */

let XA[21,7] = 2.7	3.9		7.7		-10	12.7		182.8	44.9
		      2.9	3.2		3.9		-9	12.4		182.6	45.6
		      2.9	2.8		4.7		-8	16.9		184.5	50.1
		      3.1	3.5		3.8		-7	18.4		189.7	57.2
		      3.2	3.3		5.5		-6	19.4		192.7	57.1
		      3.3	3.3		7		-5	20.1		197.8	61
		      3.6	4		6.7		-4	19.6		203.4	64
		      3.7	4.2		4.2		-3	19.8		207.6	64.4
		      4	4.1		4		-2	21.1		210.6	64.5
		      4.2	5.2		7.7		-1	21.7		215.7	67
		      4.8	5.9		7.5		0	15.6		216.7	61.2
		      5.3	4.9		8.3		1	11.4		213.3	53.4
		      5.6	3.7		5.4		2	7		207.1	44.3
		      6	4		6.8		3	11.2		202		45.1
		      6.1	4.4		7.2		4	12.3		199		49.7
		      7.4	2.9		8.3		5	14		197.7	54.4
		      6.7	4.3		6.7		6	17.6		199.8	62.7
		      7.7	5.3		7.4		7	17.3		201.8	65
		      7.8	6.6		8.9		8	15.3		199.9	60.9
		      8	7.4		9.6		9	19		201.2	69.5
		      8.5	13.8		11.6		10	21.1		204.5	75.7;

x = ones(21,1)~XA;

x1 = x[.,1]~x[.,6];
x2 = x[.,1]~x[.,6]~x[.,7];
x3 = x[.,1]~x[.,5]~x[.,8];

let Y[21,7] = 	41.9		-0.2		25.5		45.6		12.4		182.6	28.2
		    	45		1.9		29.3		50.1		16.9		184.5	32.2
			49.2		5.2		34.1		57.2		18.4		189.7	37
			50.6		3		33.9		57.1		19.4		192.7	37
			52.6		5.1		35.4		61		20.1		197.8	38.6
			55.1		5.6		37.4		64		19.6		203.4	40.7
			56.2		4.2		37.9		64.4		19.8		207.6	41.5
			57.3		3		39.2		64.5		21.1		210.6	42.9
			57.8		5.1		41.3		67		21.7		215.7	45.3
			55		1		37.9		61.2		15.6		216.7	42.1
			50.9		-3.4		34.5		53.4		11.4		213.3	39.3
			45.6		-6.2		29		44.3		7		207.1	34.3
			46.5		-5.1		28.5		45.1		11.2		202		34.1
			48.7		-3		30.6		49.7		12.3		199		36.6
			51.3		-1.3		33.2		54.4		14		197.7	39.3
			57.7		2.1		36.8		62.7		17.6		199.8	44.2
			58.7		2		41		65		17.3		201.8	47.7
			57.5		-1.9		38.2		60.9		15.3		199.9	45.9
			61.6		1.3		41.6		69.5		19		201.2	49.4
			65		3.3		45		75.7		21.1		204.5	53
			69.7		4.9		53.3		88.4		23.5		209.4	61.8;

/*	COMPUTATION OF 2SLS ESTIMATES  */

/*	Compute the reduced-form estimates by LS  */
@following had olsqr@
b1 = (y[.,1]/ x);
b2 = (y[.,2]/ x);
b3 = (y[.,3]/ x);

/*	Compute the predicted values  */

yh1 = x*b1;
yh2 = x*b2;
yh3 = x*b3;

yh4 = yh1+yh2+x[.,3];
yh5 = yh4-x[.,4]-yh3;
yh6 = x[.,7]+yh2;
yh7 = yh3+x[.,2];

zh1 = yh5~yh7~x1;
zh2 = yh5~x2;
zh3 = yh4~x3;

/*	Compute the 2SLS estimates  */
/*
b2s1 = olsqr(y[.,1], zh1);
b2s2 = olsqr(y[.,2], zh2);
b2s3 = olsqr(y[.,3], zh3);
b2sls = b2s1|b2s2|b2s3;
*/

b2s1 = (y[.,1]/ zh1);
b2s2 = (y[.,2]/ zh2);
b2s3 = (y[.,3]/ zh3);
b2sls = b2s1|b2s2|b2s3;

/*	Compute the standard errors for 2SLS estimates  */
z1 = y[.,5]~y[.,7]~x1;
z2 = y[.,5]~x2;
z3 = y[.,4]~x3;

sh11 = (y[.,1] - z1*b2s1)'(y[.,1] - z1*b2s1)/21;
vh2sls1 = sh11*inv(zh1'zh1);
s2sls1 = sqrt(diag(vh2sls1));

sh22 = (y[.,2] - z2*b2s2)'(y[.,2] - z2*b2s2)/21;
vh2sls2 = sh22*inv(zh2'zh2);
s2sls2 = sqrt(diag(vh2sls2));

sh33 = (y[.,3] - z3*b2s3)'(y[.,3] - z3*b2s3)/21;
vh2sls3 = sh33*inv(zh3'zh3);
s2sls3 = sqrt(diag(vh2sls3));

/*	COMPUTATION OF LIML ESTIMATES */

/*	Compute the W0 and W1 matrices */
y01 = y[.,1]~y[.,5]~y[.,7];
y1 = y[.,5]~y[.,7];
y02 = y[.,2]~y[.,5];
y2 = y[.,5];
y03 = y[.,3]~y[.,4];
y3 = y[.,4];
I = eye(21);

W01 = y01'*(I - x1*inv(x1'x1)*x1')*y01;
W11 = y01'*(I - x*inv(x'x)*x')*y01;
W02 = y02'*(I - x2*inv(x2'x2)*x2')*y02;
W12 = y02'*(I - x*inv(x'x)*x')*y02;
W03 = y03'*(I - x3*inv(x3'x3)*x3')*y03;
W13 = y03'*(I - x*inv(x'x)*x')*y03;

/*	Compute the Eigenvalues and find lam1 */
D1 = inv(W11)*W01;
lam1 = eig(D1);
lam11 = minc(lam1);

D2 = inv(W12)*W02;
lam2 = eig(D2);
lam12 = minc(lam2);

D3 = inv(W13)*W03;
lam3 = eig(D3);
lam13 = minc(lam3);

/*	Compute the LIML estimates */
gliml1 = inv(W01[2:3,2:3] - lam11*W11[2:3,2:3])*(W01[1,2:3]' - lam11*W11[1,2:3]');
beliml1 = inv(x1'x1)*x1'(y[.,1] - y1*gliml1);
bliml1 = gliml1|beliml1;

gliml2 = inv(W02[2,2] - lam12*W12[2,2])*(W02[1,2] - lam12*W12[1,2]);
beliml2 = inv(x2'x2)*x2'(y[.,2] - y2*gliml2);
bliml2 = gliml2|beliml2;

gliml3 = inv(W03[2,2] - lam13*W13[2,2])*(W03[1,2] - lam13*W13[1,2]);
beliml3 = inv(x3'x3)*x3'(y[.,3] - y3*gliml3);
bliml3 = gliml3|beliml3;

/*	Compute the predicted values	*/

yhliml1 = x*inv(x'x)*x'y1;
yhliml2 = x*inv(x'x)*x'y2;
yhliml3 = x*inv(x'x)*x'y3;

zhliml1 = yhliml1~x1;
zhliml2 = yhliml2~x2;
zhliml3 = yhliml3~x3;

/*	Compute the standard errors for LIML estimates  */
shliml11 = (y[.,1] - z1*bliml1)'(y[.,1] - z1*bliml1)/21;
vhliml1 = shliml11*inv(zhliml1'zhliml1);
sliml1 = sqrt(diag(vhliml1));

shliml22 = (y[.,2] - z2*bliml2)'(y[.,2] - z2*bliml2)/21;
vhliml2 = shliml22*inv(zhliml2'zhliml2);
sliml2 = sqrt(diag(vhliml2));

shliml33 = (y[.,3] - z3*bliml3)'(y[.,3] - z3*bliml3)/21;
vhliml3 = shliml33*inv(zhliml3'zhliml3);
sliml3 = sqrt(diag(vhliml3));

/* 	Report  the results */

print "2SLS and LIML Parameter Estimates";
print;
print "				2SLS					Std Error				LIML					Std Error";
print;
print "Equation 1";
print;
print "Constant		" b2s1[3]~s2sls1[3]~bliml1[3]~sliml1[3];
print "P			" b2s1[1]~s2sls1[1]~bliml1[1]~sliml1[1];
print "P(-1)			" b2s1[4]~s2sls1[4]~bliml1[4]~sliml1[4];
print "Wp+Wg		" b2s1[2]~s2sls1[2]~bliml1[2]~sliml1[2];
print;
print "Equation 2";
print;
print "Constant		" b2s2[2]~s2sls2[2]~bliml2[2]~sliml2[2];
print "P			" b2s2[1]~s2sls2[1]~bliml2[1]~sliml2[1];
print "P(-1)			" b2s2[3]~s2sls2[3]~bliml2[3]~sliml2[3];
print "K(-1)			" b2s2[4]~s2sls2[4]~bliml2[4]~sliml2[4];
print;
print "Equation 3";
print;
print "Constant		" b2s3[2]~s2sls3[2]~bliml3[2]~sliml3[2];
print "X			" b2s3[1]~s2sls3[1]~bliml3[1]~sliml3[1];
print "X(-1)			" b2s3[4]~s2sls3[4]~bliml3[4]~sliml3[4];
print "A			" b2s3[3]~s2sls3[3]~bliml3[3]~sliml3[3];
print;

end;