Modeling longitudinal binary data is challenging but common in practice. Existing methods on modeling of binary responses take no account of the fact that the correlation coefficient of binary responses must have an upper bound which is smaller than one. Ignoring this fact can lead to incorrect statistical inferences for longitudinal binary data. A novel method is proposed to model the mean and within-subject correlation coefficients for longitudinal binary data, simultaneously, by taking into account the constraints of the upper bounds. By introducing latent normally distributed random variables, the correlation coefficients of binary responses are connected to those for the latent variables, of which the correlation coefficients are modeled accordingly. A joint generalized estimating equation (GEE) method is developed for this purpose and the resulting correlation coefficients are shown to satisfy the constraints. Asymptotic normality of the parameter estimators is derived and simulation studies are made under various scenarios, showing that the proposed joint GEE method works very well even if the working covariance structures are misspecified. For illustration, the proposed method is applied to two real data practices to assess the effects of covariates on the mean and within-subject correlation coefficients.