There is no "best way" to deal with selection bias, only appropriate ones given the right circumstances.
PSM would allow you to find the heterogeneous causal effect of a treatment on a group by submitting it into sub-groups. It assumes that, given that the correct division into neighborhoods was made, then being treated in each group (of homogenous subjects, some of which experienced treatment and some were not), is as good as random. In some sense, it is a customization of an RNN design, where you have "good as randomness" only in a small neighborhood of your space, for example, getting SAT scores around the minimum requirements of getting into college. Since your scores on the SAT are mainly driven by knowledge and preparation, but a small portion of it is driven by random noise (say, the air conditioning didn't work during one's test), then getting scores within some small neighborhood is almost random.
Using the pre-defined function in STATA or not doesn't mean anything, but the question is whether you did it right. From what you wrote, I couldn't find how you actually got the scores of the 1st stage - what covariates did you use to score the propensity of a subject getting treated given covariates' values?
After having scores, you need to match them, there are multiple methods (KNN, matching, distances), and there are some tests you need to run to make your your setup is actually valid for the next step, of comparing means. Again, you can do all of that manually, but it is important that you will follow the steps of the method in order for you to be sure about the treatment effects you'd find. I am copying here the steps, as written in Wikipedia.
1. Estimate propensity scores, e.g. with logistic regression.
Dependent variable: Z = 1, if unit participated (i.e. is member of the treatment group); Z = 0, if unit did not participate (i.e. is member of the control group).
Choose appropriate confounders (variables hypothesized to be associated with both treatment and outcome)
Obtain an estimation for the propensity score: predicted probability p or log[p/(1 − p)].
2. Match each participant to one or more nonparticipants on propensity score, using one of these methods:
- Nearest neighbor matching
- Optimal full matching: match each participants to unique non-participant(s) so as to minimize the total distance in propensity scores between participants and their matched non-participants. This method can be combined with other matching techniques.
- Caliper matching: comparison units within a certain width of the propensity score of the treated units get matched, where the width is generally a fraction of the standard deviation of the propensity score
- Mahalanobis metric matching in conjunction with PSM
- Stratification matching
- Difference-in-differences matching (kernel and local linear weights)
Exact matching
3. **Check that covariates are balanced across treatment and comparison groups within strata of the propensity score.**
- Use standardized differences or graphs to examine distributions
- If covariates are not balanced, return to steps 1 or 2 and modify the procedure
4. Estimate effects based on new sample
Typically: a weighted mean of within-match average differences in outcomes between participants and non-participants.
Use analyses appropriate for non-independent matched samples if more than one nonparticipant is matched to each participant
Good luck
teffects psmatch
and include all variables you want to adjust for in the propensity score model. $\endgroup$