This paper examines upper and lower bounds of the effective elastic modulus of unidirectional short-fiber composites. The short-fibers are modeled by aligned ellipsoidal inclusions of the same aspect ratio but not necessarily the same size. We adopt a perturbation expansion of the composite local strain field by using the Green function tensor. Explicit expressions of the effective elastic modulus are derived up to the third-order term by use of the information on the correlation functions. The variational method is then employed to optimize the bounds of the effective modulus in a closed form. Numerical examples of the bounds as functions of the fiber aspect ratio and the fiber volume fraction are given for a glass/epoxy system. The present approach predicts narrower bounds than those of Hashin and coworkers for the limiting cases of spherical particles and continuous fibers since their bounds corresponds to a model that take the correlation functions up to the second order into account.