Q:

Is it important to model heteroscedasticity during multiple regression?

Given a multiple linear regression (eg. using a GLS procedure) between a response variable and several predictive variables with different, heteroscedastic relationships with the response variable and among them. How important it is to introduce in the model estimations of the different heteroscedasticities among model variables? e.g. between response and predictor 1, or among predictors one and 2...

A:

If I understand your question correctly, then you are asking how important it is to take heteroskedascticity into account when modelling in a multivariate linear regression setting.
I will use the following model for reference:
$$ y = X \beta + \varepsilon $$
where $y$ is $(n \times 1)$, $X$ is $(n \times k)$, and $\beta$ is $(k \times 1)$. So you have $k$ parameters including the constant and $n$ observations and $\varepsilon$ is the residual with same dimensions as $y$.
This representation is equivalent to:
$$y_i = \beta_0 + \beta_1 x_{i,1} + ... + \beta_{k-1} x_{i,k-1} + \varepsilon_i$$
So for ordinary least squares to be the best linear unbiased estimator (BLUE) you need the Gauss-Markov (GM) assumptions to hold:

$\mathit{E}(\varepsilon_i) = 0$ i.e. residuals have zero mean
$\mathit{Var}(\varepsilon_i) = \sigma^2 < \infty$ i.e. constant bounded residual variance
$\mathit{Cov}(\varepsilon_i,\varepsilon_j) = 0$ for $i \neq j$ i.e. homoskedascticity

You also need:

$X'X$ is full rank i.e. no multicollinearity

In the case conditions 2 or/and 3 are invalid OLS is still a consistent estimator, and furthermore your estimate for $\beta$ will be unbiased, thus valid. However you will not be able to test the significance of $\beta$ (i.e. t-test) as your standard errors will be incorrect. If this is important for you, then you need to either use GLS/FGLS (which will give a different estimate for $\beta$) or use a post-correction model e.g. HAC (heteroskedascticity-consistent standard errors).
GLS is a simple extension to OLS:
The estimate for OLS is:
$$\hat{\beta} = (X'X)^{-1}X'y$$
For GLS:
$$\hat{\beta} = (X'\Omega^{-1}X)^{-1}X'\Omega^{-1}y$$
where $\Omega$ characterises your heteroskedascticity structure i.e. $\mathit{E}(\varepsilon\varepsilon') = \Omega$
Typically, we do not know $\Omega$ so we must estimate it using an FGLS (Feasable GLS) procedure which is a multi-step process:

Estimate $\varepsilon$ using OLS
$$ \hat{\varepsilon} = y - X\hat{\beta}_{OLS}$$
Calculate $\hat{\Omega}_{OLS}$ i.e. $\hat{\Omega}_{OLS} = \frac{\hat{\varepsilon}\hat{\varepsilon}'}{n-1}$
Use $\hat{\Omega}_{OLS}$ to calculate $\hat{\beta}_{FGLS1}$
Repeat step 1 and 2 using $\hat{\beta}_{FGLS1}$ instead of $\hat{\beta}_{OLS}$ to obtain $\hat{\Omega}_{FGLS1}$
Finally use $\hat{\Omega}_{FGLS1}$ with the GLS procedure to calculate a final estimate for $\beta$

To answer your question, importance depends on whether you want to test significance of parameters or not and how big your sample size is. With large samples ($n > 200$) the efficiency gains between OLS and GLS will start to become marginal.

