The expectation values of $\hat{L}_{k}$ for $k \in \{x, y, z\}$ are indeed always real, since they are hermitian operators.
The commutator of $\hat{L}_{j}$ and $\hat{L}_{k}$ for $j \neq k$ on the other hand is not hermitian, but anti hermitian: \begin{align} [\hat{L}_{j}, \hat{L}_{k}]^{\dagger} &= (\hat{L}_{j}\hat{L}_{k})^{\dagger} - (\hat{L}_{k}\hat{L}_{j})^{\dagger}\\ &= \hat{L}_{k}^{\dagger}\hat{L}_{j}^{\dagger} - \hat{L}_{j}^{\dagger}\hat{L}_{k}^{\dagger} \\ &= \hat{L}_{k}\hat{L}_{j} - \hat{L}_{j}\hat{L}_{k} \\ &= [\hat{L}_{k},\hat{L}_{j}] \\ &= -[\hat{L}_{j},\hat{L}_{k}] \,. \end{align} Another way of seeing that is directly is by using the angular momentum algebra: \begin{align} [\hat{L}_{j}, \hat{L}_{k}] = i \sum_{l = 1}^{3} \epsilon_{jkl} \hat{L}_{l} \end{align}\begin{align} [\hat{L}_{j}, \hat{L}_{k}] = i \hbar \sum_{l = 1}^{3} \epsilon_{jkl} \hat{L}_{l} \end{align} If you take the hermitian conjugate of this equation you get an extra minus on the right hand side due to the implicit complex conjugation of the imaginary unit, while $\hat{L}_{l}$ is mapped to itself.
The eigenvalues of anti hermitian operators must be purely imaginary. Therefore, your result is perfectly valid.
I would suggest to try calculating the matrix elements of $\hat{L}_{x}$ and $\hat{L}_{y}$ by first evaluating all matrix elements of the Ladder operators \begin{align} \hat{L}_{\pm} &= \hat{L}_{x} \pm i \hat{L}_y \end{align} which fulfill the relations [Source: https://en.wikipedia.org/wiki/Ladder_operator#Angular_momentum]: \begin{align} \hat{L}_{\pm}Y_{l,m} &= \hbar \sqrt{l(l+1) - m(m \pm 1)} Y_{l,m \pm 1}\,. \end{align} I hope this was helpful to you!