A higher order time differencing method for the spatially nonhomogeneous Boltzmann equation is derived from the integral form of the equation along its characteristic line. Similar to the splitting method, which solves the collisionless equation in the convection step and the spatially homogeneous Boltzmann equation in the collision step, the present method consists of two steps, one of which is the same as the convection step in the splitting method. The difference from the splitting method is in the other step, where not only the collision term but also its variation along the characteristic line is taken into account correctly. The truncation error of this method per time step t is O( t 3 ) and its higher order accuracy is demonstrated numerically in the shock propagation problem using the BGK model equation. It is shown that such accuracy is never realized in the framework of the conventional splitting formulation, which is contrary to Bogomolov's result (U.S.S.R. Comput. Math. Math. Phys. 28, 79 (1988)). The other higher order methods based on the integral form are also presented. Furthermore, the extension to the stochastic approach, modification of the conventional direct simulation Monte-Carlo procedure, is proposed.