A comprehensive discussion of the stochastic master equation approach to chemical reactions is presented. After a review of the various system size expansion methods that have been developed for an asymptotic solution of stochastic master equations, a new technique, the Poisson representation method is introduced. The technique is based on an expansion of the probability distribution in Poisson distributions and enables chemical master equations to be transformed into exact Fokker-Planck equations. The Fokker-Planck equations obtained using the Poisson representation are used to investigate some general features of chemical reaction systems whose steady state is thermodynamic equilibrium. The solutions for such systems are discussed from grand canonical and canonical ensemble points of view. Basic formulae for the two time correlation functions are derived using the Poisson representation method. The formulae for the chemical systems in thermodynamic equilibrium are shown to relate directly to the fluctuation-dissipation theorems which may be derived from equilibrium statistical mechanical considerations. For non-equilibrium systems, the formulae are shown to be generalisations of these fluctuation-dissipation theorems, but containing an extra term which arises entirely from the non equilibrium nature of the system. In a chemical reaction mechanism, if only bimolecular reactions occur, the Fokker-Planck equations derived using the Poisson representation method involve no higher than second order derivatives. The mathematical equivalence between Fokker-Planck equations and stochastic differential equations is used to derive stochastic differential equations for chemical reaction systems. When trimolecular or higher reactions are present, the Fokker-Planck equations involve third or higher order derivatives. It is shown how the concept of a stochastic differential equation may be generalised in such cases. The stochastic differential equations lend themselves to a straightforward perturbative expansion in the inverse powers of the system size. These perturbative methods are applied to some interesting examples of non equilibrium reaction systems. Some novel features of the Poisson representation quasi probability which arise in these cases are discussed. A chemical reaction model which exhibits a second order phase transition behaviour is investigated in detail (including the process of spatial diffusion) and the first two terms in a perturbative expansion for the spatial and two time correlation functions are calculated. A stochastic model of a chemical reaction which exhibits spatial and temporal oscillations is analysed below the instability thresholds. The same model is used to show how the Poisson representation method may be applied to calculate the spatial and two time correlation functions in finite systems with boundary conditions.