This study simulates strong ground motions in the Tokyo metropolitan area during the 1923 Kanto earthquake using the stochastic Green’s function method. Source characteristics were modeled using seismic intensity inversion analysis, and path and site characteristics were modeled using inhomogeneous attenuation structure and empirical amplification factors. The results of these simulations were consistent with seismic intensities estimated based on the collapse rate of wooden houses. The distribution of pseudovelocity response spectra averaged at periods of 1–2 s was large: ∼200 cm/s in southern Kanagawa and southern Chiba prefectures, ∼100–200 cm/s in eastern Tokyo, and ∼50–100 cm/s in eastern Saitama prefecture despite its distance from strong‐motion generation areas (SMGAs). The simulation results were regressed on site characteristics and fault distance, and the residuals were interpolated using the Kriging method to estimate detailed maps of seismic intensity and response spectra on an ∼250 m mesh reflecting site‐specific characteristics. The following conclusions can be made: (1) all SMGAs, other than those in northern Tokyo Bay, were located near large slip areas based on coseismic geodetic and seismic waveform data. Although the SMGAs in the northern part of Tokyo Bay exerted little influence on the southern part of the Kanto region, their consideration was required to reproduce the seismic intensity at the northwest Tokyo and Saitama; (2) intense strong motion in central Tokyo occurred at the back marsh, delta, coastal lowlands, and filled lands, whereas low levels of strong motion were determined at terraces covered with volcanic ash soil. Combined with building distribution, this indicates areas of high seismic risk; (3) the seismic intensity and response spectra in the Tokyo metropolitan area obtained through this simulation were larger than those obtained from seismic records of the 2011 Tohoku earthquake—the most recent megathrust earthquake.