In this paper, we introduce a natively positive approximation method based on the Feynman-Kac representation using random walks, to approximate the solution to the one-dimensional parabolic Anderson model of Skorokhod type, with either a flat or a Dirac delta initial condition. Assuming the driving noise is a fractional Brownian sheet with Hurst parameters $H \geq \frac{1}{2}$ and $H_* \geq \frac{1}{2}$ in time and space, respectively, we also provide an error analysis of the proposed method. The error in $L^p (Ω)$ norm is of order \[ O \big(h^{\frac{1}{2}[(2H + H_* - 1) \wedge 1] - ε}\big), \] where $h > 0$ is the step size in time (resp. $\sqrt{h}$ in space), and $ε> 0$ can be chosen arbitrarily small. This error order matches the Hölder continuity of the solution in time with a correction order $ε$, making it `almost' optimal. Furthermore, these results provide a quantitative framework for convergence of the partition function of directed polymers in Gaussian environments to the parabolic Anderson model.