In the framework of JAXA's MMX mission to explore the Martian moon of Phobos, an analysis of the stability of three-dimensional quasi-satellite orbits in the Mars-Phobos circular restricted three-body problem was conducted. For this analysis, notions of co-orbital motion, interpreted as the slow motion of the guiding center of the trajectory along the disturbing potential of Phobos, were used. After identifying and analyzing different regions of stability for three quasi-satellite orbits at 100, 50 and 30 km from the center of Phobos, several conclusions were drawn regarding the dynamics of the ballistic escape of the spacecraft, interpreted in terms of co-orbital motion. By making use of these insights, a novel methodology to find periodic quasi-satellite orbits able to reach high latitudes over the surface of Phobos was derived. This methodology consists of two steps: a multi-objective minimization using co-orbital parameters as target functions, to isolate regions with potential periodic orbits; followed by a shooting algorithm to arrive at the final periodic orbit. As a result of this new methodology, two periodic orbits were found at 50 and 30 km from the center of Phobos, able to reach latitudes as high as 54 deg and 32 deg respectively. This result represents an important contribution to both the operations and the scientific return of the Phobos proximity phase within MMX.