The Milestoning method has achieved great success in the calculation of equilibrium kinetic properties such as rate constants from molecular dynamics simulations. The goal of this work is to advance Milestoning into the realm of non-equilibrium statistical mechanics, in particular, the calculation of time correlation functions. In order to accomplish this, we introduce a novel methodology for obtaining flux through a given milestone configuration as a function of both time and initial configuration, and build upon it with a novel formalism describing autocorrelation for Brownian motion in a discrete configuration space. The method is then applied to three different test systems: a harmonic oscillator, which we solve analytically, a two well potential, which is solved numerically, and an atomistic molecular dynamics simulation of alanine dipeptide.