Energy amplification in channel flows of Oldroyd-B fluids is studied from an input-output point of view by analyzing the ensemble-average energy density associated with the velocity field of the linearized governing equations. The inputs consist of spatially distributed and temporally varying body forces that are harmonic in the streamwise and spanwise directions and stochastic in the wall-normal direction and in time. Such inputs enable the use of powerful tools from linear systems theory that have recently been applied to analyze Newtonian fluid flows. It is found that the energy density increases with a decrease in viscosity ratio (ratio of solvent viscosity to total viscosity) and increase in Reynolds number and elasticity number. In most of the cases, streamwise constant perturbations are most amplified and the location of maximum energy density shifts to higher spanwise wavenumbers with an increase in Reynolds number and elasticity number and decrease in viscosity ratio. For similar parameter values, the maximum in the energy density occurs at a higher spanwise wavenumber for Poiseuille flow, whereas the maximum energy density is higher for Couette flow. At low Reynolds numbers, the energy density decreases monotonically when the elasticity number is sufficiently small, but shows a maximum when the elasticity number becomes sufficiently large, suggesting that elasticity can amplify disturbances even when inertial effects are weak.