We study randomly shifted rank-1 lattice rules for computing $\mathbb{E}[g(\mathbf{X}_T)]$, where $\mathbf{X}_T$ is the terminal state of a stochastic differential equation (SDE) with constant diffusion matrix. Discretizing the SDE by the Euler–Maruyama scheme with $m$ uniform time steps turns the expectation into an integral whose dimension grows with $m$ and whose integrand has a nested composition structure; both complicate the lattice rule error analysis. We generate the Brownian path by the Brownian bridge construction and handle the nested structure using the multivariate Faà di Bruno formula and a rooted-tree expansion. Under suitable regularity and growth conditions, we prove that an $N$-point randomly shifted lattice rule with a uniformly random generating vector achieves a mean absolute sampling error of $O(N^{-1+\varepsilon})$ for every $\varepsilon\in(0,1/2)$ and every prime $N$, with the implied constant independent of $m$. To improve stability and control higher error moments, the median lattice estimator is used. The estimator takes the median of $R=2k-1$ independent lattice estimates and achieves an $L^q$ error of $O(N^{-1+\varepsilon})$ for every $1\le q<k/(1-\varepsilon)$. We further establish error bounds independent of both the state dimension and the Brownian path dimension under an additional summability condition on the columns of the diffusion matrix. Numerical experiments show that median lattice estimators converge faster than Monte Carlo estimators and remain stable as the number of time steps $m$ grows. For models satisfying the additional summability condition, stability is also observed as the state and Brownian path dimensions increase.