In this talk I consider a single-server queue with Levy input, and in particular its workload process $Q(t)$, focusing on its correlation structure. With the correlation function defined as $r(t) := {\\rm Cov}(Q(0), Q(t))/{\\rm Var} Q(0)$ (assuming the workload process is in stationarity at time 0), we first study its transform $\\int_0^\\infty r(t)e^{-\\theta t} dt,$ both for the case that the Levy process has positive jumps, and that it has negative jumps. These expressions allow us to prove that $r(t)$ is positive, decreasing, and convex, relying on the machinery of completely monotone functions. For the light-tailed case, we estimate the behavior of $r(t)$ for $t$ large. We then focus on techniques to estimate $r(t)$ by simulation. Naive simulation techniques require roughly $1/r(t)^2$ runs to obtain an estimate of a given precision, but we develop a coupling technique that leads to substantial variance reduction (required number of runs being roughly $1/r(t)$). If this is augmented with importance sampling, it even leads to a logarithmically efficient algorithm.