The matrix diagonal represents the self-potential. Normally, when you solve for Moran's-I you remove the diagonal of the matrix. However, I cannot find anything in the GeoDa documentation that explains the default behavior of the statistic. I have never seen options for deriving interzonal weights so, I would imagine that, since it is 0, the diagonal is removed. You may need to contact the authors for a definitive answer.
I know that the ArcGIS implementation has an option for including self-potential in the univariate case (sorry, no bivariate implementation). However, they do not have a correct test statistic and the z-value and p-values are not stable when faced with non-normal distributions.
Likely the only way you will be able to implement this is to code it yourself in something like PySal or R. The univariate interzonal weight is calculated as: dij = 0.5 * [(Aij / π)**0.5] but you will have to figure out the bivariate adjustment. If you implement this you will have to think long and hard about the meaning. I am not sure that it is going to provide you with what your are thinking.
You could consider using a scan statistic that is better suited for spatial time-series analysis under specific distributional assumptions. This would provide you with a more suitable framework for hypothesis testing. I would also look into Crimestat. As I recall there is some flexibility in defining the behavior of contingency and there is a bivariate Moran's-I/LISA.
You may also consider an autoregressive model (Li et al., 2007). This method re-scales the measure by a function of the eigenvalues of the spatial weights matrix and provides a much more robust measure of the spatial dependence.
Li, H., C.A. Calder and N. Cressie. (2007). Beyond Moran’s I: Testing for spatial
dependence based on the spatial autoregressive model. Geographical Analysis
39:357–375.
Sorry for not being able to provide a more definitive answer.
Getting a p value is not that easy. Marcon and Lang (Testing randomness of spatial point patterns with the Ripley statistic, http://arxiv.org/abs/1006.1567) have demonstrated that (K1, K2..Kn) where K1, K2..Kn are the values of the Ripley s K function at different distances is a Gaussian vector. They have then computed its mean and covarianace matrix in the square and used a chi 2 test for spatial randomness. You can use the same idea by computing empirically (with Monte-Carlo simulations) the mean/covariance matrix in your domain.
Regards
Thibault Lagache
Best Answer
After much searching in the back corners of ESRI documentation, I've concluded that there is no reasonable way run a bivariate Ripley's K function in Arcpy/ArcGIS. However, I have found a solution using R: