The oscillatory motions of shocks in highly underexpanded jets with nozzle pressure ratios of 5.60, 7.47, 9.34, and 11.21 are quantitatively studied by using large eddy simulation. Two types of shock oscillations are observed: one is the Mach disk oscillation in the streamwise direction and the other is the shock oscillation in the radial direction. It is found that the Mach disk moves quickly in the middle of the oscillatory region but slowly at the top or bottom boundaries. The oscillation cycles of Mach disk are the same for different cases, and are all dominated by an axisymmetric mode of 5.298 kHz. For the oscillation in the radial direction, the shocks oscillate more toward the jet centerline but less in the jet shear layer, and the oscillation magnitude is an increasing function of screech amplitude. The cycles of the radial shock oscillation switch randomly between the two screech frequencies for the first two cases. However, the oscillation periodicity is more complex for the jets with high nozzle pressure ratios of 9.34 and 11.21 than for the jets with the low nozzle pressure ratios of 5.6 and 7.47. In addition, the shock oscillation characteristics are also captured by coarse mesh and Smagorinsky model, but the coarse mesh tends to predict a slower and weaker shock oscillation.