Arbitrarily high order discontinuous Galerkin method in time combined with discontinuous Galerkin method with Lagrange multiplier (DGLM) in space is proposed to approximate the solution to hyperbolic conservation laws with boundary conditions. Stability of the approximate solution is proved in a broken $L^2(L^2)$ norm and also in an $l^\infty(L^2)$ norm. Error estimates of ${\mathcal{O}}(h^{r+\frac12}+k_n^{q+\frac12})$ with $P_r(E)$ and $P_q(J_n)$ elements $(r, q\ge \frac{d+1}2)$ are derived in a broken $L^2(L^2)$ norm, where $h$ and $k_n$ are the maximum diameters of the elements and the time step of $J_n$, respectively, $J_n$ is the time interval, and $d$ is the dimension of the spatial domain. An explanation on algorithmic aspects is given. $P_0$ time and space subcell limiting processes are applied to resolve the shocks. It is numerically shown that the high order DG-DGLM method is well-suited for long time integrations. Several numerical experiments for advection, shallow water, and compressible Euler equations are presented to show the performance of the high order DG-DGLM with $P_0$ time and space subcell limiting processes.