randomaccess;
{Dimension arrays}
	Dim b1(2);
	Dim y(10);
	Dim xy1(10,5);
	Dim xx1(5,5);
	Dim xy2(10,4);
	Dim xx2(4,4);
	Dim Predicted(10);
	Dim Residuals(10);
	Dim table (10,10); {for Scrollsheet}

	Dim work1(10);     {scratch}
	Dim work2(10,2);
	Dim work3(2,2);

{Stage 1}
	MatrixFill (1, xy1, 1, 1, Ncases, 1);          {constant}
	MatrixCopy (data, 1, 2, Ncases, 1, xy1, 1, 2); {Income}
	MatrixCopy (data, 1, 4, Ncases, 1, xy1, 1, 3); {Lagprice}
	MatrixCopy (data, 1, 5, Ncases, 1, xy1, 1, 4); {Instrum}
	MatrixCopy (data, 1, 3, Ncases, 1, xy1, 1, 5); {Price}
	MatrixCrossProductOfDev (xy1, 0, xx1);	     {Crossproducts}
	MatrixSweep (xx1, 1, 4, 1);

{Predicted Price values for stage 2}
	for i:=1 to Ncases do begin
	 for j:=1 to 4 do
	  xy2(i,2):=xy2(i,2)+xx1(5,j)*xy1(i,j);
	end;
{Stage 2}
	MatrixFill (1, xy2, 1, 1, Ncases, 1);          {constant}
	MatrixCopy (data, 1, 2, Ncases, 1, xy2, 1, 3); {Income}
	MatrixCopy (data, 1, 1, Ncases, 1, xy2, 1, 4); {Demand}
	MatrixCOpy (data, 1, 1, Ncases, 1, y  , 1, 1); {copy also to array y}
	MatrixCrossProductOfDev (xy2, 0, xx2);	     {Crossproducts}
	MatrixSweep (xx2, 1, 3, 1);                    {Coefficients are in xx2(4,1..3)}
	MatrixDisplay (xx2, 'Sweep Matrix');

{Predicted and residual values;
 NOTE: replace predicted values from stage 1 with original values for PRICE}
	MatrixCopy (data, 1, 3, Ncases, 1, xy2, 1, 2);
	for i:=1 to Ncases do begin
	 for j:=1 to 3 do predicted(i):=predicted(i)+xx2(4,j)*xy2(i,j);
	 residuals(i):=y(i)-predicted(i);
	end;
{copy to data file}
	MatrixCopy (predicted, 1, 1, Ncases, 1, data, 1, 6);
	MatrixCopy (residuals, 1, 1, Ncases, 1, data, 1, 7);
	MatrixSumOfSquares (residuals, sse);                  {SS residual}
	mse:=sse/(Ncases-3);
{regression summary}
	line01$:='Coefficients';
	line01$:=line01$+'|MSE='+Str(mse, 10, 7);
      kname2$:='Coeff | Std.Err. | t(7) | p |-95% Cnf|+95% Cnf';
	t:=VStudent (.975, 7);
	for i:=1 to 3 do begin;
	  table(i,1):=xx2(4,i);                     		{regression coefficients }
	  table(i,2):=sqrt(-mse*xx2(i,i));   		      {standard error of coeff.}
	  table(i,3):=table(i,1)/table(i,2);                  {t-value                 }
	  table(i,4):=(1-IStudent (abs(table(i,3)), ncases-3))*2;   {p-value}
	  table(i,5):=table(i,1)-t*table(i,2);
	  table(i,6):=table(i,1)+t*table(i,2);
	end;
	kname$:='Interc|'+VarName (3)+'|'+VarName(2);
	shandle:=NewScrollsheet (3, 6, table, line01$, kname$, kname2$);
{highlight significant factors}
	for i:=1 to 3 do begin
	 if table(i,4)<.05 then begin
         for j:=1 to 6 do ScrollsheetSetHilite (shandle, i, j, 1);
	 end;
	end;


{Plot of observed versus residuals}
	line01$:='Observed vs. Residual Values';
	kname$:='Observed Values: '+VarName(1);
	ghandle:=NewGraph (SCATTERPLOT, line01$, 'Residuals',kname$,
	                   Ncases, y, residuals);
	GraphSetPlotFitting (ghandle, 1, FIT_LINEAR);

{Plot of case number versus residuals}
	line01$:='Case No. vs. Residual Values';
	kname$:='Case Number';
	for i:=1 to ncases do work1(i):=i;
	ghandle:=NewGraph (SCATTERPLOT, line01$, 'Residuals',kname$,
				 Ncases, work1, residuals);
	GraphSetPlotFitting (ghandle, 1, FIT_LINEAR);

2300 East 14th Street, Tulsa, OK 74104
Phone: (918) 749-1119; Fax: (918) 749-2217
e-mail: info@statsoft.com

©Copyright StatSoft, Inc., 1984-2004.
StatSoft, StatSoft logo, STATISTICA, SEWSS, SEDAS, Data Miner, SEPATH and GTrees are trademarks of StatSoft, Inc.