We investigate the sub-crustal stress in the Taiwan region. A tectonic configuration in this region is dominated by a collision between the Philippine oceanic plate and the Eurasian continental margin. The horizontal components of the sub-crustal stress are computed based on the modified Runcorn’s formulae in terms of the stress function with a subsequent numerical differentiation. This modification increases the (degree-dependent) convergence domain of the asymptotically-convergent series and consequently allows evaluating the stress components to a spectral resolution, which is compatible with currently available global crustal models. Moreover, the solution to the Vening Meinesz-Moritz’s (VMM) inverse isostasy problem is explicitly incorporated in the stress function definition. The sub-crustal stress is then computed for a variable Moho geometry, instead of assuming only a constant Moho depth. The regional results reveal that the Philippine plate subduction underneath the Eurasian continental margin generates the shear sub-crustal stress along the Ryukyu Trench. Some stress anomalies associated with this subduction are also detected along both sides of the Okinawa Trough. A tensional stress along this divergent tectonic plate boundary is attributed to a back-arc rifting. The sub-crustal stress, which is generated by a (reverse) subduction of the Eurasian plate under the Philippine plate, propagates along both sides of the Luzon (volcanic) Arc. This stress field has a prevailing compressional pattern.