We present a numerical model of the subducting lithosphere that provides an alternative explanation for stresses causing deep earthquakes. Our model lithosphere is composed of α olivine, ¿ spinel, γ spinel, and perovskite + magnesiow¿stite. The heat conduction equation is solved to determine temperature conditions in the slab and to locate the equilibrium phase transformations in pressure-temperature space. Volumetric strains in the subducting lithosphere are calculated from the density of individual phases and from the heat released or consumed in the phase changes. These strains are used as sources of stress in the subducting lithosphere. Dislocation creep and Peierls stress creep laws are included in the viscoelastic rheology. Volumetric reductions due to equilibrium phase transformations cause high shear stress in the transition zone because of the variable viscosity inside the subducting slab. Aspects of the model shear stresses are in agreement with observations of high seismic activity in the Tonga-Wadati Benioff zone. Compression is oriented along the dip of the slab, and extension is oriented in the plane perpendicular to the compression axis. Since our model stresses agree with the seismic observations, and because the model stresses are larger than those caused by buoyancy forces, our model provides a possible explanation of stresses causing deep earthquakes. Also, our model does not need metastability of olivine to explain the occurrence of high shear stress in the transition zone.