(*: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[k20, 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[]