The stresses and the elastic constants of bcc sodium are calculated by molecular dynamics (MD) for temperatures to T = 340 K. The total adiabatic potential of a system of sodium atoms is represented by a pseudopotential model. The resulting expression has two terms: a large, strictly volume-dependent potential, plus a sum over ion pairs of a small, volume-dependent two-body potential. The stresses and the elastic constants are given as strain derivatives of the Helmholtz free energy. The resulting expressions involve canonical ensemble averages (and fluctuation averages) of the position and volume derivatives of the potential. An ensemble correction relates the results to MD equilibrium averages. Evaluation of the potential and its derivatives requires the calculation of integrals with infinite upper limits of integration, and integrand singularities. Methods for calculating these integrals and estimating the effects of integration errors are developed. A method is given for choosing initial conditions that relax quickly to a desired equilibrium state. Statistical methods developed earlier for MD data are extended to evaluate uncertainties in fluctuation averages, and to test for symmetry. The fluctuation averages make a large contribution to the elastic constants, and the uncertainties in these averages are the dominant uncertainties in the elastic constants. The strictly volume-dependent terms are very large. The ensemble correction is small but significant at higher temperatures. Surprisingly, the volume derivatives of the two-body potential make large contributions to the stresses and the elastic constants. The effects of finite potential range and finite system size are discussed, as well as the effects of quantum corrections and electronic excitations. The agreement of theory and experiment is very good for the magnitudes of C₁₁ and C₁₂. The magnitude of C₄₄ is consistently small by ∼9 kbar for finite temperatures. This discrepancy is most likely due to the neglect of three-body contributions to the potential. The agreement of theory and experiment is excellent for the temperature dependences of all three elastic constants. This result illustrates a definite advantage of MD compared to lattice dynamics for conditions where classical statistics are valid. MD methods involve direct calculations of anharmonic effects; no perturbation treatment is necessary.