This paper considers a continuous three‐phase polynomial regression model with two threshold points for dependent data with heteroscedasticity. We assume the model is polynomial of order zero in the middle regime, and is polynomial of higher orders elsewhere. We denote this model by ℳ2$$ {\mathcal{M}}_2 $$, which includes models with one or no threshold points, denoted by ℳ1$$ {\mathcal{M}}_1 $$ and ℳ0$$ {\mathcal{M}}_0 $$, respectively, as special cases. We provide an ordered iterative least squares (OiLS) method when estimating ℳ2$$ {\mathcal{M}}_2 $$ and establish the consistency of the OiLS estimators under mild conditions. When the underlying model is ℳ1$$ {\mathcal{M}}_1 $$ and is false(d0prefix−1false)$$ \left({d}_0-1\right) $$th‐order differentiable but not d0$$ {d}_0 $$th‐order differentiable at the threshold point, we further show the Opfalse(Nprefix−1false/false(d0+2false)false)$$ {O}_p\left({N}^{-1/\left({d}_0+2\right)}\right) $$ convergence rate of the OiLS estimators, which can be faster than the Opfalse(Nprefix−1false/false(2d0false)false)$$ {O}_p\left({N}^{-1/\left(2{d}_0\right)}\right) $$ convergence rate given in Feder when d0≥3$$ {d}_0\ge 3 $$. We also apply a model‐selection procedure for selecting ℳκ$$ {\mathcal{M}}_{\kappa } $$; κ=0,1,2$$ \kappa =0,1,2 $$. When the underlying model exists, we establish the selection consistency under the aforementioned conditions. Finally, we conduct simulation experiments to demonstrate the finite‐sample performance of our asymptotic results.