program drop _all
mata: mata clear
*! suregub CFBaum 1.0.0 23oct2007
* 1.0.0: initial release
program suregub, eclass
version 9.2
// from reg3.ado 6.4.1
// not supported v9 local cmdline : copy local 0
local cmdline `0'
// disable options (such as constraints)
if strpos("`0'",",") {
di as err _n "Options not permitted"
exit 198
}
/* Parse (y1 x1 x2) (y2 y1 x2 x3) structure.
* Pick up full varlist (flist), y-array (y`i'), independent
* variables (ind`i') left-hand sides (lhslist),
* equation names (eqnm`i')
*/
local neq = 0
/* Parse the equations ([eqname:] y1 [y2 y3 =] x1 x2 x3)
* and fill in the structures required for estimation */
gettoken fulleq 0 : 0, parse(" ,[") match(paren)
IsStop `fulleq'
while `s(stop)' == 0 {
if "`paren'" != "(" { /* must be eq */
eq ?? `fulleq'
local fulleq "`r(eqname)': `r(eq)'"
}
parseEqn `fulleq'
/* Set-up equation bookeeping structures */
local flist `flist' `depvars' `indvars'
tokenize `depvars'
local i 1
while "``i''" != "" {
local neq = `neq' + 1
local y`neq' ``i''
local lhslist `lhslist' ``i''
local ind`neq' `indvars'
local cons`neq' = ("`constan'" == "")
nameEq "`eqname'" "``i''" "`eqlist'" `neq'
local eqnm`neq' = "`s(eqname)'"
local eqlist `eqlist' `s(eqname)'
local i = `i' + 1
}
gettoken fulleq 0 : 0, parse(" ,[") match(paren)
IsStop `fulleq'
}
local 0 `"`fulleq' `0'"'
if `neq' < 1 {
di in red "equation(s) required"
exit 198
}
// generate residual series
local minn = 999999
local maxn = 0
forv i=1/`neq' {
local dv : word `i' of `eqlist'
local eq`i' = "`dv' `ind`i''"
qui {
reg `dv' `ind`i''
tempvar touse`i' es eps`i'
g byte `es' = e(sample)
predict double `eps`i'' if `es', resid
g byte `touse`i'' = cond(e(sample),1,.)
su `eps`i'', meanonly
local maxn = max(`maxn',r(N))
local minn = min(`minn',r(N))
}
}
tempname sigma
matrix `sigma' = J(`neq',`neq',0)
// generate pairwise correlation matrix of resids;
// for comparison with sureg, use divisor N
local neq1 = `neq'-1
forv i = 1/`neq1' {
forv j = 2/`neq' {
qui correlate `eps`i'' `eps`j'', cov
mat `sigma'[`i',`i'] = r(Var_1)*(r(N)-1)/(r(N))
mat `sigma'[`j',`j'] = r(Var_2)*(r(N)-1)/(r(N))
mat `sigma'[`i',`j'] = r(cov_12)*(r(N)-1)/(r(N))
mat `sigma'[`j',`i'] = `sigma'[`i',`j']
}
}
mata: mm_suregub(`neq',"`eqlist'","`sigma'")
di _n "Seemingly unrelated regressions for an unbalanced panel"
di _n "Min obs per unit = `minn'"
di "Max obs per unit = `maxn'"
mat b = r(b)
mat V = r(V)
eret clear
// mat list b
// mat list V
eret post b V
eret local cmd "suregub"
eret local minobs `minn'
eret local maxobs `maxn'
eret display
end
version 9.2
mata:
void mm_suregub(real scalar neq, string scalar eqlist, string scalar ssigma)
{
pointer (real matrix) rowvector eq
pointer (real matrix) rowvector xx
pointer (real matrix) rowvector yy
eq = J(1,neq,NULL)
xx = J(1,neq,NULL)
yy = J(1,neq,NULL)
isigma = invsym(st_matrix(ssigma))
// isigma
nrow = 0
ncol = 0
string rowvector coefname, eqn
string matrix mstripe
for(i=1;i<=neq;i++) {
lt = "touse"+strofreal(i)
touse = st_local(lt)
st_view(tt,.,touse)
le = "eq"+strofreal(i)
eqv = st_local(le)
vars=tokens(eqv)
v = vars[|1,.|]
// pull in full matrix, including missing values
st_view(eqq,.,v)
eq[i] = &(tt :* eqq)
// matrix eq[i] is [y|X] for ith eqn
eqname=v[1]
stripe = v[2::cols(v)], "_cons"
coefname = coefname, stripe
eqn = eqn, J(1,cols(v),eqname)
// form X, assuming constant term
nrow = nrow + rows(*eq[i])
iota = J(rows(*eq[i]),1,1)
xx[i] = &((*eq[i])[| 1,2 \ .,. |],iota)
ncol = ncol + cols(*xx[i])
// form y
yy[i] = &(*eq[i])[.,1]
}
XX = J(ncol,ncol,0)
YY = J(ncol,1,0)
ii = 0
for(i=1;i<=neq;i++) {
i2 = cols(*xx[i])
xi = *xx[i]
jj = 0
for(j=1;j<=neq;j++) {
xj = *xx[j]
j2 = cols(*xx[j])
yj = *yy[j]
XX[| ii+1,jj+1 \ ii+i2,jj+j2 |] = isigma[i,j] :* cross(xi,xj)
YY[| ii+1, 1 \ ii+i2, 1 |] = YY[| ii+1, 1 \ ii+i2, 1 |] + isigma[i,j] :* cross(xi,yj)
jj = jj + j2
}
ii = ii + i2
}
// compute SUR beta (X' [Sigma^-1 # I] X)^-1 (X' [Sigma^-1 # I] y)
vee = invsym(XX)
beta = vee*YY
st_matrix("r(b)",beta')
mstripe=eqn',coefname'
st_matrixcolstripe("r(b)",mstripe)
st_matrix("r(V)",vee)
st_matrixrowstripe("r(V)",mstripe)
st_matrixcolstripe("r(V)",mstripe)
}
end
program define parseEqn
version 6
/* see if we have an equation name */
gettoken token uu : 0, parse(" =:") /* rare, pull twice if found */
gettoken token2 : uu, parse(" =:") /* rare, pull twice if found */
if index("`token2'", ":") != 0 {
gettoken token 0 : 0, parse(" =:") /* sic, to set 0 */
gettoken token2 0 : 0, parse(" =:") /* sic, to set 0 */
c_local eqname `token'
}
else c_local eqname
/* search just for "=" */
gettoken token 0 : 0, parse(" =")
while "`token'" != "=" & "`token'" != "" {
local depvars `depvars' `token'
gettoken token 0 : 0, parse(" =")
}
if "`token'" == "=" {
tsunab depvars : `depvars'
syntax [varlist(ts)] [ , noConstant ]
}
else { /* assume single depvar */
local 0 `depvars'
syntax varlist(ts) [ , noConstant ]
gettoken depvars varlist : varlist
}
c_local depvars `depvars'
c_local indvars `varlist'
c_local constan `constan'
end
program define IsStop, sclass
version 6
if `"`0'"' == "[" /*
*/ | `"`0'"' == "," /*
*/ | `"`0'"' == "if" /*
*/ | `"`0'"' == "in" /*
*/ | `"`0'"' == "" {
sret local stop 1
}
else sret local stop 0
end
/* Drop all duplicate tokens from list */
program define DropDup /* <newlist> : <list> */
version 6
args newlist /* name of macro to store new list
*/ colon /* ":"
*/ list /* list with possible duplicates */
gettoken token list : list
while "`token'" != "" {
local fixlist `fixlist' `token'
local list : subinstr local list "`token'" "", word all
gettoken token list : list
}
c_local `newlist' `fixlist'
end
/* Remove all tokens in dirt from full */
* Returns "cleaned" full list in cleaned */
program define Subtract /* <cleaned> : <full> <dirt> */
version 6
args cleaned /* macro name to hold cleaned list
*/ colon /* ":"
*/ full /* list to be cleaned
*/ dirt /* tokens to be cleaned from full */
tokenize `dirt'
local i 1
while "``i''" != "" {
local full : subinstr local full "``i''" "", word all
local i = `i' + 1
}
c_local `cleaned' `full' /* cleans up extra spaces */
end
/* Returns tokens found in both lists in the macro named by matches.
* Duplicates must be duplicated in both lists to be considered
* matches a 2nd, 3rd, ... time. */
program define Matches
version 6
args matches /* macro name to hold cleaned list
*/ colon /* ":"
*/ list1 /* a list of tokens
*/ list2 /* a second list of tokens */
tokenize `list1'
local i 1
while "``i''" != "" {
local list2 : subinstr local list2 "``i''" "", /*
*/ word count(local count)
if `count' > 0 {
local matlist `matlist' ``i''
}
local i = `i' + 1
}
c_local `matches' `matlist'
end
/* Find all occurances in List of tokens in FindList and replace with
* corresponding token from SubsList. Assumes FindList and SubstList
* have same number of elements.
*/
program define Subst, sclass /* <NewList> : <List> <FindList> <SubstList> */
version 6
args newname /* macro name to hold list after replacements
*/ colon /* ":"
*/ list /* varlist with tokens to be replaced
*/ fndList /* list of tokens to be replaced
*/ subList /* varlist with replacement tokens */
tokenize `fndList'
local i 1
while "``i''" != "" {
gettoken repltok subList : subList
local list : subinstr local list "``i''" "`repltok'", word all
local i = `i' + 1
}
c_local `newname' `list'
end
/* determine equation name */
program define nameEq, sclass
version 6
args eqname /* user specified equation name
*/ depvar /* dependent variable name
*/ eqlist /* list of current equation names
*/ neq /* equation number */
if "`eqname'" != "" {
if index("`eqname'", ".") {
di in red "may not use periods (.) in equation names: `eqname'"
}
local eqlist : subinstr local eqlist "`eqname'" "`eqname'", /*
*/ word count(local count) /* overkill, but fast */
if `count' > 0 {
di in red "may not specify duplicate equation names: `eqname'"
exit 198
}
sreturn local eqname `eqname'
exit
}
local depvar : subinstr local depvar "." "_", all
if length("`depvar'") > 32 {
local depvar "eq`neq'"
}
Matches dupnam : "`eqlist'" "`depvar'"
if "`dupnam'" != "" {
sreturn local eqname = substr("`neq'`depvar'", 1, 32)
}
else {
sreturn local eqname `depvar'
}
end