I want to model a variable $y$, which is known to be affected by a set of 5 variables: {$x_1$, $x_2$, $x_3$, $x_4$, $x_5$}. These variables are known to be correlated with one another and $y$ to some extent. From the literature it also seems probable that some polynomial (e.g. $x_1^2$) and interaction effects (e.g. $x_1x_2$) may also be present, but I have no way of telling beforehand which (if any) of these 5 variables exhibits this behaviour.
What I want to do is to infer the effect of $x_1$ on $y$; specifically I want to know the coefficient and its standard error. However, with the potential polynomial and interaction terms from all 5 variables I don't think that a simple OLS fit will work here.
From what I understand the ideal model of $y$ should be both parsimonious and have a high prediction score in order for any information about $x_1$ can be derived. I've seen in the past that LASSO, Ridge, and Elastic Net regularisations have been proposed to perform feature selection on multicollinear data, but I can't calculate standard errors or p-values using these techniques because these are biased estimators. A lot of material online tackles multicollinearity from the perspective of prediction-only models, but I've seen little discussion on how inference can be performed using such data.
What should my approach here be? Can regularisation methods still be used for inference somehow? What are the best practices for statistical inference using multicollinear data?