  Celestial operator product expansions (OPEs) arise from the collinear limit
of scattering amplitudes and play a vital role in celestial holography. In this
paper, we derive the celestial OPEs of massless fields in string theory from
the worldsheet. By studying the worldsheet OPEs of vertex operators in
worldsheet CFT and further examining their behaviors in the collinear limit, we
find that new vertex operators for the massless fields in string theory are
generated and become dominant in the collinear limit. Mellin transforming to
the conformal basis yields exactly the celestial OPEs in celestial CFT. We also
derive the celestial OPEs from the collinear factorization of string amplitudes
and the results derived in these two different methods are in perfect agreement
with each other. Our final formulae of celestial OPEs are applicable to general
dimensions, corresponding to Einstein-Yang-Mills theory supplemented by some
possible higher derivative interactions. Specializing to 4D, we reproduce all
the celestial OPEs for gluon and graviton in the literature. We consider
various string theories, including the open and closed bosonic string, as well
as the closed superstring theory with $\mathcal N=1$ and $\mathcal N=2$
worldsheet supersymmetry. In the case of $\mathcal N=2$ string, we also derive
all the $\overline{SL (2,\mathbb R)}$ descendant contributions in the celestial
OPE; the soft sector of such OPE just yields the $w_{1+\infty}$ algebra after
rewriting in terms of chiral modes. Our stringy derivation of celestial OPEs
thus initiates the first step towards the realization of celestial holography
in string theory.
