/* Smat - returns the estimated covariance matrix S Pietro Tebaldi , Harvard School of Public Health and Bocconi University July 2010 */ program define Smat , rclass version 10.1 syntax , x(varname) y(varname) g(varname) [bins(integer 20)] tempname d ld D N N1 N2 S _diff , x(`x') y(`y') g(`g') bins(`bins') matrix `d' = r(d) local N = r(N) local N1 = r(N1) local N2 = r(N2) local ld = rowsof(`d') matrix dissimilarity `D' = `x' `y' forvalues k = 1/`ld' { mata: Ind_`k' = J(`N',`N',.) forvalues i = 1/`N' { forvalues j = 1/`N' { if `D'[`i',`j']<=`d'[`k',1] { mata: Ind_`k'[`i',`j'] = 1 } else { mata: Ind_`k'[`i',`j'] = 0 } } } forvalues i = 1/`N' { mata: Ind_`k'[`i',`i']=0 } } mata: S = J(`ld',`ld',.) forvalues s = 1/`ld' { // product matrices for each of the entries of S forvalues k = 1/`ld' { mata: Prod`s'_`k' = Ind_`s'*Ind_`k' mata: prob`s'_`k' = (J(1,`N',1)*Prod`s'_`k'*J(`N',1,1))/(2*comb(`N',2)+(`N'*(`N'-1)*(`N'-2))) mata: S[`s',`k'] = 4*(`N'/`N1'*prob`s'_`k'+`N'/`N2'*prob`s'_`k')/`N' mata: mata drop Prod`s'_`k' mata: mata drop prob`s'_`k' } } forvalues k = 1/`ld' { // clear matrices mata: mata drop Ind_`k' } mata: st_matrix("`S'",S) return matrix S = `S' end