We single out the Starobinsky model and its extensions among generic $f(R)$ gravity as attractors at large field values for chaotic inflation. Treating a $R^3$ curvature term as a perturbation of the Starobinsky model, we impose the phenomenological bounds on the additional term satisfying the successful inflationary predictions... We find that the scalar spectral index can vary in both the red or blue tilted direction, depending on the sign of the coefficient of the $R^3$ term, whereas the tensor-to-scalar ratio is less affected in the Planck-compatible region. We also discuss the role of higher order curvature term for stability and the reheating dynamics for the unambiguous prediction for the number of efoldings up to the $R^3$ term. (read more)