Simulations in the warm dense matter regime using finite temperature Kohn-Sham density functional theory (FT-KS-DFT), while frequently used, are computationally expensive due to the partial occupation of a very large number of high-energy KS eigenstates which are obtained from subspace diagonalization. We have developed a stochastic method for applying FT-KS-DFT, that overcomes the bottleneck of calculating the occupied KS orbitals by directly obtaining the density from the KS Hamiltonian. The proposed algorithm, scales as $O\left(NT^{-1}\right)$ and is compared with the high-temperature limit scaling $O\left(N^{3}T^{3}\right)$ of the deterministic approach, where $N$ is the system size (number of electrons, volume etc.) and $T$ is the temperature. The method has been implemented in a plane-waves code within the local density approximation (LDA); we demonstrate its efficiency, statistical errors and bias in the estimation of the free energy per electron for a diamond structure silicon. The bias is small compared to the fluctuations, and is independent of system size. In addition to calculating the free energy itself, one can also use the method to calculate its derivatives and obtain the equations of state.