(* $Id: Econometrics.m,v. 5.1 2008/03/22 *)

(*:Version: Mathematica 6.0 *)

(*:Name: Econometrics *)

(*:Title: Econometric Estimation Package *)

(*:Legal:
  Copyright (c) 1987-2008, David A. Belsley
*)

(*:Summary:
  Ordinary least squares, block regression, mixed-
estimation, two-stage least squares, and instrumental-variables estimation.  BKW collinearity diagnostics.  Facility for White heteroskedasticity-
consistent estimated variance-covariance matrix.
*)

(*:Keywords: Econometrics, simultaneous systems, collinearity diagnostics  *)

(*:Requirements: No special system requirements. *)

(*:Warning: Be sure not to turn off warnings Inverse::invc and Inverse::matsq *)

(*:Sources: Econometrics texts, Principles of Econometrics, Wiley (1971) by H. Theil, and Conditioning Diagnostics: Collinearity and Weak Data in Regression Wiley (1991) by David A. Belsley. *)

History:
© David A. Belsley, 1989-2008. Version 5.1 (3/22/08).  Department of Economics,
 Boston College and Center for Computational Research in Economics and 
 Management Science, MIT.  Version 2.0 compatible.

(1.0) basic package, published in first issue of Mathematica Journal 10/88.
(1.1) adds Mixed Estimation (1/90)
(1.1.1) adds cosmetic output change, y's and x's for 2SLS and IV  (12/11/90)
(2.0) Mathematica 2.0 compatible.  Uses new formatting in varDCom (10/14/91)
(2.0.1) Includes White heteroskedasticity procedure, and uses QRDecomposition
         linked to SVD when conditioning gets very bad.  (10/16/91).
(2.0.2) Allows SVD to be chosen through "algorithm" option and cleans up calculations
        of varDCom and VarCovMatrix under different conditions. (10/18/91)
(2.1) Makes all Blocks completely modular (including checkInput[]) and cleans up all
         parameter passes.  Could use Module instead of Block, but does not since it
         runs more slowly. (10/20/91)
(2.1.1) Corrects default labeling problem when using IV[]. (3/1/92)
(2.2) Includes labeling facility, new utilty functions Name[], Call[], and UnCall[], and 
        adds labeling to the lag[] and delta[] functions.  Removes Mean[] and Correlation[]
        functions from being viewed outside package. (4/9/92)
(2.2.1) Wholesale cleaning and editing of local variates and parameters.  Extraneous
        parameters and variates removed. (4/15/92)
(2.2.2) Adds second definitions of Name[] and Call[] to allow labeling of matrices and
        tensors. (4/25/92)
(2.2.3) Editing and cleanup of code.  Addition of v/c Matrix type in display. (5/4/92)
(2.2.4) Incorporates option labels into $PreRead, and protects various symbols. (5/6/92)
(2.2.5) Further code tightening in doSVDbeta, doStats, doVarDCom, labelmaker, do2SData,
        missingObservations, checkInput, and doMixedEstData.  Fixes: handling of "const" as
        dependent variate in do2SData; calculation of InvRootSigX and xData in 
        doMixedEstData; and matrix check in checkInput. (5/14/92)
(2.2.6) Code tightening in restore[]. As published in "Economic and Financial Modeling
        with Mathematica," H. Varian, Editor, Springer-Verlag (1993) (5/20/92)
(2.2.7) Adds UnName to "usage" lists. (2/8/93)
(2.3)  Alters White HCCM to avoid nxn matrix, and applies blow-up factor of n/(n-p).
					(3/21/93)
(2.3.1)  Change For's changed to Do's in displayOutput[], checkInput[], Name[], Call[]. 
         Alter restore[] to use Indeterminate instead of "na" when filling missing observations
					(2/26/94)
(2.3.2) Fix correlation[] to return indeterminate with a constant variate argument, and
          simplify Select[] statement in missingObservations[] (8/15/94)
(2.3.3) Changed preread QR to QRFac so it would not conflict with built-in system
					symbols (2/17/95)
(2.3.4) Removes option names from PreRead and makes them visable through ::usages. 
					QRFac again becomes simply QR.  Removes FortranForm from VarDCom output
					format for singular values. (3/27/95)
(2.3.5) Adds StringTruncate to format VarDCom condition-number output in place of 
					FortranForm, etc. (6/29/95)
(2.3.6) Changes mean[] to use dot rather than Plus[]. (8/9/95)
(3.0)  Miriad changes made to make compatible with Mma. 3.0 (OutputForms needed in DisplayOut[];
		HoldForms required in Name[], lag[], delta[], Call[] (11/23/96)
(4.0) Makes 4.0 compatible; adds occasional Nulls in missing parameter places (,, goes to ,Null,); removes extra , in Protect statement
(4.0.1) Adds "Scaled" to Condition Indexes in legend of varDCom output. (3/23/03)
(4.0.2) makes minor changes to displayOut[] (D.o.F line needed shifting) to accommodate 6.0
(5.0) Replaces SingularValues[] with SingularValueDecomposition[] since the former function is to be phased out in future releases. 
           This is accomplished by using a user defined function singularValues[] (small s) that calls SingularValueDecomposition[]. (3/22/08)
(5.1) Notebook converted to 6.0 styles (code vs. input for code cells); otherwise identical to (5.0). (3/22/08)

BeginPackage["Econometrics`"]

Reg::usage = "Reg[y_,x_,opts___] runs a regression of y on
x.  y is a vector or list of vectors, x is a matrix, a list
of vectors or matrices, or a tensor of lists of vectors.
Include \"const\" (quotes are essential here) in x list if
there is to be an intercept.  y and x may contain missing
observations (any non-numeric entries).\n
Options[Reg] = {displayOutput -> True, varDCom -> False,
displayDigits -> 3, varCovMatrix -> OLS (White), 
algorithm -> QR (SVD), returnValues -> Level1 (Level2)} \n
returnValues Options
  Level1 = {beta,err,yhat,sig2}  (default)\n
  Level2 = {beta,err,yhat,sig2,se,vcMat}"

TwoStage::usage = "TwoStage[y_,yincluded_,xincluded_,
xexlcuded_,opts___] does two stage least squares.  y must
be a vector (or 1xn matrix).\n
Options same as Reg[]"

IV::usage = "IV[y_,includedVariables_,instrumentalVariables_,
opts___] does IV estimation.  If number of instVariables
equals yincluded, instVariables become instruments, else
the instruments are the yhats of includedy regressed on
includedx's and instVariables. y must be vector or 1xn
matrix.\n
Options same as Reg[]"

MixedEstimation::usage = "MixedEstimation[y_,X_,c_,R_,SigR_,
SigX_:Null,opts___] does Theil's mixed-estimation.  Basic 
model is y = Xb + e with constraints c = Rb + v.  Var(e) = 
SigX and Var(v) = SigR.  SigX defaults to sig2*I, where sig2
is the estimated variance of e from the OLS estimation of the
basic model. If R is rectangular, it may be entered with any
orientation. If it is square it must be r x p, where r is the
number of constraints and p the number of variates. All 
standard Reg opts apply."

lag::usage = "lag[list_,n_,m_] creates lags (leads) of list.
Negative values denote leads. lag[list] is list lagged once.
lag[list,n] is list lagged n times. lag[list,n,m] creates a
matrix of list lagged n through m>n times. \n\n
Options[labels -> True, lagValue -> Indeterminate]. \n
Attributes: HoldFirst."

delta::usage = "delta[x_,n_] forms x-lag[x,n], for the vector x.\n\n
Options[labels -> True]. \n
Attributes: HoldFirst."

Name::usage = "Name[x_] assigns the Mma. symbol name of x as 
its label. If x is a matrix or tensor, sequential labels based
on the symbol x are attached.\n
Attributes: HoldFirst."

UnName::usage = "UnName[x_] strips a label, if any, from the 
beginning of a data series."

Call::usage = "Call[x_,name__] labels data array x with name(s).
Do not include prefex underscore in name. x may be a vector,
matrix,or tensor of data series. There must be the same
number of names as data series in x. The names need not be
enclosed in quotes. \n
Attributes: HoldRest."

displayOutput::usage = "Option: default true. False supresses
screen display."

varDCom::usage = "Option: false.  True causes pie matrix of
conditioning diagnostics to be displayed"

returnValues::usage = "Option: default Level1."

Level1::usage = "Option for returnValues giving beta, err, yhat, 
and sig2"

Level2::usage = "Option for returnValues giving beta, err, yhat,
sig2, se, and vcMat"

displayDigits::usage = "Option: default 3. Number of digits
displayed beyond decimal point in screen display."

varCovMatrix::usage = "Option for chosing estimator for the 
variance/covariance matrix: default OLS."

White::usage = "Option for heterokedasticity consistent esitmator
of varCovMatrix"

OLS::usage = "Option for usual estimator of varCovMatrix"

algorithm::usage = "Option for chosing matrix decomposition 
procedure used for solving LS problem: default QR. "

QR::usage = "Option causing algorithm to use QR decomposition."

SVD::usage = "Option causing algorith to use SVD"

labels::usage = "Option: default true.  False prevents lag[] and 
delta[] from applying labels to variates."

lagValue::usage = "Option: default Indeterminate.  Allows lag
value to be set in lag[] and delta[]"

(*$PreRead = StringReplace[#,{"@" -> "//Name"}]&;*)

Off[NumberForm::sigz,General::spell,General::spell1];

Begin["`Private`"]

Reg[y_,x_,opts___] :=
Block[{beta,sig2,xx,se,n,g,p,missObs,err,yhat,ydata,xdata,
				retainedCases,temp,xlabels,ylabels},
	
	If[TrueQ[!(temp=doRegData[y,x])],
		Return[False],
		{ydata,xdata,n,p,g,missObs,retainedCases,ylabels,
											xlabels} = temp
	];
	
	If[TrueQ[!(temp = doRegCalculations[ydata,xdata,n,p,g,missObs,
							retainedCases,ylabels,xlabels,opts])],
		Return[False],
		{beta,err,yhat,sig2,se,xx} = temp
	];

	If[TrueQ[missObs!=0],
		err = Map[restore[#,retainedCases,n,missObs]&,err];
		yhat = Map[restore[#,retainedCases,n,missObs]&,yhat] 
	];
								
	doReturn[opts]	(* must not have semicolon *)
]

Reg::vecx = "Independent variables are not of proper form";
Reg::vecy = "Dependent variable(s) not of proper form";
Reg::doff = "Inadequate degrees of freedom";
Reg::con2 = "Constant used twice";
Reg::coll = "X matrix is essentially perfectly collinear; using SVD.";
Reg::ggt1 = "Only 1 dependent variable allowed with White varCovMatrix";

Options[Reg] := {displayOutput -> True, varDCom -> False,
	returnValues -> Level1, displayDigits->3, 
	varCovMatrix -> OLS, algorithm -> QR};
	
Options[doBeta] := {output -> 2};
	
Options[SimEqns] := {estimator -> TwoStage};

doReturn[opts___] :=
	Switch[returnValues/.{opts}/.Options[Reg],
		Level2,
			Return[Map[reduceRank,{beta,err,yhat,sig2,se,xx}]],
		_,
			Return[Map[reduceRank,{beta,err,yhat,sig2}]]
	]

doRegCalculations[yMat_,xMat_,n_,p_,g_,missObs_,retainedCases_,
									ylabels_,xlabels_,opts___] :=
Block[{beta,psinvX,r,err,yhat,xx,sig2,se,rsqr,rbar2,runcen},

	{beta,psinvX,r} = doBeta[yMat,xMat,opts];
	{err,yhat} = doErr[yMat,xMat,beta];
	
	If[TrueQ[!(xx = doVarCovMatrix[r,psinvX,err,g,opts])],
		Return[False]
	];

	{sig2,se,rsqr,rbar2,runcen} = 
		doStats[yMat,yhat,err,xx,n,p,g,opts];
	
	If[displayOutput/.{opts}/.Options[Reg],
		displayOut[beta,se,err,sig2,n,p,g,missObs,retainedCases,
						rsqr,rbar2,runcen,ylabels,xlabels,opts]
	];  
	
	If[varDCom/.{opts}/.Options[Reg],
		doVarDCom[r,p,xlabels]
	];
	
	Return[{beta,err,yhat,sig2,se,xx}]
		(* Note: all values returned here are matrices -- even
			if row matrices -- not vectors. reduceRank is used
			at doReturn to reduce ranks if needed. *)
]

doBeta[yMat_,xMat_,opts___] :=
Block[{q,r,beta},
	
	Switch[algorithm/.{opts}/.Options[Reg],
		SVD,
			{r,q} = doSVDbeta[xMat,opts],
		_,	
			{q,r} = QRDecomposition[Transpose[xMat]];
			If[TrueQ[Check[q = Inverse[r].q,True,Inverse::invc,
						Inverse::matsq]],
				Message[Reg::coll];
				{r,q} = doSVDbeta[xMat,opts]
			]
	];		
	beta = yMat.Transpose[q];
	
	Switch[output/.{opts}/.Options[doBeta],
		1,
			Return[{beta,{},{}}],
		2,
			Return[{beta,q,r}]
	] 
]
(* 
Returns:
If SVD used:
	q is unscaled pseudo inverse of xMat.
	r = {v,sv,normx}, where v and sv are scaled if varDCom is
	   opted, unscaled otherwise, and normx is the scaling (set
	   to {} when there is no scaling).
	   
If QR used:
	q is unscaled pseudo inverse of xMat,
	r is unscaled R matrix from QR.
	
If QR cannot return an R matrix of full rank, a message is given,
and the SVD procedure is automatically invoked.
	
Uses:
	q is fed to doVarCovMatrix[] for computation of "White"
	  VC Matrix when opted.
	r is fed to doVarCovMatrix[] for computation of "OLS"
	  VC Matrix.
	  * when r contains SVD triple, doVarCovMatrix[] rescales
	    to obtain proper VC matrix.
	r is also fed to doVarDCom[] for collinearity diagnostics.
	  * when r is R matrix, it is scaled there and used in SVD.
	  * when r is SVD triple, v and sv are already scaled and
	    are used directly in calculating collinearity diagnostics.
*)

doSVDbeta[xMat_,opts___] :=
Block[{u,sv,v,normx,S,scaling = varDCom/.{opts}/.Options[Reg]},

	If[TrueQ[scaling],
		normx = Map[EuclideanNorm,xMat];
		S = DiagonalMatrix[1.0/normx],
		normx = {}
	];
	
	{v,sv,u} = singularValues[
					If[TrueQ[scaling],
						xMat/normx,
						xMat
					],
					tolerance -> 10.^-30];
				
	u = Transpose[v].DiagonalMatrix[1/sv].u;  (* pseudo inv *)
	
	If[TrueQ[scaling],
		u =	S.u
	];
	
	Return[{{v,sv,normx},u}] (* note: u is the unscaled pseudo inverse now,
		but v and sv are "scaled" if varDCom has been chosen *)
]

doErr[yMat_,xMat_,beta_] :=
Block[{yhat,err},
	yhat = beta.xMat;
	err = yMat-yhat;
	Return[{err,yhat}]
]

doVarCovMatrix[r_,psinvX_,err_,g_,opts___] :=

(* note: r can be either the R matrix if QR was used in doBeta[],
or the tensor {v,sv,normx} if the SVD was used. See note 
accompanying doBeta[] for further information *)

Block[{xx,S,p,n},
	Switch[varCovMatrix/.{opts}/.Options[Reg],
		White,
			If[TrueQ[g!=1],
				Message[Reg::ggt1];
				Return[False],
				{p,n} = Dimensions[psinvX];
				p = n/(n-p);
				xx = p*psinvX.(reduceRank[err]^2*Transpose[psinvX])
			],
		_,	
			If[!MatrixQ[r],
				If[Length[r[[3]]]>0,
					S = DiagonalMatrix[1.0/r[[3]]];
					xx = S.Transpose[r[[1]]].
							DiagonalMatrix[1/r[[2]]^2].r[[1]].S,
					xx = Transpose[r[[1]]].
							DiagonalMatrix[1/r[[2]]^2].r[[1]]
				],
				xx = Inverse[r];
				xx = xx.Transpose[xx]
			]
	];
	Return[xx]
]

doStats[yMat_,yhat_,err_,xx_,n_,p_,g_,opts___] :=
Block[{i,sigEE,sig2,se,rsqr,rbar2,runcen},
	sigEE = err.Transpose[err]/(n-p);		
	sig2 = Table[sigEE[[i,i]],{i,g}];
	se = Sqrt[Table[xx[[i,i]],{i,p}]];
	Switch[varCovMatrix/.{opts}/.Options[Reg],
		White,
			se = {se},
		_,
			se = DiagonalMatrix[Sqrt[sig2]].Table[se,{g}]
	];
	rsqr = Table[correlation[yhat[[i]],yMat[[i]]]^2,{i,g}];
	rbar2 = 1.0-((n-1)/(n-p))(1.0-rsqr);
	runcen = Table[yhat[[i]].yhat[[i]]/yMat[[i]].yMat[[i]],{i,g}];
	Return[{sig2,se,rsqr,rbar2,runcen}]
]

displayOut[beta_,se_,err_,sig2_,n_,p_,g_,missObs_,retainedCases_,
			rsqr_,rbar2_,runcen_,ylabels_,xlabels_,opts___] :=
Block[{i,dw,er = err},
	Do[ Print[StringJoin["Dependent variable is ",ylabels[[i]]] ];
		If[TrueQ[missObs!=0],
			er = Map[restore[#,retainedCases,n,missObs]&,er];
			dw = Apply[Plus,Select[(Drop[er[[i]],1]-
				Drop[er[[i]],-1])^2,NumberQ]]/Apply[Plus,
				Select[er[[i]]^2,NumberQ]],
			dw = Apply[Plus,(Drop[er[[i]],1]-
				Drop[er[[i]],-1])^2]/((n-p)*sig2[[i]])
		];
	
		Print[StringForm["RSquared = ``  RBarSquared = ``",
			rsqr[[i]],rbar2[[i]]]];
		Print[StringForm["R2uncentered = ``  SER = ``",
			runcen[[i]],Sqrt[sig2[[i]]] ]];
		Print[StringForm["Num of Observations = ``  Degrees of Freedom = ``",n,n-p] ];
		Print[StringForm["dw = `` with `` missing obs.",
			dw, missObs] ];
		If[(varCovMatrix/.{opts}/.Options[Reg])==White,
			Print["Variance-Covariance Matrix: White"] ];
		Print[" "];
		
		Print[PaddedForm[TableForm[{beta[[i]],se[[i]],
			beta[[i]]/se[[i]]},TableDirections->{Row,Column},
			TableSpacing->{0,0},TableHeadings->{{"       coef.",
			"     st. err.","         t"},xlabels}],
			{10,displayDigits/.{opts}/.Options[Reg]}]];		
		Print[" "];
		If[g<Length[ylabels], (* only if 2SLS or IV *)
			Print["Instruments:",Drop[ylabels,1]]
		];
		Print[" "],
	{i,g}]
]

doVarDCom[r_,p_,labels_] :=

(* note: this function can receive its parameter r in two forms:
r can either be the R matrix from the QR decomposition if this
was used in doBeta[], or the tensor {v,sv,normx} if the SVD was 
used. It is first determined below which of these prevails. 

See note accompanying doBeta[] for further information *)

Block[{temp,sv,v,vMat,condind},
	If[!MatrixQ[r],
		v = r[[1]];
		sv = r[[2]],
		temp = Transpose[r];
		temp = temp/Map[EuclideanNorm,temp];
		{v,sv,temp} = singularValues[temp,tolerance->10^-30]
	];

	vMat = (Transpose[v].DiagonalMatrix[1.0/sv])^2;
	vMat = Truncate[Map[(#/Apply[Plus,#])&,vMat],4];
	
(* Form condition indexes and output matrix with
	condition indexes in first row and the variance-
	decomposition proportions in the next p rows *)
		 		
 	PrependTo[vMat,condind = Max[sv]/sv];
	vMat = Transpose[Drop[Transpose[
					Reverse[Sort[Transpose[vMat]]]],1]];
	condind = StringTruncate[Reverse[Sort[condind]]];

	Print[" "];
	Print["Variance-decomposition proportions"];
	Print[" "];
	Print["           Scaled condition indexes"];
	Print[" "];
	
	Print[PaddedForm[TableForm[vMat,TableSpacing->{1,0},
		TableHeadings->{condind,labels},
		TableDirections->{Row,Column}],{4,3},
		NumberPadding->{" ","0"}]]
]

doRegData[yblock_,xblock_] :=
Block[{ydata,xdata,n,p,g,xCon,yCon,missObs,retainedCases,
	temp,ny,nx,ylabels,xlabels},
	If[TrueQ[!(temp = checkInput[yblock])],
		Message[Reg::vecy];
		Return[False],
		{ydata,ny,g,yCon,ylabels} = temp
	];
	
	If[TrueQ[!(temp = checkInput[xblock])],
		Message[Reg::vecx];
		Return[False],
		{xdata,nx,p,xCon,xlabels} = temp
	];
	
	If[yCon && xCon,
		Message[Reg::con2];
		Return[False]
	];

	n = Min[nx,ny];
	
	If[yCon,
		PrependTo[ydata,Table[1.0, {n}]];
		PrependTo[ylabels,"Const"]
	];
	If[xCon,
		PrependTo[xdata,Table[1.0, {n}]];
		PrependTo[xlabels,"Const"]
	];
	
	ylabels = labelMaker[ylabels,"y"];
	xlabels = labelMaker[xlabels,"X"];

	If[TrueQ[!(temp = missingObservations[Join[ydata,xdata],n,p])],
		Message[Reg::doff];
		Return[False],
		ydata = Take[temp[[1]],{1,g}];
		xdata = Drop[temp[[1]],{1,g}];
		{retainedCases,missObs,n} = Drop[temp,1]
	];
	
	Return[{ydata,xdata,n,p,g,missObs,retainedCases,ylabels,xlabels}]
]

checkInput[dataBlock_] :=

(*  Recursive routine to build a data matrix from an input of 
any combination of lists, matrices or irregular arrays.  On 
normal return, "data" is possibly an irregular array, nobs is 
the length of the shortest row in data, varCount is the number 
of rows, and Con is set True if a "const" is encountered and 
otherwise False.   checkInput returns False if dataBlock has
Depth 1 or if two "constants" are encountered. *)

Block[{varCount=0,data={},Con=False,nobs=Infinity,temp,i,
			labels={},tempLabels,tempNobs},

	If[MemberQ[{"const","Const","constant","Constant"},dataBlock],
		varCount++;
		Con = True;
		Return[{data,nobs,varCount,Con,labels}]
	];
	
	If[TrueQ[Depth[dataBlock] == 1],
		Return[False]  
	];
	
	If[VectorQ[dataBlock],
		{labels,data,nobs} = labelStrip[dataBlock];
		varCount++;
		Return[{{data},nobs,varCount,Con,{labels}}]
	];
	
	If[MatrixQ[dataBlock],
		{nobs,varCount} = Dimensions[dataBlock];
		If[nobs < varCount,
			{nobs,varCount} = {varCount,nobs};
			data = dataBlock,
			data = Transpose[dataBlock]
		];
		Do[ {tempLabels,data[[i]],tempNobs}=labelStrip[ data[[i]] ];
			nobs = Min[nobs,tempNobs];
			AppendTo[labels,tempLabels],
		{i,varCount}];
		Return[{data,nobs,varCount,Con,labels}]
	];
	
	Do[ If[TrueQ[!(temp = checkInput[dataBlock[[i]]])],
			Return[False],
			If[Con&&temp[[4]],
				Message[Reg::con2];
				Return[False],
				Con = Con||temp[[4]]
			];
			data = Join[data,temp[[1]]];
			nobs = Min[nobs,temp[[2]]];
			varCount += temp[[3]];
			labels = Join[labels,temp[[5]]]
		],
	{i,Length[dataBlock]}];
	
	Return[{data,nobs,varCount,Con,labels}]
]

labelStrip[dataVector_] :=
Block[{firstPiece,underscore,nobs},
	firstPiece = ToString[dataVector[[1]]];
	underscore = StringPosition[firstPiece,"_",1];
	nobs=Length[dataVector];

	If[underscore=={{1,1}},
		Return[{StringDrop[firstPiece,1],Drop[dataVector,1],nobs-1}]
	];
	
	If[Length[underscore]>0,
		Return[{firstPiece,Drop[dataVector,1],nobs-1}]
	];
	
	Return[{"**",dataVector,nobs}]
]

labelMaker[labels_,character_] :=
	 Table[	If[labels[[i]]=="**",
				StringJoin[character,"[",ToString[i],"]"],
				labels[[i]]
			],
			{i,Length[labels]}]

missingObservations[data_,n_,p_] :=
Block[{tempData,numMissing,nobs,retainedCases},
	tempData = Map[Take[#,n]&,data]//N;
	PrependTo[tempData,Range[n]];
	tempData = Transpose[tempData];
	
	tempData = Select[tempData,VectorQ[#,NumberQ]&]; (* remove non-num'l obs *)

	nobs = Length[tempData];
	If[nobs < p, Return[False] ]; (* negative degrees of freedom *)
	numMissing = n - nobs;
	tempData = Transpose[tempData];
	retainedCases = tempData[[1]];
	tempData = Drop[tempData,1]; (* data matrix suitable for regr. *)
	
	Return[{tempData,retainedCases,numMissing,nobs}]
]

restore[series_,index_,n_,misObs_] :=
Block[{temp},
	temp = Table[Indeterminate,{n+misObs}];
	Map[(temp[[ index[[#]] ]] = series[[#]])&,
		Range[ Length[index] ] ];
	Return[temp]
]

IV[y_,includedVariates_,instrumentalVariables_,opts___]:=
	TwoStage[y,includedVariates,{},instrumentalVariables,opts,
		estimator -> IV]

TwoStage[y_,includedy_,includedx_,excludedx_,opts___] :=
Block[{yhat,beta,xx,ydata,y1data,x1data,x2data,n,g,g1,k1,k2,
	Con,retainedCases,err,missObs,sig2,se,rsqr,rbar2,runcen,
	temp,psinv,r,ylabels,y1labels,x1labels,x2labels},
	
	(* Prepare data *)
	If[TrueQ[!(temp=do2SData[y,includedy,includedx,excludedx,opts])],
		Return[False],
		{ydata,y1data,x1data,x2data,n,g,g1,k1,k2,Con,retainedCases,
		missObs,ylabels,y1labels,x1labels,x2labels} = temp
	];
	
	Switch[estimator/.{opts}/.Options[SimEqns],
		IV,
			If[k2>g1,
				temp = doBeta[y1data,x2data,opts,output->1];
				x2data = temp[[1]].x2data
			];
			{beta,xx,err,yhat} = doIV[ydata,y1data,x2data,opts];
			x1labels = y1labels;
			ylabels = Join[ylabels, x2labels],
			
		_,  	(* Default is TwoStage  *)
				(* form first stage *)
			temp = doBeta[y1data,Join[x1data,x2data],opts,output->1];
			yhat = temp[[1]].Join[x1data,x2data];
			
			ylabels = Join[ylabels,x2labels];
				(* Adjust constant to first row if present *)
			If[Con,
				PrependTo[yhat,x1data[[1]]];
				PrependTo[y1data,x1data[[1]]];
				x1data = Drop[x1data,1];
				x1labels = Join[Take[x1labels,1],y1labels,Drop[x1labels,1]],
				x1labels = Join[y1labels,x1labels]
			];
			
				(* Form 2nd stage estimator using 1st stage yhat  *)
			{beta,psinv,r} = doBeta[ydata,Join[yhat,x1data],opts];
			
				(* Form 2nd stage err and yhat *)
			{err,yhat} = doErr[ydata,Join[y1data,x1data],beta];
			If[TrueQ[!(xx = doVarCovMatrix[r,psinv,err,g,opts])],
				Return[False]
			]
	];
	

	(* second-stage statistics *)
	{sig2,se,rsqr,rbar2,runcen} = 
		doStats[ydata,yhat,err,xx,n,g1+k1,g,opts];
		
	If[displayOutput/.{opts}/.Options[Reg],
		displayOut[beta,se,err,sig2,n,k1+g1,g,missObs,retainedCases,
							rsqr,rbar2,runcen,ylabels,x1labels,opts]
	];

	If[varDCom/.{opts}/.Options[Reg],
		doVarDComIV[xx,reduceRank[se],reduceRank[sig2],g1+k1,x1labels]
	];
	
	If[TrueQ[missObs!=0],
		err = Map[restore[#,retainedCases,n,missObs]&,err];
		yhat = Map[restore[#,retainedCases,n,missObs]&,yhat] 
	];
	
	doReturn[opts] (* must not have semicolon *)
]

doVarDComIV[xx_,se_,sig2_,p_,labels_] :=
Block[{v,sv,S,vMat},
	S = DiagonalMatrix[1.0/se];
	vMat = S.(sig2*xx).S;
	{v,sv,vMat} = singularValues[vMat,tolerance -> 10.^-30];
	doVarDCom[{v,Sqrt[sv]},p,labels]
]

do2SData[y_,y1_,x1_,x2_,opts___] :=
Block[{ydata,y1data,x1data,x2data,n,g,g1,k1,k2,Con,yCon,y1Con,
	x1Con,x2Con,temp,ny1,nx1,nx2,retainedCases,missObs,ylabels,
	y1labels,x1labels,x2labels},
	
	If[TrueQ[!(temp = checkInput[y])],
		Message[TwoStage::ydat];
		Return[False],
		{ydata,n,g,yCon,ylabels} = temp
	];
	
	If[g!=1,
		Message[TwoStage::ydat];
		Return[False]
	];
	
	If[TrueQ[!(temp = checkInput[y1])],
		Message[TwoStage::y1dat];
		Return[False],
		{y1data,ny1,g1,y1Con,y1labels} = temp
	];
	
	If[TrueQ[!(temp = checkInput[x2])],
		Message[TwoStage::x2dat];
		Return[False],
		{x2data,nx2,k2,x2Con,x2labels} = temp
	];
	
	If[k2<g1,
		Message[TwoStage::unid];
		Return[False]
	];
	
	Switch[estimator/.{opts}/.Options[SimEqns],
		IV,
			x1data = {};
			nx1 = Infinity;
			k1 = 0;
			If[(yCon&&y1Con)||(yCon&&x2Con),
				Message[TwoStage::con2];
				Return[False]
			];
			x1Con = False;
			Con = y1Con,
		_,
			If[TrueQ[!(temp = checkInput[x1])],
				Message[TwoStage::x1dat];
				Return[False],
				{x1data,nx1,k1,x1Con,x1labels} = temp
			];
			If[Count[{yCon,y1Con,x1Con,x2Con},True] > 1,
				Message[TwoStage::con2];
				Return[False]
			];
			Con = x1Con
	];

	n = Min[n,ny1,nx1,nx2];
	
	If[yCon,
		PrependTo[ydata,Table[1.0, {n}]];
		PrependTo[ylabels,"Const"]
	];
	If[y1Con,
		PrependTo[y1data,Table[1.0, {n}]];
		PrependTo[y1labels,"Const"]
	];
	If[x1Con,
		PrependTo[x1data,Table[1.0, {n}]];
		PrependTo[x1labels,"Const"]
	];
	If[x2Con,
		PrependTo[x2data,Table[1.0, {n}]];
		PrependTo[x2labels,"Const"]
	];
	If[TrueQ[!(temp = missingObservations[Join[ydata,
								y1data,x1data,x2data],n,k1+k2])],
			Message[TwoStage::doff];
			Return[False],
			{temp,retainedCases,missObs,n} = temp
	];
	
	ydata  = Take[temp, {1,1}];
	y1data = Take[temp, {2,g1+1}];
	x1data = Take[temp, {g1+2,g1+1+k1}];
	x2data = Take[temp, {-k2,-1}];
		
	ylabels  = labelMaker[ylabels,"y"];
	y1labels = labelMaker[y1labels,"Y"];
	x1labels = labelMaker[x1labels,"X"];
	x2labels = labelMaker[x2labels,"Z"];
	
	Return[{ydata,y1data,x1data,x2data,n,g,g1,k1,k2,Con,
		retainedCases,missObs,ylabels,y1labels,x1labels,x2labels}]
]

doIV[y_,z_,w_,opts___] :=
Block[{delta,zwInv,psinvZW,xx,err,yhat},
	zwInv = Inverse[z.Transpose[w]];
	psinvZW = Transpose[zwInv].w;
	delta = y.Transpose[psinvZW];
	yhat = delta.z;
	err = y-yhat;
	Switch[varCovMatrix/.{opts}/.Options[Reg],
		White,
			xx = psinvZW.DiagonalMatrix[reduceRank[err]^2]
				.Transpose[psinvZW],
		_,
			xx = psinvZW.Transpose[psinvZW]
	];
				
	Return[{delta,xx,err,yhat}]
]

TwoStage::ydat := "y variate incorrect."
TwoStage::y1dat := "Included variables or included
endogenous variates incorrect."
TwoStage::x1dat := "Included exogenous variates incorrect."
TwoStage::x2dat := "Instrumental variables or excluded
exogenous variates incorrect."
TwoStage::unid := "Inadequate instruments -- underidentified."
TwoStage::con2 := "Constant used twice."
TwoStage::doff := "Inadequate degrees of freedom."
MixedEstimation::ydat := "Dependent variate must be a vector."
MixedEstimation::cdat := "c must be a vector of length r."
MixedEstimation::Rdat := "R must be a matrix of size r x p."
MixedEstimation::sigR := "SigR must be a symmetric matrix of size r."
MixedEstimation::sigX := "SigX must be a symmetric matrix of size n."

MixedEstimation[y_,X_,c_,R_,SigR_,SigX_:Null,opts___] :=
Block[{beta,sig2,xx,se,n,g,p,missObs,err,yhat,rsqr,rbar2,runcen,
	ydata,xdata,retainedCases,temp,s,ylabels,xlabels,r,psinvX},

	If[TrueQ[!(temp=doRegData[y,X])],
		Return[False],
		{ydata,xdata,n,p,g,missObs,retainedCases,
			ylabels,xlabels} = temp
	];
	
	If[TrueQ[!(temp=doMixedEstData[ydata,xdata,n,p,g,c,R,
											SigR,SigX,opts])],
		Return[False],
		{ydata,xdata,s}=temp (*s is the no.of prior constraints*)
	];	

	temp = n+missObs;
	retainedCases = Join[retainedCases,Range[temp+1,temp+s]];
	n = n+s;
	
	(* ready to do Mixed Estimation.*)
	
 	{beta,psinvX,r} = doBeta[ydata,xdata,opts];
 	{err,yhat} = doErr[ydata,xdata,beta];
	
	If[TrueQ[!(xx = doVarCovMatrix[r,psinvX,err,g,opts])],
		Return[False]
	];

	{sig2,se,rsqr,rbar2,runcen} = 
		doStats[ydata,yhat,err,xx,n,p,g,opts];
	
	se = se/Sqrt[reduceRank[sig2]];
	sig2 = {1.0};
	
	If[displayOutput/.{opts}/.Options[Reg],
		displayOut[beta,se,err,sig2,n,p,g,missObs,retainedCases,
							rsqr,rbar2,runcen,ylabels,xlabels,opts]
	];  
	
	If[varDCom/.{opts}/.Options[Reg],
		doVarDCom[r,p,xlabels]
	];

 	If[TrueQ[missObs!=0],
		err = Map[restore[#,retainedCases,n,missObs]&,err];
		yhat = Map[restore[#,retainedCases,n,missObs]&,yhat] 
	];
								
	doReturn[opts]	(* must not have semicolon *)
]

doMixedEstData[y_,X_,n_,p_,g_,c_,R_,SigR_,SigX_,opts___] :=
Block[{ydata,xdata,InvRootSigR,InvRootSigX,RR = R,s,sp,a,b,v,sv,u},

	If[g!=1,
		Message[MixedEstimation::ydat];
		Return[False]
	];
	
	If[!MatrixQ[RR],
		Message[MixedEstimation::Rdat];
		Return[False],
		{s,sp} = Dimensions[RR];
		If[sp!=p && s!=p,
			Message[MixedEstimation::Rdat];
			Return[False]
		];
		If[sp!=p,
			RR = Transpose[RR];
			s = sp
		]
	];  
	
	(*Note: A problem arises as to the orientation R (RR) should
	have if s = sp = p.  It is assumed in this case that s is the
	number of constraints, and sp the number of variates. *)
		
	If[!VectorQ[c] || s!=Length[c],
		Message[MixedEstimation::cdat];
		Return[False]
	];
		
	If[!SymmetricQ[SigR] || !TrueQ[Length[SigR] == s],
		Message[MixedEstimation::sigR];
		Return[False],
		{v,sv,u} = singularValues[SigR//N,tolerance -> 10.^-30];
		InvRootSigR = -DiagonalMatrix[1.0/Sqrt[sv]].u
	];
	
	If[TrueQ[SigX==Null],
		a = doBeta[y,X,opts,output->1];
		b = doErr[y,X,a[[1]]]; (* a[[1]] is "beta" *)
		b = reduceRank[b[[1]]]; (* b[[1]] is "err" *)
		b = 1.0/Sqrt[b.b/(n-p)];(* this is 1/sig2 *)
		ydata = b reduceRank[y];
		xdata = b X,
		If[!SymmetricQ[SigX] || !TrueQ[Length[SigX] == n],
			Message[MixedEstimation::sigX];
			Return[False],
			{v,sv,u} = singularValues[SigX//N,tolerance -> 10.^-30];
			InvRootSigX=-Transpose[u].DiagonalMatrix[1.0/Sqrt[sv]];
			ydata = reduceRank[y].InvRootSigX;
			xdata = X.InvRootSigX
		]
	]; 

	
	ydata = Join[ydata,InvRootSigR.c];
	xdata = Adjoin[xdata,Transpose[InvRootSigR.RR]];

	Return[{{ydata},xdata,s}]
]

(* Utility Programs for Reg *)

(*mean[list_List] := N[Apply[Plus,list]/Length[list]]*)

mean[list_] := With[{n = Length[list]},
					(Table[1.,{n}].list)/n]

correlation[x1_List,x2_List] :=
	Block[{Z = {x1,x2}},
		Z = Z-mean[Transpose[Z]];
		Z = Z.Transpose[Z];
		If[Z[[1,1]]== 0. || Z[[2,2]] == 0.,Indeterminate,
		Z[[1,2]]/Sqrt[Z[[1,1]] Z[[2,2]]]
		]
	]/; Length[x1] == Length[x2]

EuclideanNorm[list_] := N[Sqrt[list.list]]

Truncate[x_,n_] := N[Round[x 10.^n] 10.^-n]

StringTruncate[x_] :=
Block[{digits, num},
	{digits, num} = RealDigits[x+.05//N];
	digits = Map[ToString, digits];
	If[num<=5,
		If[Length[digits]>num,
			digits = Insert[digits,".",num+1];
			Return[StringJoin[Take[digits,num+2]]],
			AppendTo[digits,".0"];
			Return[StringJoin[digits]]
		],
		digits = Insert[digits,".",2];
		digits[[6]] = "e";
		digits[[7]] = ToString[num-1];
		Return[StringJoin[Take[digits,7]]]
	]
]
SetAttributes[StringTruncate,Listable]

SymmetricQ[mat_] := TrueQ[MatrixQ[mat] && mat == Transpose[mat]]

reduceRank[x_] := If[MatrixQ[x]&&(Min[Dimensions[x]]==1)||
	VectorQ[x]&&(Length[x]==1),
	Return[x[[1]]],Return[x] ]

lag[list_,n_:1,m_:Null,opts___] :=
Block[{i,r,temp,label,heldList=HoldForm[list],
										first=StringMatchQ[ToString[list[[1]]],"*_*"]},
	If[TrueQ[n==Null],r=1,r=n];
	If[TrueQ[m<r],
		Message[lag::mltn];
		Return[False]
	];
	If[!VectorQ[list],
		Message[lag::vec];
		Return[False]
	];
	If[TrueQ[labels/.{opts}/.Options[lag]],
		If[first,
			label = ToString[list[[1]]];
			temp = Drop[list,1],
			label = StringJoin["_",ToString[heldList]];
			temp = list
		],
		If[first,
			temp = Drop[list,1],
			temp = list
		];
		label = {}
	];
		
	If[TrueQ[m==Null],
		Which[
			r<0,
				If[TrueQ[labels/.{opts}/.Options[lag]],
					label = {StringJoin[label,"(",ToString[-r],")"]}
				];
				temp = Join[label,Drop[temp,-r],Table[lagValue/.{opts}/.Options[lag],{-r}]],
			r==0,
				If[TrueQ[labels/.{opts}/.Options[lag]],
						label = {label}
				];
				temp = Join[label,temp],
			r>0,
				If[TrueQ[labels/.{opts}/.Options[lag]],
						label = {StringJoin[label,"(",ToString[-r],")"]}
				];
				temp = Join[label,Table[lagValue/.{opts}/.Options[lag],{r}],Drop[temp,-r]]
		],
		temp = Table[lag[list,i,Null,opts], {i,r,m}]
	];
	Return[temp]
]

Options[lag] := {labels->True,lagValue->Indeterminate};
SetAttributes[lag,HoldFirst];
lag::mltn := "End of range before beginning.";
lag::vec := "First argument must be a vector.";

delta[x_,n_:1,opts___] := 
Block[{r,temp},
	If[TrueQ[n==Null],r=1,r=n];
	temp = lag[x,0,Null,opts] - lag[x,r,Null,opts];
	If[TrueQ[labels/.{opts}/.Options[lag]],
		temp[[1]] = StringJoin["Delta(",ToString[-r],")",ToString[temp[[1,1]]]]
	];
	Return[temp]
]
SetAttributes[delta,HoldFirst];

Adjoin[x_,y_] := 
Block[{nx,ny,m},
	Adjoin::nMat = "inputs must be matrices";
	Adjoin::wrSz = "inputs must have same number of rows";
	If[MatrixQ[x] && MatrixQ[y],
		{nx,m} = Dimensions[x];
		{ny,m} = Dimensions[y],
		Message[Adjoin::nMat];
		Return[False]
	];
	If[nx != ny,
		Message[Adjoin::wrSz];
		Return[False],
		Map[Join[x[[#]],y[[#]]]&,Range[ny]]
	]
]

(* Note: this function requires DataPrep[] and the messages
from the Call[] function *)

Name[x_] :=
	With[{heldx=HoldForm[x]},
		If[StringMatchQ[ToString[First[x]],"*_*"],
			Prepend[Drop[x,1],StringJoin["_",ToString[heldx]]],
			Prepend[x,StringJoin["_",ToString[heldx]]]
		]
	]/;VectorQ[x]
SetAttributes[Name,HoldAll];

Name[x_] :=
Block[{heldx=HoldForm[x],temp = x,p,i,base},
	If[MatrixQ[temp]&&Apply[Greater,Dimensions[temp]],
		temp = Transpose[temp]
	];
	p = Length[temp];
	base = ToString[heldx];
	Do[ If[!VectorQ[temp[[i]]],
			Message[Call::nvec,i];
			Return[False]
		];
		temp[[i]] = DataPrep[temp[[i]], 
					  StringJoin[base,"[",ToString[i],"]"]],
	{i,p}];
	Return[temp]
]

UnName[x_] :=
	If[StringMatchQ[ToString[First[x]],"*_*"],
		Drop[x,1],
		x
	]/;VectorQ[x]

Call[x_,label_] := 
	DataPrep[x,ToString[HoldForm[label]]]/;VectorQ[x]
	
Call[x_,label__] :=
Block[{data = x,p,i,commas,ids=ToString[HoldForm[{label}]]},
	p = Length[commas=Join[{{0,0}},StringPosition[ids,","],
		{{StringLength[ids],0}}]]-1;
	ids = Table[StringTake[ids,{commas[[i,1]]+2,commas[[i+1,1]]-1}],{i,p}];
	If[MatrixQ[data]&&Apply[Greater,Dimensions[data]],
		data = Transpose[data]
	];
	If[(Length[data] != p)||VectorQ[data],
		Message[Call::ulen];
		Return[False]
	];
	Do[ If[!VectorQ[data[[i]]],
			Message[Call::nvec,i];
			Return[False]
		];
		data[[i]] = DataPrep[data[[i]],ids[[i]]],
	{i,p}];
	Return[data]
]

SetAttributes[Call,HoldRest];

Call::ulen := "Number of series and names not equal.";
Call::nvec := "Data series `1` not a vector.";

DataPrep[x_,label_] :=
	If[StringMatchQ[ToString[First[x]],"*_*"],
		Prepend[Drop[x,1],StringJoin["_",label]],
		Prepend[x,StringJoin["_",label]]
	]

(* This function added to bridge the old SingularValues[], being phased out my Mma., 
with the new SingularValueDecomposition[] to allow the original function calls to 
remain transparently intact. *)

Options[singularValues] := {tolerance -> 0};

singularValues[x_,opts___] := 
Block[{k,u,d,v},
	k = Min[Dimensions[x]];
	{u,d,v} = SingularValueDecomposition[x,k,Tolerance -> tolerance/.{opts}/.Options[singularValues]];
	d = Map[d[[#,#]]&,Range[k]];
	Return[{Transpose[u],d,Transpose[v]}]
]

Protect[displayOutput,varDCom,returnValues,Level1,
	Level2,displayDigits,varCovMatrix,White,OLS,
	algorithm,QR,SVD,labels,lagValue];

End[];

Protect[Reg,TwoStage,IV,MixedEstimation,lag,delta,
	UnName,Name,Call];

EndPackage[]