Using the value of the precession of Earth's axis in 25,800 years and the value of Earth's radius, mass, moment of inertia etc., would it be possible to calculate/estimate that Earth's axis must be tilted approx. 23.5°?
The question can probably be posed in different ways. E.g. is the precession period related to the tilting angle?
So out of curiosity, I decided to go through some more fundamental models and their reasonable simplifications in deriving the principle behind Earth's rotational axis' precession. If you want, read on. Otherwise, you can see a final formula at the end. The long-time averaged precession can be interpreted as the result of the averaged toques that the Sun and the Moon exert on Earth, due to the slightly ellipsoidal shape of Earth (if Earth were a perfect sphere, there would be no such torque).
First, let us fix an inertial coordinate system $O\, \vec{e}_x \, \vec{e}_y \, \vec{e}_z$ at the center of the Sun, assuming that the Sun is stationary with respect to this coordinate system. The position of the Earth in $O\, \vec{e}_x \, \vec{e}_y \, \vec{e}_z$ is described by the position vector $$\vec{r} = x \vec{e}_x + y \vec{e}_y + z \vec{e}_z$$ while the position of the Moon in $O\, \vec{e}_x \, \vec{e}_y \, \vec{e}_z$ is described by the position vector $$\vec{r}_m = x_m \vec{e}_x + y_m \vec{e}_y + z_m \vec{e}_z$$ Consequently, the difference vector $\vec{\rho} = \vec{r}_m - \vec{r}$ describes the position of the Moon relative to the Earth. Furthermore, assume that the Earth has the shape of a rotational ellipsoid with a unit vector $\vec{E}_z$ aligned with its geometric ellipsoidal axis of symmetry. There is a chosen Earth-fixed coordinate system $O_E\, \vec{E}_x \, \vec{E}_y \, \vec{E}_z$, whose $z-$axis is aligned with the unit vector $\vec{E}_z$. From the point of view of the inertial frame of reference $O\, \vec{e}_x \, \vec{e}_y \, \vec{e}_z $, the Earth-fixed coordinate system $O_E\, \vec{E}_x \, \vec{E}_y \, \vec{E}_z$ moves and rotates with time and consequently the axis $\vec{E}_z = \vec{E}_z(t)$ also changes with time relative to the inertial system $O\, \vec{e}_x \, \vec{e}_y \, \vec{e}_z$.
Let $M$ be the mass of the Sun, $m$ be the mass of the Earth and $m_m$ be the mass of the Moon. According to the so called $J_2$ gravitational geopotential model, the gravitational acceleration exerted by the Earth on any mass-point $p$ with position vector $\vec{r}_{p}$ in the inertial coordinate system $O\, \vec{e}_x\, \vec{e}_y\, \vec{e}_z$, and consequently with relative position vector $\vec{\rho}_{p} = \vec{r}_{p} - \vec{r}$ in the inertial coordinate system $O\, \vec{e}_x \, \vec{e}_y \, \vec{e}_z$ is
\begin{align}
\vec{g}\big(\vec{\rho}_{p}, \, \vec{E}_z\big) \, = \, - \, \frac{Gm}{|{\rho}_{p}|^3} \, \, \vec{\rho}_{p} \, - \, \frac{3 J_2}{2} \, \frac{\big( 5 \, (\vec{\rho}_p \cdot \vec{E}_z)^2 \, - \, |\rho_p|^2\big)}{|\rho_p|^7} \,\, \vec{\rho}_p \, + \, 3J_2\, \frac{(\vec{\rho}_p \cdot \vec{E}_z)}{|\rho_p|^5}\,\,\vec{E}_z
\end{align}
which reflects the fact the Earth is mostly a sphere, but it does have a bit of an extra radius around the equator, making it a very mild rotational ellipsoid.
That is why $\vec{E}_z$ is featured in this model, to account for the fact that this geopotential field is rotationally symmetric around the axis determined by $\vec{E}_z$ because the Earth is a rotational ellipsoid with symmetry axis $\vec{E}_z$.
Observe that the Earth exerts force $ M \, \vec{g}(-\, \vec{r}, \, \vec{E}_z)$ on the Sun, so by Newton's third law, the Sun exerts a force on the Earth $- \, M \, \vec{g}( - \vec{r}, \, \vec{E}_z)$. The point at which the force $-\, M\, \vec{g}\big(-\vec{r}, \, \vec{E}_z\big)$ is applied is the center of the Sun $O$, because the Sun is the one acting on the Earth. Therefore, since the vector pointing from the center of the Earth $O_E$ to the center of the Sun $O$ is $- \, \vec{r}$, the torque exerted by the Sun on the Earth is $$ \text{Torque}_{\text{ Sun}} \, = \, M\, \big( - \vec{r} \,\big) \times \big( - \vec{g}\big(- \vec{r}, \, \vec{E}_z\,\big) \,\big)\, = \, M\, \vec{r} \,\times\, \vec{g}\big(-\vec{r}, \, \vec{E}_z\,\big)$$ Analogously, the Earth exerts force $ m_m \, \vec{g}(\vec{\rho}, \, \vec{E}_z)$ on the Moon, so by Newton's third law, the Moon exerts a force on the Earth $- \, m_m \, \vec{g}( \vec{\rho}, \, \vec{E}_z)$. The point at which the force $-\, m_m\, \vec{g}\big(\vec{\rho}, \, \vec{E}_z\big)$ is applied is the center of the Moon, call it $O_m$, because the Moon is the one acting on the Earth. Therefore, since the vector pointing from the center of the Earth $O_E$ to the center of the Moon $O_m$ is $\vec{\rho}$, the torque exerted by the Moon on the Earth is $$ \text{Torque}_{\text{ Moon}} \, = \, m_m\, \vec{\rho} \times \big( - \vec{g}\big(\,\vec{\rho}, \, \vec{E}_z\,\big) \,\big) \, = \, - \, m_m\,\, \vec{\rho} \, \times\, \vec{g}\big(\,\vec{\rho}, \, \vec{E}_z\,\big) $$ The force that the Sun exerts on the Moon is assumed to be the usual symmetric Newtonian force, which assumes that both the Sun and the Moon are spherical bodies with uniformly distributed mass \begin{align} - \, \frac{Gm_m M}{|r_m|^3} \, \,{\vec{r}_m} \end{align} and because we already defined $\vec{\rho} = \vec{r}_m - \vec{r}$, we have $$\vec{r}_m = \vec{r} + \vec{\rho}$$ and therefore the latter force is \begin{align} - \, \frac{Gm_m M}{|\vec{r} + \vec{\rho}|^3} \, \,\big(\, \vec{r} + \vec{\rho} \,\big) \end{align} Consequently, the system of ordinary differential equations of motion of this system, for the unknown trajectory vector functions $$\big(\,\vec{r} = \vec{r}(t), \,\, \vec{\rho} = \vec{\rho}(t), \,\, \vec{E}_z = \vec{E}_z(t) \,\big)$$ which interpret the Earth as a rotational ellipsoid rather than as a uniform mass sphere, are \begin{align} m\, &\frac{d^2 \vec{r}}{dt^2} \, = \, - \, M\,\vec{g}\big(-\vec{r}, \, \vec{E}_z\,\big) \, - \,m_m \, \vec{g}\big(\,\vec{\rho}, \, \vec{E}_z\,\big) \\ m_m\, & \frac{d^2}{dt^2} \big(\,\vec{r} + \vec{\rho}\,\big) \, = \, \,m_m \, \vec{g}\big(\,\vec{\rho}, \, \vec{E}_z\,\big) \, - \, \frac{Gm_m M}{|\vec{r} + \vec{\rho}|^3} \, \,\big(\, \vec{r} + \vec{\rho} \,\big)\\ \frac{d}{dt} &I\, \vec{\omega} \, = \, M\,\, \vec{r} \,\times\, \vec{g}\big(-\vec{r}, \, \vec{E}_z\,\big) \,\, - \,\, m_m\,\, \vec{\rho} \, \times\, \vec{g}\big(\,\vec{\rho}, \, \vec{E}_z\,\big) \end{align}
where $\vec{\omega}$ is the angular velocity of Earth and $I$ is its inertia tensor. For the purpose of this analysis we are going to observe that in terms of an inertia tensor, Earth is pretty much a sphere and if one assumes that it rotates around its axis with the same averaged constant rate $\omega_0$, one could assume that $I$ is a constant number and $$\vec{\omega} \, = \, \omega_0 \, \vec{E}_z$$ Thus, the equations of motion become \begin{align} m\, &\frac{d^2 \vec{r}}{dt^2} \, = \, - \, M\,\vec{g}\big(-\vec{r}, \, \vec{E}_z\,\big) \, - \,m_m \, \vec{g}\big(\,\vec{\rho}, \, \vec{E}_z\,\big) \\ m_m\, &\,\frac{d^2\vec{\rho}}{dt^2} + m_m \,\frac{d^2\vec{r}}{dt^2} \, = \, \,m_m \, \vec{g}\big(\,\vec{\rho}, \, \vec{E}_z\,\big) \, - \, \frac{Gm_m M}{|\vec{r} + \vec{\rho}|^3} \, \,\big(\, \vec{r} + \vec{\rho} \,\big)\\ I\omega_0 &\frac{d\vec{E}_z}{dt} \, = \, M\,\, \vec{r} \,\times\, \vec{g}\big(-\vec{r}, \, \vec{E}_z\,\big) \,\, - \,\, m_m\,\, \vec{\rho} \, \times\, \vec{g}\big(\,\vec{\rho}, \, \vec{E}_z\,\big) \end{align} Since $\vec{r} \times (\lambda \, \vec{r}) = \vec{0}$ for any scalar number $\lambda$, the equations of motion can be written as follows: \begin{align} m\, &\frac{d^2 \vec{r}}{dt^2} \, = \, - \, M\,\vec{g}\big(-\vec{r}, \, \vec{E}_z\,\big) \, - \,m_m \, \vec{g}\big(\,\vec{\rho}, \, \vec{E}_z\,\big) \\ m_m\, &\,\frac{d^2\vec{\rho}}{dt^2} + m_m \,\frac{d^2\vec{r}}{dt^2} \, = \, \,m_m \, \vec{g}\big(\,\vec{\rho}, \, \vec{E}_z\,\big) \, - \, \frac{Gm_m M}{|\vec{r} + \vec{\rho}|^3} \, \,\big(\, \vec{r} + \vec{\rho} \,\big)\\ I\omega_0 &\frac{d\vec{E}_z}{dt} \, = \, - \, 3MJ_2\, \frac{(\vec{r} \cdot \vec{E}_z)}{|r|^5}\, \Big(\, \vec{r} \,\times\, \vec{E}_z \, \Big)\, - \, 3m_mJ_2\, \frac{(\vec{\rho} \cdot \vec{E}_z)}{|\rho|^5}\, \Big(\, \vec{\rho} \,\times\, \vec{E}_z \, \Big) \end{align} Let us simplify our model a bit further. If you look at the forces produced by an acceleration of the form $\vec{g}(\vec{\rho}_p, \, \vec{E}_z)$ you will notice that the only term without a $J_2$ constant is the usual classical Newtonian force. Now, the rest of the terms are divided by at least $|\rho_p|^4$. Therefore, the $J_2$ effects on Earth's and the Moon's orbital motions are so weak, that they are practically zero, so the equations of motion simplify to \begin{align} m\, &\frac{d^2 \vec{r}}{dt^2} \, = \, - \, \frac{GmM}{|r|^3} \, \vec{r} \, + \, \frac{Gm_mm}{|\rho|^3} \, \vec{\rho} \\ m_m\, &\,\frac{d^2\vec{\rho}}{dt^2} + m_m \,\frac{d^2\vec{r}}{dt^2} \, = \, -\,\frac{Gm_mm}{|\rho|^3} \, \vec{\rho} \, - \, \frac{Gm_m M}{|\vec{r} + \vec{\rho}|^3} \, \,\big(\, \vec{r} + \vec{\rho} \,\big)\\ I\omega_0 &\frac{d\vec{E}_z}{dt} \, = \, - \, 3MJ_2\, \frac{(\vec{r} \cdot \vec{E}_z)}{|r|^5}\, \Big(\, \vec{r} \,\times\, \vec{E}_z \, \Big)\, - \, 3m_mJ_2\, \frac{(\vec{\rho} \cdot \vec{E}_z)}{|\rho|^5}\, \Big(\, \vec{\rho} \,\times\, \vec{E}_z \, \Big) \end{align} Observe that after this realistic model simplification, the first two vector ordinary differential equations are decoupled from the last rigid-body vector differential equation, which on the other hand, is coupled to the first two equations.
Let us first rewrite the latter system as follows: \begin{align} &\frac{d^2 \vec{r}}{dt^2} \, = \, - \, \frac{GM}{|r|^3} \, \vec{r} \, + \, \frac{Gm_m}{|\rho|^3} \, \vec{\rho} \\ &\,\frac{d^2\vec{\rho}}{dt^2} \, = \, -\,\frac{Gm}{|\rho|^3} \, \vec{\rho} \, - \, \frac{G M}{|\vec{r} + \vec{\rho}|^3} \, \,\big(\, \vec{r} + \vec{\rho} \,\big) \, - \, \frac{d^2\vec{r}}{dt^2}\\ I\omega_0 &\frac{d\vec{E}_z}{dt} \, = \, - \, 3MJ_2\, \frac{(\vec{r} \cdot \vec{E}_z)}{|r|^5}\, \Big(\, \vec{r} \,\times\, \vec{E}_z \, \Big)\, - \, 3m_mJ_2\, \frac{(\vec{\rho} \cdot \vec{E}_z)}{|\rho|^5}\, \Big(\, \vec{\rho} \,\times\, \vec{E}_z \, \Big) \end{align} Next, use the first equation in the second one \begin{align} &\frac{d^2 \vec{r}}{dt^2} \, = \, - \, \frac{GM}{|r|^3} \, \vec{r} \, + \, \frac{Gm_m}{|\rho|^3} \, \vec{\rho} \\ &\,\frac{d^2\vec{\rho}}{dt^2} \, = \, -\,\frac{Gm}{|\rho|^3} \, \vec{\rho} \, - \, \frac{G M}{|\vec{r} + \vec{\rho}|^3} \, \,\big(\, \vec{r} + \vec{\rho} \,\big) \, + \, \frac{GM}{|r|^3} \, \vec{r} \, - \, \frac{Gm_m}{|\rho|^3} \, \vec{\rho}\\ I\omega_0 &\frac{d\vec{E}_z}{dt} \, = \, - \, 3MJ_2\, \frac{(\vec{r} \cdot \vec{E}_z)}{|r|^5}\, \Big(\, \vec{r} \,\times\, \vec{E}_z \, \Big)\, - \, 3m_mJ_2\, \frac{(\vec{\rho} \cdot \vec{E}_z)}{|\rho|^5}\, \Big(\, \vec{\rho} \,\times\, \vec{E}_z \, \Big) \end{align} and after regrouping \begin{align} &\frac{d^2 \vec{r}}{dt^2} \, = \, - \, \frac{GM}{|r|^3} \, \vec{r} \, + \, \frac{Gm_m}{|\rho|^3} \, \vec{\rho} \\ &\,\frac{d^2\vec{\rho}}{dt^2} \, = \, -\,\frac{G m}{|\rho|^3} \, \vec{\rho} \, + \, \left(\frac{G M}{|{r}|^3} \, - \, \frac{G M}{|\vec{r} + \vec{\rho}|^3}\right) \, \vec{r} \, - \, \frac{GM}{|\vec{r} + \vec{\rho}|^3} \, \vec{\rho} \, -\,\frac{Gm_m}{|\rho|^3} \, \vec{\rho}\\ I\omega_0 &\frac{d\vec{E}_z}{dt} \, = \, - \, 3MJ_2\, \frac{(\vec{r} \cdot \vec{E}_z)}{|r|^5}\, \Big(\, \vec{r} \,\times\, \vec{E}_z \, \Big)\, - \, 3m_mJ_2\, \frac{(\vec{\rho} \cdot \vec{E}_z)}{|\rho|^5}\, \Big(\, \vec{\rho} \,\times\, \vec{E}_z \, \Big) \end{align}
As it stand up to now, the first two equations from the latter system model a type of a planar three-body problem. However, we will make another reasonable simplification, ignoring small perturbation-like terms.
For example, we can safely take the coefficient $\frac{G M}{|{r}|^3} \, - \, \frac{G M}{|\vec{r} + \vec{\rho}|^3} \approx 0$, because $|\vec{r} + \vec{\rho}| \approx |\vec{r}|$ based on the fact that the distance $|\rho|$ between the Moon and the Earth is miniscule compared to the distance between the Sun and the Earth $|r|$. For the same reason $\frac{GM}{|\vec{r} + \vec{\rho}|^3}\, \vec{\rho} \approx 0$. Finally, the acceleration $\frac{G m_m}{|\rho|^3} \, \vec{\rho} \, \approx \, 0$ compared to the other accelerations, because the Moon's mass is much smaller than the mass of the Earth and the Sun (Moon's mass is 0.012300 of Earth's). As a result of all of these simplifications of our model, we end up with the following easier to analyse system \begin{align} &\frac{d^2 \vec{r}}{dt^2} \, = \, - \, \frac{GM}{|r|^3} \, \vec{r} \\ &\,\frac{d^2\vec{\rho}}{dt^2} \, = \, -\,\frac{G m}{|\rho|^3} \, \vec{\rho}\\ I\omega_0 &\frac{d\vec{E}_z}{dt} \, = \, - \, 3MJ_2\, \frac{(\vec{r} \cdot \vec{E}_z)}{|r|^5}\, \Big(\, \vec{r} \,\times\, \vec{E}_z \, \Big)\, - \, 3m_mJ_2\, \frac{(\vec{\rho} \cdot \vec{E}_z)}{|\rho|^5}\, \Big(\, \vec{\rho} \,\times\, \vec{E}_z \, \Big) \end{align} Observe now that the first equation is completely decoupled from the second and the third equation, and the second equation is decoupled from the first and the third equation. Both, the first and the second equations describes the very classical Kepler's problem, where its solutions $\vec{r} = \vec{r}(t)$ and $\vec{\rho} = \vec{\rho}(t)$ follow motions along elliptical orbits. Therefore, we could be able to solve the first and the second equations and find the motions along the elliptical orbits. After that we can plug them in the third equation and focus on its solution.
As a sanity check, observe that \begin{align} \left( \vec{E}_z \cdot \frac{d \vec{E}_z}{dt}\right) \, =&\, \left(- \, \frac{3MJ_2}{I\omega_0} \,\, \frac{(\vec{r} \cdot \vec{E}_z)}{|r|^5} \right) \, \Big( \vec{E}_z \cdot \big(\, \vec{r} \,\times\, \vec{E}_z \, \big)\, \Big)\, \\& \,\,\, + \, \left(- \, \frac{3m_mJ_2}{I\omega_0} \,\, \frac{(\vec{\rho} \cdot \vec{E}_z)}{|\rho|^5} \right) \, \Big( \vec{E}_z \cdot \big(\, \vec{\rho} \,\times\, \vec{E}_z \, \big)\, \Big) \, = \, 0 \end{align} Therefore $$\frac{d}{dt} \big| \vec{E}_z\big|^2 = 2 \, \left( \vec{E}_z \cdot \frac{d \vec{E}_z}{dt}\right) \, = \, 0$$ which means that $\big| \vec{E}_z(t) \big|^2 \, = \, \big| \vec{E}_z(0) \big|^2 \, = \, 1 $ for all $t$, i.e. the vector $\vec{E}_z(t)$ stays always unit.
To understand the long-term motion of Earth's rotational axis $\vec{E}_z = \vec{E}_z(t)\,$, let $\vec{r} = \vec{r}(t)$ be the time-evolution of the Earth along its elliptical orbit around the Sun, where the Sun's center coincides with the origin of the inertial coordinate system $O \,\vec{e}_x\, \,\vec{e}_y\,\,\vec{e}_z$ and is at the same time one of the two foci of the orbital ellipse traversed by $\vec{r} = \vec{r}(t)$. Let $\vec{\rho} = \vec{\rho}(t)$ be the solution to the second equation. The vector $\vec{\rho}$ points from the Earth to the Moon and, again by the fact that it is a solution to Kepler's problem, one can interpret it as a position vector in a coordinate system $O_E \,\vec{e}_z \, \vec{e}_y \, \vec{e}_z$ that has the center of the Earth $O_E$ as origin $O_E$ and axes always parallel to the axes of the inertial coordinate coordinate system $O \,\vec{e}_z \, \vec{e}_y \, \vec{e}_z$. Then, we are left with analysing the third differential equation \begin{align} &\frac{d \vec{E}_z}{dt} \, = \, \left(- \, \frac{3MJ_2}{I \omega_0}\right) \, \frac{\big(\,\vec{r}(t) \cdot \vec{E}_z\,\big)}{|r(t)|^5}\, \big(\, \vec{r}(t) \,\times\, \vec{E}_z \, \big) \,\, + \,\, \left(- \, \frac{3m_mJ_2}{I \omega_0}\right) \, \frac{\big(\,\vec{\rho}(t) \cdot \vec{E}_z\,\big)}{|\rho(t)|^5}\, \big(\, \vec{\rho}(t) \,\times\, \vec{E}_z \, \big) \end{align} or written in another form \begin{align} \frac{d \vec{E}_z}{dt} \, = \,& \left(- \, \frac{3MJ_2}{I \omega_0}\right) \, \frac{1}{|r(t)|^3}\,\left( \frac{\big(\,\vec{r}(t) \cdot \vec{E}_z\,\big)}{|r(t)|^2}\, \vec{r}(t) \right) \times \vec{E}_z \\ &\, + \, \left(- \, \frac{3m_mJ_2}{I \omega_0}\right) \, \frac{1}{|\rho(t)|^3}\,\left( \frac{\big(\,\vec{\rho}(t) \cdot \vec{E}_z\,\big)}{|\rho(t)|^2}\, \vec{\rho}(t) \right) \times \vec{E}_z \end{align} Observe that the expressions $$P(\vec{r}) \, \vec{E}_z \, = \, \frac{\big(\,\vec{r} \cdot \vec{E}_z\,\big)}{|r|^2}\,\, \vec{r} \,\,\,\,\, \text{ and } \,\,\,\,\,P(\vec{\rho}) \, \vec{E}_z \, = \, \frac{\big(\,\vec{\rho} \cdot \vec{E}_z\,\big)}{|\rho|^2}\,\, \vec{\rho}$$ define the orthogonal projection operators of an arbitrary vector $\vec{E}_z$ onto the axes determined by the vectors $\vec{r}$ and $\vec{\rho}$ respectively. Hence, the equation of motion becomes \begin{align} \frac{d \vec{E}_z}{dt} \, = \,& \left(- \, \frac{3MJ_2}{I \omega_0}\right) \, \frac{1}{|r(t)|^3}\,\left( P\big(\vec{r}(t)\big) \, \vec{E}_z \right) \times \vec{E}_z \\ &+ \, \left(- \, \frac{3m_mJ_2}{I \omega_0}\right) \, \frac{1}{|\rho(t)|^3}\,\left( P\big(\vec{\rho}(t)\big) \, \vec{E}_z \right) \times \vec{E}_z \end{align} Since the Earth's motion around the Sun takes approximately $T_e = 365.256$ days to complete one revolution around the Sun, and the Moon's motion around the Earth (in the coordinate system $O_E \, \vec{e}_x \, \vec{e}_y \, \vec{e}_z$) takes approximately $T_m = 27.322$ days to complete one revolution around the Earth, the latter equation of motion can be considered $T-$periodic, where $T = n_m\, T_m = n_e \, T_e$ for two appropriate positive integers $n_m$ and $n_e$. According to the averaging principle (which for one single averaged variable, $t$ in our case, is actually a theorem), the long-time-evolution of Earth's rotational axis $\vec{E}_z = \vec{E}_z(t)$ follows very closely the solution of the averaged equation \begin{align} &\frac{d \vec{E}_z}{dt} \, = \, \left(- \, \frac{3MJ_2}{I \omega_0}\right) \, \left( Q \, \vec{E}_z \right) \times \vec{E}_z \,\, + \,\, \left(- \, \frac{3m_mJ_2}{I \omega_0}\right) \, \left( Q_m \, \vec{E}_z \right) \times \vec{E}_z \end{align} oscillating around it with very small amplitude. In the averaged equation \begin{align} Q &= \frac{1}{T}\int_{0}^{T} \frac{1}{|r(t)|^3} \, P\big(\vec{r}(t)\big) \, dt\\ Q_e &= \frac{1}{T}\int_{0}^{T} \frac{1}{|\rho(t)|^3} \, P\big(\vec{\rho}(t)\big) \, dt \end{align} In order to calculate the averaged linear operators $Q$ and $Q_m$, we are going to look more closely at the solutions of the Kepler's problems $$\frac{d^2 \vec{r}}{dt^2} \, = \, - \, \frac{GM}{|r|^3} \, \vec{r} \,\,\,\,\, \text{ and } \,\,\,\,\, \frac{d^2 \vec{\rho}}{dt^2} \, = \, - \, \frac{Gm}{|\rho|^3} \, \vec{\rho}$$ which are the first and the second equation of motion in our system, each of them decoupled from the rest of the system.
Without loss of generality, we can assume that the inertial coordinate system $O\, \vec{e}_x\, \vec{e}_y\, \vec{e}_z$ is chosen so that the Earth's motion $\vec{r} = \vec{r}(t)$, which is a solution to the first differential equation above, is on the $O\, \vec{e}_x\, \vec{e}_y$ plane. In order to simplify the problem a bit further, let us assume that the Moon's orbit is also in the same plane as the Earth's (i.e. both of them move on the $O\, \vec{e}_x\, \vec{e}_y$ plane, the ecliptic). In reality, this is not quite the case. Moon's orbit is $5.1^{\circ}$ tilted with respect to the ecliptic $O\, \vec{e}_x\, \vec{e}_y$. However, this tilt is small enough to be ignored for this calculation, so both position vectors $\vec{r} = \vec{r}(t)$ and $\vec{\rho} = \vec{\rho}(t)$ traverse ellipses in the $O\, \vec{e}_x\, \vec{e}_y$ ecliptic plane.
In what follows, we are going to demonstrate the calculation of the operator $Q$ for the Keplerian motion $\vec{r} = \vec{r}(t)$. The calculation of the other operator $Q_m$ for the Keplerian motion $\vec{\rho} = \vec{\rho}(t)$ are absolutely analogous.
Without loss of generality, we can assume that the coordinate axis $O\, \vec{e}_x$ is aligned with the major semi-axis of Earth's elliptical orbit. Such a choice of coordinate system is sometimes called perifocal coordinate system.
Let us express any arbitrary position vector $\vec{r}$ from the $O\, \vec{e}_x\, \vec{e}_y$ plane in polar coordinates $$\vec{r} \, = \, r \cos(\varphi) \, \vec{e}_x \, + \, r \sin(\varphi) \, \vec{e}_y$$ where $r = |\vec{r}|$. Then Kepler's problem $$\frac{d^2 \vec{r}}{dt^2} \, = \, - \, \frac{GM}{|r|^3} \, \vec{r}$$ reduces to \begin{align} \frac{d\varphi}{dt} \, &= \, \frac{2\, \pi b}{T_ea(1-e^2)^2} \,\, \big( 1 \, + \, e \, \cos(\varphi)\big)^2\\ r \, &= \, \frac{a(1-e^2)}{1 \, + \,e\, \cos(\varphi)} \end{align} (this is a standard derivation when solving Kepler's problem). Here, $a$ is the major semi-axis of the elliptical orbit, $b$ is the minor semi-axis, $e$ is the eccentricity of the orbit, $b = a \sqrt{1-e^2}$, and $T_e = 2\,\pi \, \sqrt{\frac{a^3}{GM}}$ is the orbit's period.
For our purposes, in order to to calculate the averaged torque acting on Earth's axis of rotation $\vec{E}_z$, we can rewrite the latter pair of equations as \begin{align} \frac{d\varphi}{dt} \, &= \, \frac{2\, \pi ab}{T_e} \, \frac{1}{r^2}\\ r \, &= \, \frac{a(1-e^2)}{1 \, + \,e\, \cos(\varphi)} \end{align} In fact, the first equation is the conservation of the angular momentum of Earth's orbital motion, while the second one comes from adding to the latter the conservation of the Runge-Lenz vector of the orbit. Since $\varphi = \varphi(t)$ is a montonic function, it can be inverted as $t = t(\varphi)$ and thus $\vec{r} = \vec{r}(\varphi)$ and ${r} = {r}(\varphi)$. This allows us to change the integration variable in the averaging integral \begin{align} Q \, &= \, \frac{1}{T}{\int_{0}}^{T} \, \frac{1}{|r(t)|^3} \, P\big(\vec{r}(t)\big) \, dt \, = \, \frac{1}{m_e T_e}\int_{0}^{m_e T_e} \, \frac{1}{|r(t)|^3} \, P\big(\vec{r}(t)\big) \, dt\\ &= \, \frac{m_e}{m_e T_e}{\int_{0}}^{T_e} \, \frac{1}{|r(t)|^3} \, P\big(\vec{r}(t)\big) \, dt \, = \, \frac{1}{T_e}\int_{0}^{T_e} \, \frac{1}{|r(t)|^3} \, P\big(\vec{r}(t)\big) \, dt\\ &= \, \frac{1}{T_e}\int_{0}^{2\pi} \, \frac{1}{r(\varphi)^3} \, P\big(\vec{r}(\varphi)\big) \, \left({\frac{d\varphi}{dt}}\right)^{-1} d\varphi\\ &= \, \frac{1}{T_e}\int_{0}^{2\pi} \, \frac{1}{r(\varphi)^3} \, P\big(\vec{r}(\varphi)\big) \, \left({\frac{T_e \,\, r(\varphi)^2}{2 \pi \, ab}}\right) d\varphi\\ &= \, \frac{1}{2 \pi \, ab}\int_{0}^{2\pi} \, \frac{1}{r(\varphi)} \, P\big(\vec{r}(\varphi)\big) \, d\varphi \end{align} Let us write the $\varphi$-family of projection $P\big(\vec{r}(\varphi)\big)$ operators more explicitly. For any arbitrary vector $\vec{v} = v_x \vec{e}_x + v_y \vec{e}_y + v_z \vec{e}_z$ \begin{align} P\big(\vec{r}(\varphi)\big) \, \vec{v} \, &= \, \frac{\big(\,\vec{r}(\varphi)\cdot \vec{v}\,\big)}{r(\varphi)^2} \,\, \vec{r}(\varphi)\\ &= \, \frac{\big(\, {r}(\varphi)\cos(\varphi)\vec{e}_x + {r}(\varphi)\sin(\varphi)\vec{e}_y\,\big)\cdot \vec{v}\,}{r(\varphi)^2} \, \, \vec{r}(\varphi) \\ &= \, \frac{\cos(\varphi)\big(\vec{e}_x \cdot \vec{v}\big) + \sin(\varphi)\big(\vec{e}_y \cdot \vec{v}\big)\,}{r(\varphi)} \, \, \big(\, {r}(\varphi)\cos(\varphi)\vec{e}_x + {r}(\varphi)\sin(\varphi)\vec{e}_y\,\big) \\ &= \, \big( v_x \cos(\varphi) + v_y \sin(\varphi)\,\big) \, \, \big(\, \cos(\varphi)\vec{e}_x + \sin(\varphi)\vec{e}_y\,\big) \end{align} so we can write explicitly the $\varphi-$family of projectors as \begin{align} P\big(\vec{r}(\varphi)\big) \, \vec{v} \, &= \, \big( v_x \cos^2(\varphi) + v_y \cos(\varphi) \sin(\varphi)\,\big) \vec{e}_x \, + \, \big( v_x \cos(\varphi)\sin(\varphi) + v_y \sin^2(\varphi)\,\big) \vec{e}_y \end{align} From this latter expression, we can see that there are three integrals that need to be computed: \begin{align} Q_{11} &= \frac{1}{2\pi \, ab} \, \int_{0}^{2\pi} \, \frac{\cos^2(\varphi)}{r(\varphi)} \, d\varphi \\ Q_{12} &= \frac{1}{2\pi \, ab} \, \int_{0}^{2\pi} \, \frac{\cos(\varphi)\sin(\varphi)}{r(\varphi)} \, d\varphi \\ Q_{22} &= \frac{1}{2\pi \, ab} \, \int_{0}^{2\pi} \, \frac{\sin^2(\varphi)}{r(\varphi)} \, d\varphi \end{align} Since $r = r(\varphi) = \frac{a(1-e^2)}{1 \, + \, e \cos(\varphi)}$ the integrals explicitly look like \begin{align} Q_{11} &= \frac{1}{2\pi \, a^2b\,(1-e^2)} \, \int_{0}^{2\pi} \, \cos^2(\varphi)\big( \, 1 \, + \, e \cos(\varphi) \, \big) \, d\varphi \\ Q_{12} &= \frac{1}{2\pi \, a^2b\,(1-e^2)} \, \int_{0}^{2\pi} \, \cos(\varphi)\sin(\varphi)\big( \, 1 \, + \, e \cos(\varphi) \, \big) \, d\varphi \\ Q_{22} &= \frac{1}{2\pi \, a^2b\,(1-e^2)} \, \int_{0}^{2\pi} \, \sin^2(\varphi)\big( \, 1 \, + \, e \cos(\varphi) \, \big)\, d\varphi \end{align} and can be easily integrated \begin{align} Q_{11} &= \frac{1}{2\pi \, b^3} \, \int_{0}^{2\pi} \, \cos^2(\varphi)\big( \, 1 \, + \, e \cos(\varphi) \, \big) \, d\varphi = \frac{\pi}{2\pi \, b^3}\\ Q_{12} &= \frac{1}{2\pi \, b^3} \, \int_{0}^{2\pi} \, \cos(\varphi)\sin(\varphi)\big( \, 1 \, + \, e \cos(\varphi) \, \big) \, d\varphi = 0\\ Q_{22} &= \frac{1}{2\pi \, b^3} \, \int_{0}^{2\pi} \, \sin^2(\varphi)\big( \, 1 \, + \, e \cos(\varphi) \, \big) \, d\varphi = \frac{\pi}{2\pi \, b^3} \end{align} and thus \begin{align} Q \, \vec{v} \, &= \, \frac{1}{T}\int_{0}^{T} \, \frac{1}{|r(t)|^3} \, P\big(\vec{r}(t)\big) \, \vec{v} \, dt \\ &= \, \frac{1}{2 \pi \, ab}\int_{0}^{2\pi} \, \frac{1}{r(\varphi)} \, P\big(\vec{r}(\varphi)\big) \,\vec{v} \, d\varphi\\ &= \, \frac{1}{2b^3} \big( \, v_x \vec{e}_x \, + \, v_y \vec{e}_y \, \big) \end{align} So if we denote by $P$ the orthogonal projector onto the coordinate plane $O \, \vec{e}_z \, \vec{e}_y$ $$P \, \vec{v} = P \big( v_x \vec{e}_x + v_y \vec{e}_y + v_z \vec{e}_z \big) = v_x \vec{e}_x + v_y \vec{e}_y$$ we can write down the averaged operator as $$Q \, = \, \frac{1}{2b^3} \, P$$ As already stated earlier, using analogous calculations $$Q_m \, = \, \frac{1}{2b_m^3} \, P$$ where $b_m$ is the semi-minor axis of Moon's orbit. Consequently, the averaged differential equation for Earth's rotational axis is \begin{align} &\frac{d \vec{E}_z}{dt} \, = \, \left(- \, \frac{3MJ_2}{2b^3\, I \omega_0}\right) \, \left( P \, \vec{E}_z \right) \times \vec{E}_z \, + \, \left(- \, \frac{3m_mJ_2}{2b_m^3\, I \omega_0}\right) \, \left( P \, \vec{E}_z \right) \times \vec{E}_z \end{align} or after regrouping \begin{align} &\frac{d \vec{E}_z}{dt} \, = \, - \, \nu_0 \left( P \, \vec{E}_z \right) \times \vec{E}_z \end{align} where $$\nu_0 \, = \, \frac{3MJ_2}{2b^3\, I \omega_0} \, + \, \frac{3m_mJ_2}{2b_m^3\, I \omega_0} \, = \, \frac{3\, J_2}{2\, I \omega_0}\left(\, \frac{M}{b^3} \, + \, \frac{m_m}{b_m^3}\,\right)$$ Set $\vec{E}_{x,y} = v_x \vec{e}_x + v_y \vec{e}_y$. Then $\vec{E}_z = \vec{E}_{x,y} + E \vec{e}_z$ and by the definition of an orthogonal projector, $P \, \vec{E}_z = \vec{E}_{x,y}$. The averaged equation is \begin{align} \frac{d \vec{E}_{x,y}}{dt} \, + \, \frac{dE}{dt}\,\vec{e}_z \, &= \, -\, \nu_0 \, \left( \vec{E}_{x,y} \right) \times \big(\, \vec{E}_{x,y} \, + \, E \, \vec{e}_z\,\big)\\ &= \, - \nu_0 \, \left( \vec{E}_{x,y} \right) \times \big(\, E \, \vec{e}_z\,\big)\\ &= \, E\, \nu_0 \, \big(\, \vec{e}_z \times \vec{E}_{x,y} \,\big) \end{align} The cross-product $\vec{e}_z \times \vec{E}_{x,y}$ is simply the counter-clockwise $90^{\circ}-$rotation of $\vec{E}_{x,y}$, where $\vec{E}_{x,y}$ is rotated in the $O \, \vec{e}_x \, \vec{e}_y$ plane and is therefore still orthogonal to $\vec{e}_z$. Consequently, \begin{align} \frac{d \vec{E}_{x,y}}{dt} \, &= \, E\, \nu_0 \, \big(\, \vec{e}_z \times \vec{E}_{x,y} \,\big) \\ \frac{dE}{dt} \, &= \, 0 \end{align} Therefore, the tilt of Earth's rotational axis is measured as the angle $\varepsilon$ between the unit axis-vector $\vec{E}_z$ and the unit vector $\vec{e}_z$. Therefore $E = \cos(\varepsilon)$ is constant and thus $$\vec{E}_z(t) \, = \, \vec{E}_{x, y}(t) \, + \, \cos(\varepsilon) \,\vec{e}_z$$ where $\vec{E}_{x, y} = \vec{E}_{x, y}(t)$ is the 2D vector in $O \, \vec{e}_x \, \vec{e}_y$ that satisfies the equation \begin{align} \frac{d \vec{E}_{x,y}}{dt} \, &= \, \nu_0 \cos(\varepsilon) \,\, \big(\, \vec{e}_z \times \vec{E}_{x,y} \,\big) \end{align} This is a linear vector equation with constant coefficients, so its solutions can be readily found $$\vec{E}_{x, y} \, = \, \Big(\, E^0_x \, \cos\big(\nu_{\varepsilon} \, t \big) \, - \, E^0_y \, \sin\big(\nu_{\varepsilon} \, t \big) \,\Big) \vec{e}_x \,\, + \,\, \Big(\, E^0_x \, \sin\big(\nu_{\varepsilon} \, t \big) \, + \, E^0_y \, \cos\big(\nu_{\varepsilon} \, t \big) \,\Big) \vec{e}_y$$ where $$\nu_{\varepsilon} \, = \, \nu_0 \cos(\varepsilon)\, = \, \frac{3\, J_2}{2\, I \omega_0}\left(\, \frac{M}{b^3} \, + \, \frac{m_m}{b_m^3}\,\right)\cos(\varepsilon) $$ and $E^0_x, \, E^0_y$ are initial constants. Finally, the preiod $T_{\text{prec}} = \frac{2\pi}{\nu_{\varepsilon}}$ of the solution, which is the period of precession, is related to the angle of tilt of Earth's rotation axis $\vec{E}_z$ relative to the ecliptic plane $O \, \vec{e}_x \, \vec{e}_y$ of Earth's orbit around the Sun, by the equation $$\frac{2\pi}{T_{\text{prec}}} = \frac{3\, J_2}{2\, I \omega_0}\left(\, \frac{M}{b^3} \, + \, \frac{m_m}{b_m^3}\,\right)\cos(\varepsilon) $$ The coefficient $$J_2 = G\,(I - I_{\text{eq}})$$ where $I \approx 0.330\, m \, R_{\text{mean}}^2$ is Earth's moment of inertia along the axis determined by $\vec{E}_z$ and $I_{\text{eq}} \approx 0.329\, m \, R_{\text{mean}}^2$ is the moment of inertia along an axis in Earth's equtorial plane. $R_{\text{mean}}$ is Earth's mean radius. $G$ is the gravitational constant. The equation above becomes $$\frac{2\pi}{T_{\text{prec}}} = \frac{3\, G (I - I_{\text{eq}})}{2\, I \omega_0}\left(\, \frac{M}{a^3(1-e^2)^{\frac{3}{2}}} \, + \,\frac{m_m}{a_m^3(1-e_m^2)^{\frac{3}{2}}}\,\right)\cos(\varepsilon) $$ Now, after having figured out how the averaged precession of Earth's rotation axis works, we simply have to start plugging in constants and calculating the quantities in the latter formula to calculate the period $T_{\text{prec}}$ from the angle of tilt $\varepsilon$ or vice versa.
The precession of the earth's tilt is caused primarily by the torque from the moon pulling on the equatorial bulge (with periodic help from the sun). The pull of each is stronger on the near side bulge. The torque varies throughout a month and year and will be strongest when the earth has its maximum tilt toward the moon and sun. This gets complicated because the moon's orbit is tilted relative to that of the earth.
First some general remarks:
In real life astronomy it is the moment of inertia that is known to the least level of precision.
Celestial bodies don't have a uniform density. If the goal is to calculate the moment of inertia taking the the density as a function of distance-to-the-center into account, then you need those density data.
In real life astronomy: whenever necessary the moment of inertia is inferred from the rate of precession.
Journal article:
Information on internal structure from shape, gravity field and rotation
As to your question specifically: Let's use the name 'C' for any celestial body. The rate of precession arises from the magnitude of the torque that C is subject to. The magnitude of the tilt does affect the magnitude of the torque that C is subject to.
A possible snag: what if there are, for a given rate of precession, two different tilts that will give rise to that rate of precession? Just of the top of my head I can't exclude that possibility.
With that caveat: it is the case that the magnitude of the tilt affects the torque that C is subject to.
