Proof that p-value in OLS regression is symmetric OLS regression is not symmetric, meaning that it produces different relationships if you flip the dependent and independent variables; however, it would seem odd if the p-values were different and indeed empirical testing shows them to be identical:
import numpy as np
import statsmodels.api as sm
xs = np.random.normal(0, 1, 100)
ys = np.random.normal(0, 1, 100)
sm.OLS(ys, xs).fit().pvalues == sm.OLS(xs, ys).fit().pvalues

What is the proof of this fact?

This is an over-simplified answer regarding only simple linear regression. I assume you know in advance assumptions of the simple linear regression model I imply and understand notations I use as I am short in time to reply in detail. The Pearson's correlation coefficient is:
$r = \frac{\sum_1^n(x_i - \bar{x})(y_i - \bar{y})}{\sqrt{\sum_1^n (x_i - \bar{x})^2\sum_1^n (y_i - \bar{y})^2}}$
It is a fixed quantity in both models. Therefore, $R^2 = r^2$ is also a fixed quantity no matter the
equation is $y \sim x$ or $x \sim y$
Assuming linear independence hypothesis.
Therefore, $f = \frac{R^2}{1-R^2}(n-2)$ is a fixed quantity, so is $|t|$ as $f = t^2$.
$P_{value} = 2Pr(T_{n-2} >= |t|)$ will be a fixed quantity no matters what regression equation you choose.
There is definitely a better solution by deriving the t-statistics directly from the formula $t = \frac{b_1}{\sqrt{MSE}}$, but I am short in time to write it down. It is not possible to show what I showed above without going from basic formulas, so my solution is definitely not a good one, but I hope you find it convincing.
