Why does examining fitted vs residuals plot help us determine whether there is heteroskedasticity or not? Could someone give me a detailed theoretical rationale for this test? Does the randomness of residuals around the horizontal line y = 0 guarantee homoskedasticity, or does it just support the presence of it?
In my opinion "examining fitted vs residuals" is a graphical alternative to the Box-Cox test When (and why) should you take the log of a distribution (of numbers)? which can help to evaluate the hypothesis that the variability of residuals is or is not linearly related (sympathetic) to the level i.e.the fitted values of the series.
IFF you are analyzing time series or spatial data there is an additional scheme available....
"examining the plot of the residuals over time" is a graphical way of assessing possible deterministic change points in the error variance which suggest periods of time where the variance of the errors changes significantly. See a formal test for this here http://docplayer.net/12080848-Outliers-level-shifts-and-variance-changes-in-time-series.html
Finally if you don't treat one-time anomalies as just that you can often mis-analyze the data as was done with the classic INTERNATIONAL AIRLINE PASSENGER DATA example presented by Box and Jenkins in 1964 https://autobox.com/cms/index.php/blog . A few anomalies at the highest level of the series skewed the Box-COX test to conclude the higher levels meant higher model error variance which is not true in general. https://autobox.com/pdfs/vegas_ibf_09a.pdf has some very illuminating and specific details on the AIRLINE SERIES data.
Box & Jenkins conclusion (without suggesting viable alternatives) that you need to take logarithms was roundly criticized in the respected literature by other time series statisticians particularly Professor Chris Chatfield from Bath University.
Heteroskedasticity means quite general that the residual variance is not homogenous, i.e. depends on groups or some other variable. Note that, according to this definition,
if you find that var(res) ~ fitted, you clearly have heteroskedasticity
However, if you don't find a relationship between var(res) ~ fitted, you could still have heteroskedasticity, for example when testing var(res) ~ predictor1. In principle, you have to test against all possible predictors to exclude the possibility of heteroskedasticity.
In general, it is therefore advisable to plot residuals against all predictors and grouping variables (after all, when you model heteroskedasticity, you can make it also dependent on particular variables). This, however, means that you have to do many plots in a multiple regression. If you want to do only one plot var(res) ~ fitted is a good compromise, because fitted will be correlated with all predictors that have an effect, and are therefore of particular interest in an analysis of heteroskedasticity.