Population genomic datasets collected over the past decade have spurred interest in developing methods that can utilize massive numbers of loci for inference of demographic and selective histories of populations. The site frequency spectrum (SFS) provides a convenient framework for such analysis and accordingly much attention has been paid to predicting theoretical expectations of the SFS under a number of different models. However, to date, exact solutions for the joint SFS of two or more population under models of migration and divergence have not been found. Here we present a novel Markov chain representation of the coalescent on the state space of the joint SFS that allows for rapid, exact calculation of the joint SFS under generalized isolation with migration (IM) models. In turn, we show how our Markov chain method, in the context of composite likelihood estimation, can be used for accurate inference of parameters of the IM model using SNP data. Lastly, we apply our method to recent whole genome datasets from Drosophila melanogaster and human populations.