  Recent work by Wu {\em et al.} [arXiv:1910.11011] proposed a numerical
method, so-called matrix product operator-matrix product state (MPO-MPS)
method, by which several types of quantum many-body wave functions, in
particular, the projected Fermi sea state, can be efficiently represented as a
tensor network. In this paper, we generalize the MPO-MPS method to study
Gutzwiller projected paired states of fermions, where the maximally localized
Wannier orbitals for Bogoliubov quasiparticles/quasiholes have been adapted to
improve the computational performance. The study of $SO(3)$-symmetric spin-1
chains reveals that this new method has better performance than variational
Monte Carlo for gapped states and similar performance for gapless states.
Moreover, we demonstrate that dynamic correlation functions can be easily
evaluated by this method cooperating with other MPS-based accurate approaches,
such as the Chebyshev MPS method.
