A Monte Carlo algorithm has been implemented which couples the results of a three-dimensional magnetohydrodynamic equilibrium code with particle transport calculations similar to those of Boozer and Kuo-Petravic. The equilibrium results are in a form which lends itself to easy and accurate computation of the coefficients required to follow particle orbits in flux variables at finite pressure. In particular. the parallel current and the Clebsch potential for the current are obtained from Fourier series solutions of first-order partial differential equations with constant coefficients.Confinement times are estimated from the exponential decay of expected values of appropriately chosen functionals of the particle distribution. The observation that these functionals satisfy boundary conditions helps us to compute confinement times over a wide range of collision frequencies. including cases where losses due to particle trapping are very high. Initial conditions are chosen to optimize the Monte Carlo calculation by using known information about the expected particle distribution.Results for actual and proposed stellarator experiments are given. Electron and ion confinement times are compared, since this is relevant to the issue of the effect of ambipolar electric fields on stellarator confinement. A spectral method is given to determine the elcctric field from charge neutrality. and numerical evidence is presented to suggest that anomalous electron transport may be due to small resonant terms in the electric potential. Also, comparisons are made with a simplified theory of particle transport already incorporated in the equilibrium code. zero net current stellarators are characterized by a quiet plasma, which might be more accurately represented by the neoclassical transport theory than tokamaks, which exhibit turbulence and anomalous transport; see [13].We take full account of the complicated geometry and finite p effects by solving the equilibrium problem numerically with a computer code BETA that has been described in detail elsewhere; see [l]. This results in a solution of the MHD equations where B is the magnetic field, J is the current density and p is the fluid pressure. A basic assumption is that there exists a nested set of flux surfaces s = const., with the pressure restricted to be a function p = p ( s ) of s alone. The magnetic field can be represented in terms of flux functions s and JI and Clebsch potentials 9 and S. We have B = 0 s X vJI = V+ + ~J v s , and as a consequence J = p,VS X VS.We introduce s as an independent variable in the radial direction, together with two angle variables u and u. These variables are normalized so that u has a unit period in the poloidal direction while remaining periodic in the toroidal direction, whereas u is periodic in the poloidal direction with unit period in the toroidal direction. The functions JI, + and 5 are multi-valued on each torus s = const., with periods given by the representations are single-valued. Here the total toroidal flux is equal to s, so that the pe...