We study the average of the product of the central values of two $L$-functions of modular forms $f$ and $g$ twisted by Dirichlet characters to a large prime modulus $q$. As our principal tools, we use spectral theory to develop bounds on averages of shifted convolution sums with differences ranging over multiples of $q$, and we use the theory of Deligne and Katz to estimate certain complete exponential sums in several variables and prove new bounds on bilinear forms in Kloosterman sums with power savings when both variables are near the square root of $q$. When at least one of the forms $f$ and $g$ is non-cuspidal, we obtain an asymptotic formula for the mixed second moment of twisted $L$-functions with a power saving error term. In particular, when both are non-cuspidal, this gives a significant improvement on M.~Youngs asymptotic evaluation of the fourth moment of Dirichlet $L$-functions. In the general case, the asymptotic formula with a power saving is proved under a conjectural estimate for certain bilinear forms in Kloosterman sums.