Content

* ITSP 14.8 -

mata: mata clear

/*
// generate residual series
local minn = .
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'
                predict double `eps`i'' if e(sample), 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
*/

do mm_suregub.mata

type mm_suregub.mata
/*
version 10.1
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 = xx = yy = J(1, neq, NULL)
        isigma = invsym(st_matrix(ssigma))        
        nrow = 0
        ncol = 0
        string rowvector coefname, eqn
        string matrix mstripe
// equation loop 1
        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
// equation loop 2
        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
*/

* webuse grunfeld, clear
do mm_suregub.mata

use grunfeld, clear
drop in 75/80
drop in 41/43
drop in 18/20
keep if company <= 4
drop time
reshape wide invest mvalue kstock, i(year) j(company)
list year invest*, sep(0)

suregub (invest1 mvalue1 kstock1) (invest2 mvalue2 kstock2) ///
        (invest3 mvalue3 kstock3) (invest4 mvalue4 kstock4)