(*:Version: Mathematica 1.2 *)

(*:Name: Econometrics *)

(*:ID:            *)

(*:Title: Econometrics.m *)

(*:Author: David A. Belsley, Department of Economics, Boston College and Center for Computational Research in Economics and Management Science, MIT, 1990.*)

(*:Legal:          *)

(*:Reference:      *)

(* Copyright, David A. Belsley, 1989. Version 1.0. *)

BeginPackage["Econometrics`"]
Reg::usage = "Reg[y_,x_,opts___] runs a regression of y on
x.  y is a vector, x is a matrix or a list of vectors,
matrices or 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).\r \r
Options[Reg] = {displayOutput -> True, varDCom ->False,
returnValues -> Level1}"

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

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."

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.  This version
preserves the length of the data series, adding 'na's as needed."

Mean::usage = "Mean[list_List], the arithmetic mean of list."

Correlation::usage = "Correlation[x1_List, x2_List], gives the
correlation between any two lists of the same length."



(* Aliases make simpler names globally accessable  *)
displayOutput ::= Econometrics`private`displayOutput
varDCom ::= Econometrics`private`varDCom
returnValues ::= Econometrics`private`returnValues

Begin["`private`"]

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

	If[TrueQ[missObs!=0],
		err = Map[restore[#,retainedCases]&,err];
		yhat = Map[restore[#,retainedCases]&,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";

Options[Reg] := {displayOutput -> True, varDCom -> False,
	returnValues -> Level1};
	
Options[SimEqns] := {estimator -> TwoStage}
doReturn[opts___] :=
	Switch[ToString[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_,opts___] :=
Block[{v,sv,yhat,beta,xx,
		err,sigEE,sig2,se,rsqr,rbar2,runcen},

	{v,sv,beta,xx} = doReg[yMat,xMat,opts];
	
	{yhat,err,sigEE,sig2,se,rsqr,rbar2,runcen} = 
		doStats[yMat,xMat,beta,xx,n,p,g];
	
	If[displayOutput/.{opts}/.Options[Reg],
		displayOut[beta,se,err,sig2,xCon,n,p,g,missObs,
			retainedCases,rsqr,rbar2,runcen]
	];  
	
	If[varDCom/.{opts}/.Options[Reg],
		doVarDCom[v,sv,p]
	];
	
	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. *)
]

doReg[yMat_,xMat_,opts___] :=
Block[{u,sv,v,beta,xx,normx,S,
			scaling = varDCom/.{opts}/.Options[Reg]},
	
	If[TrueQ[scaling],
		normx = Map[EuclideanNorm,xMat];
		S = DiagonalMatrix[1.0/normx]
	];
	
	{v,sv,u} = SingularValues[
					If[TrueQ[scaling],
						xMat/normx,
						xMat
					],
				Tolerance -> 10.^-30];
				
	u = Transpose[v].DiagonalMatrix[1/sv].u;(* pseudo inv *)
	
	beta = If[TrueQ[scaling],
				yMat.Transpose[u].S,
				yMat.Transpose[u]
			];
			
	xx = If[TrueQ[scaling],
			S.Transpose[v].DiagonalMatrix[1/sv^2].v.S,
			Transpose[v].DiagonalMatrix[1/sv^2].v
		];
	Return[{v,sv,beta,xx}]
]

doStats[yMat_,xMat_,beta_,xx_,n_,p_,g_] :=
Block[{i,temp,yhat,err,sigEE,sig2,se,rsqr,rbar2,runcen},
	yhat = beta.xMat;
	err = yMat-yhat;
	sigEE = err.Transpose[err]/(n-p);		
	sig2 = Table[sigEE[[i,i]],{i,g}];
	temp = Sqrt[Table[xx[[i,i]],{i,p}]];
	se = DiagonalMatrix[Sqrt[sig2]].Table[temp,{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[{yhat,err,sigEE,sig2,se,rsqr,rbar2,runcen}]
]

displayOut[beta_,se_,err_,sig2_,xCon_,n_,p_,g_,missObs_,
		retainedCases_,rsqr_,rbar2_,runcen_] :=
Block[{i,j,f,labels=Array[f,p],er = err},
	For[i=1,i<=p,i++,
		labels[[i]] = StringJoin["X[",ToString[i],"]"]
	];
	If[xCon,labels[[1]] = "Const"];
	For[i=1,i<=g,i++,
		Print[StringForm["Dependent variable ``",i] ];
		If[TrueQ[missObs!=0],
			er = Map[restore[#,retainedCases]&,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] ];
		Print[" "];
		Print[StringForm["var \t\tcoef.\t\tst. err.\t\tt"]];
		Print[" "];
		
		For[j=1,j<=p,j++,
			Print[StringForm["``\t\t``\t\t``\t\t``",
				labels[[j]], N[beta[[i,j]],5],
					N[se[[i,j]],5],
					N[beta[[i,j]]/se[[i,j]],5] ] ] ];
		Print[" "];
		Print[" "]
	]
]

doVarDCom[v_,sv_,p_] :=
Block[{vMat,i,j,temp,tempStr},	

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

	Print[" "];
	Print["Variance-decomposition proportions"];
	Print[" "];
	Print["                  Condition indexes"];
	Print[" "];
	
	
	tempStr = "    ";
	For[i=1,i<=p,i++,
		tempStr = StringJoin[tempStr,"\t",ToString[FortranForm
			[N[vMat[[1,i]],3]]]]
	];
	Print[tempStr];
	
	Print["  "];
	For[i=1,i<=p,i++,
		tempStr = ToString[StringForm["V[``]",i]];
		For[j=1,j<=p,j++,
			tempStr = StringJoin[tempStr,"\t",ToString[StringForm["``",
				If[(temp = vMat[[i+1,j]])<1.0 && temp > 0.0, temp, 
				If[temp == 0.0, "0.000","1.000"]]] ]]
		];
		Print[tempStr];
	];
]

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

	n = Min[nx,ny];
	
	If[yCon,PrependTo[ydata,Table[1.0, {n}]] ];
	If[xCon,PrependTo[xdata,Table[1.0, {n}]] ];

	If[TrueQ[!(temp=missingObservations[Join[ydata,xdata],n,p])],
		Message[MessageName[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,yCon,xCon,missObs,retainedCases}]
]

checkInput[dataBlock_] :=

(*  Produces a data matrix from an input of any combination 
of lists, matrices or irregular arrays.  Sets globals for 
checkParts[] which is called recursively.  On return, data is 
the possibly irregular data 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 *)

Block[{varCount=0,data={},Con=False,nobs=Infinity},
	If[!checkParts[dataBlock], Return[False] ];
			
	Return[{data,nobs,varCount,Con}]
]

checkParts[el_] :=

(*  Recursive routine to coalesce divergent input forms into
a single, possibly irregular array.  Globals set by calling
routine: Con,data,nobs, and varCount. nobs is set to minimum 
of all data lengths.  Con is set True if a "const" is included 
among el. Returns True unless an empty series is encountered.*)
 
Block[{nt,pt,i},

	If[MemberQ[{"const","Const","constant","Constant"},el],
		If[Con,
			Message[MessageName[Reg,"con2"]];
			Return[False],
			varCount++;
			Con = True;
			Return[True]
		]
	];
	If[TrueQ[Depth[el] == 1],
		Return[False]  
	];
	If[VectorQ[el],
		nobs = Min[nobs,Length[el]];
		varCount++;
		AppendTo[data,el];
		Return[True]
	];
	If[MatrixQ[el]&&(Depth[el]==3), 
(* This second condition is necessary only because
MatrixQ[] works incorrectly  *)
		{nt,pt} = Dimensions[el];
		If[nt < pt,
			{nt,pt} = {pt,nt};
			data = Join[data,el],
			data = Join[data,Transpose[el]]
		];
		nobs = Min[nobs, nt];
		varCount += pt;
		Return[True]
	]; 
	For[i=1,i<=Length[el],i++,
		If[!checkParts[el[[i]]],
			Return[False]
		]
	];
	Return[True]
]

missingObservations[dataBlock_,blockLength_,testLength_] :=

(* routine to remove incomplete observations - ones with
non-numerical entries.  Globals: n, p, xdata, ydata, 
retainedCases,missObs. At end, xdata and ydata contain data 
suitable for regression, and caseRecord contains indexes of
retained observations. Returns True unless degrees of freedom 
become zero or below. *)

Block[{i=1,j,tempData,missingCount=0,nobs=blockLength,
				retainedCases},
	tempData = Map[Take[#,nobs]&,dataBlock]//N;
	PrependTo[tempData,Table[j,{j,nobs}]];
	tempData = Transpose[tempData];
	
	While[i<=nobs,
		If[TrueQ[Scan[If[!NumberQ[#], 
				Return[False]]&,tempData[[i]] ]==False],
			tempData = Drop[tempData,{i,i}];
			missingCount++;
			nobs--,
			i++
		];
		If[nobs<=testLength,
			Return[False]
		]			
	];
	tempData = Transpose[tempData];
		
	retainedCases = tempData[[1]];
	tempData = Drop[tempData,1];
	Return[{tempData,retainedCases,missingCount,nobs}]
] 

restore[series_,index_] :=
Block[{temp,i},
	temp = Table["na",{index[[-1]]}];
	For[i=1,i<=Length[index],i++,
		temp[[index[[i]]]] = series[[i]]
	];
	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,sigEE,rsqr,rbar2,
	runcen,temp},
	
	(* 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} = temp;
	];
	
	Switch[estimator/.{opts}/.Options[SimEqns],
		IV,
			If[k2>g1,
				temp = doReg[y1data,x2data];
				x2data = temp[[3]].x2data
			];
			{beta,xx} = doIV[ydata,y1data,x2data],
			
		_,  (* Default is TwoStage  *)
				(* form first stage *)
			temp = doReg[y1data,Join[x1data,x2data]];
			yhat = temp[[3]].Join[x1data,x2data];
			
				(* Adjust constant to first row if present *)
			If[Con,
				PrependTo[yhat,x1data[[1]]];
				PrependTo[y1data,x1data[[1]]];
				x1data = Drop[x1data,1];
			];
				(* Form estimator -- or second stage  *)
			temp = doReg[ydata,Join[yhat,x1data]];
			beta = temp[[3]];
			xx = temp[[4]]
	];
	

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

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

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

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

	n = Min[n,ny1,nx1,nx2];
	
	If[y1Con,PrependTo[y1data,Table[1.0, {n}]] ];
	If[x1Con,PrependTo[x1data,Table[1.0, {n}]] ];
	If[x2Con,PrependTo[x2data,Table[1.0, {n}]] ];
	
	If[TrueQ[!(temp = missingObservations[Join[ydata,
		y1data,x1data,x2data],n,k1+k2])],
			Message[MessageName[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}];
	
	Return[{ydata,y1data,x1data,x2data,n,g,g1,k1,k2,Con,
		retainedCases, missObs}]
]

doIV[y_,z_,w_] :=
Block[{delta,xx,zwInv},
	zwInv = Inverse[z.Transpose[w]];
	xx = Transpose[zwInv].
			(w.Transpose[w]).zwInv;
	delta = (y.Transpose[w]).zwInv;
	Return[{delta,xx}]
]

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."

(* Utility Programs for Reg *)

Mean[list_List] := N[Apply[Plus,list]/Length[list]]

Correlation[x1_List,x2_List] :=
Block[{x1bar,x2bar},
	x1bar = Mean[x1];
	x2bar = Mean[x2];
	(x1-x1bar).(x2-x2bar)/Sqrt[(x1-x1bar).(x1-x1bar)*
									(x2-x2bar).(x2-x2bar)]
]/; Length[x1] == Length[x2]

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

printString[str_,args__] := WriteString["stdout",StringForm[
	str, args] ]

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

reduceRank[x_] := If[MatrixQ[x]&&(Min[Dimensions[x]]==1)||
	VectorQ[x]&&(Length[x]==1),
	Return[x[[1]]],Return[x] ]
	
(* lag[list_List,n_:1,m_:undefined] :=
Block[{i,temp},
	If[TrueQ[m==undefined],
		Which[
			n<0,
				temp = Nest[Drop[#,1]&,list,-n],
			n==0,
				temp = list,
			n>0,
				temp = Drop[Nest[Prepend[#,"a"]&,list,n],-n]
		],
		temp = Table[lag[list,i], {i,n,m}]
	];
	temp = temp  (* This kludge corrects for Mathematica bug
				when Return[] is used in a function with
				optional parameters with defaults.  Remove
				temps and replace with Returns when bug
				is fixed. *)
]/;(m==undefined || m>n)  old version  *)

lag[list_List,n_:1,m_:undefined] :=
Block[{i,temp},
	If[TrueQ[m==undefined],
		Which[
			n<0,
				temp = Join[Drop[list,-n],Table["a",{-n}]],
			n==0,
				temp = list,
			n>0,
				temp = Join[Table["a",{n}],Drop[list,-n]]
		],
		temp = Table[lag[list,i], {i,n,m}]
	];
	Return[temp]
]/;(m==undefined || m>n)

End[]


EndPackage[]
