We present general arguments and construct a stress tensor operator for finite lattice spin models. The average value of this operator gives the Casimir force of the system close to the bulk critical temperature $T_c$. We verify our arguments via exact results for the force in the two-dimensional Ising model, $d$-dimensional Gaussian and mean spherical model with $2<d<4$. On the basis of these exact results and by Monte Carlo simulations for three-dimensional Ising, XY and Heisenberg models we demonstrate that the standard deviation of the Casimir force $F_C$ in a slab geometry confining a critical substance in-between is $k_b T D(T)(A/a^{d-1})^{1/2}$, where $A$ is the surface area of the plates, $a$ is the lattice spacing and $D(T)$ is a slowly varying nonuniversal function of the temperature $T$. The numerical calculations demonstrate that at the critical temperature $T_c$ the force possesses a Gaussian distribution centered at the mean value of the force $<F_C>=k_b T_c (d-1)Delta/(L/a)^{d}$, where $L$ is the distance between the plates and $Delta$ is the (universal) Casimir amplitude.