How to calculate a phase shift in the process of deriving Baade-Wesselinke distance? The above is the equation used in the Baade-Wesselink method to independently get distances to Cepheid variables (in the process of deriving Period-Luminosity relation).

According to the paper New Baade-Wesselink distances and radii for four metal-rich Galactic Cepheids, δθ in the equation accounts for a phase shift between the radial velocity curve and the angular diameter changes measured via surface-brightness relation. My question is how can I calculate this δθ value given that I have lists of radial velocity values and V-R values with the observed HJD dates (which is possible to be converted in terms of phase).
It is admirable that you are tackling this type of problem as a high school student, but I really think that some of these topics are out of your grasp currently, partly because of the highly technical nature, but also because of your lack of access to the necessary data. 

For one thing, you say you have V-R colors as a function of time, but all these equations were calibrated to use V-K colors. From your linked paper, you're going to need equation 2 for your process to calculate $p$ and equation 3 for getting $\theta_0$. Unfortunately, these equations assume V-K colors. You can redo the analysis done in this paper to calculate $p$ and this paper to get $\theta_0$ for V-R colors, but it's not going to be easy.

For another, you're going to need measurements of the angular size as a function of time (and your cited paper suggests that only exists for a handful of Cepheids and was measured using the VLTI) as well as the RV as a function of time. The $\delta \theta$ that you reference is described exactly in the paper as "a  phase shift  between  the  radial  velocity  curve  and  the  angular  diameter changes measured either interferometrically or via the SB relation." This means if you have $RV(t)$ and $\Delta \theta(t)$, you can determine a phase shift between the two curves, presumably by fitting some equation of the form $f(t) = a\ \mathrm{sin}(bt+c)$ to both equations and saying that $\delta \theta = c_{RV} - c_{\Delta \theta}$.

The analysis you're trying to do requires very specific data and it seems that the relevant papers on this subject have already done the analysis on all stars in which all necessary data existed. You cannot easily apply this process to any Cepheid without significant work.