---
title: "Feasible GLS under multiplicative heteroskedasticity"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Feasible GLS under multiplicative heteroskedasticity}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.5
)
options(hcinfer.use_emoji = FALSE)
```

Feasible generalized least squares (FGLS) complements the
heteroskedasticity-consistent (HC) estimators in `hcinfer`. The HC estimators
keep the ordinary least squares (OLS) coefficients and only robustify their
covariance, assuming nothing about the form of the variance function. Feasible
GLS instead models the conditional variance as an exponential function of
observed regressors, re-weights the regression accordingly, and can be more
efficient than OLS when that variance model is adequate. The `gls_mult()`
function offers two estimators of this model, Harvey's two-step procedure and
Gaussian maximum likelihood, and the maximum likelihood fit supplies a proper
`logLik()` together with `AIC()` and `BIC()`. The method follows Harvey (1976)
and Cribari-Neto and Pereira (2019).

## HC inference versus feasible GLS

The two approaches answer different questions and suit different situations.
Calling `hcinfer()` or `vcov_hc()` returns the OLS point estimates together with
a heteroskedasticity-consistent sandwich covariance that is valid under
heteroskedasticity of unknown form, so it makes no commitment to how the
variance depends on the regressors. Calling `gls_mult()` instead estimates an
explicit variance model, returns re-weighted coefficients that can be more
efficient than OLS, and reports a model-based covariance that is valid when the
variance model is correctly specified. As a practical rule, the HC estimators
are the natural choice when you do not wish to model the variance and want
inference that is agnostic about its form, whereas feasible GLS is attractive
when a multiplicative variance model is plausible and the efficiency gain
matters.

## The multiplicative variance model

Feasible GLS in `hcinfer` is built on a linear mean model paired with a
multiplicative model for the conditional variance. The mean model is

$$
y = X \beta + e,
\qquad E(e_t) = 0,
\qquad \operatorname{Var}(e_t) = \sigma_t^2,
$$

where $y$ is the $n \times 1$ response, $X$ is the $n \times p$ mean design
matrix, and $\beta$ is the vector of mean coefficients. The variance is modelled
as an exponential function of a second set of regressors,

$$
\sigma_t^2 = \exp(\eta_t),
\qquad \eta_t = z_t^\top \gamma,
$$

so that the log-variance $\eta_t$ is linear in the dispersion regressors $z_t$,
the rows of a dispersion design matrix $Z$, with coefficients $\gamma$. Because
the response has response units, $\exp(\eta_t)$ has squared response units and
$\exp(\eta_t/2)$ is the conditional standard deviation in response units. The
exponential link keeps the fitted variances positive for any value of $\gamma$.
The `variance` argument selects $Z$: a one-sided formula such as
`~ income_scaled` builds $Z$ from the named terms, while the default
`variance = NULL` sets $Z = X$ and reuses the mean design. This is a model for
the variance, not an assumption that the residuals or leverage complements
follow any particular distribution.

## Two estimators

Both estimators target the same mean and variance model but reach it by
different routes. The two-step estimator is a single pass through an auxiliary
regression, whereas maximum likelihood locally optimises the Gaussian profile
likelihood and is accepted only at a stationary solution. Maximum likelihood is
the default, and it is also the fit that supports the information criteria.

### Harvey two-step

The two-step estimator of Harvey (1976) begins from the OLS residuals
$\hat e_t$ and regresses $\log(\hat e_t^2)$ on the dispersion regressors $Z$.
Because the expectation of $\log \varepsilon^2$ under normality is not zero, the
intercept of this auxiliary regression is bias-corrected by adding
$c = -\operatorname{digamma}(1/2) - \log 2 \approx 1.2704$, the bias of
$\log \varepsilon^2$ under normality. The fitted log-variances yield the weights
$\hat w_t = \exp(-z_t^\top \hat\gamma)$, and a single weighted least squares pass
with these weights delivers the mean coefficients, with reported dispersion
covariance $\tfrac{\pi^2}{2}(Z^\top Z)^{-1}$.

### Maximum likelihood

Maximum likelihood optimises the Gaussian log-likelihood of the same model. The
mean coefficients $\beta$ are profiled out by weighted least squares at each
candidate $\gamma$, so a BFGS optimiser searches only over the dispersion
parameters $\gamma$, and the reported asymptotic dispersion covariance is
$2(Z^\top Z)^{-1}$. Under the model and standard regularity conditions, maximum
likelihood is asymptotically efficient. For an accepted locally optimised
stationary fit, `logLik()` returns the likelihood evaluated at that solution,
which enables `AIC()` and `BIC()` without claiming that the stationarity check
has certified a global maximum.

## Data and model

The examples use the `PublicSchools2` data, which have complete observations for
the 50 U.S. states and the District of Columbia and include a `south` indicator
that the original `PublicSchools` data lack. The `expenditure` response is
annual expenditure per student enrolled in K-12 public schools for 2025,
measured in U.S. dollars. Income is divided by ten thousand so that the
estimated coefficients are on a readable scale. The mean model regresses
expenditure on rescaled income and the regional indicator, and printing the
fitted `lm` object shows the OLS starting point that `gls_mult()` will refine.

```{r}
library(hcinfer)

schools <- PublicSchools2
schools$income_scaled <- schools$income / 10000

fit <- lm(expenditure ~ income_scaled + south, data = schools)
fit
```

## Maximum likelihood fit

With no `estimator` argument, `gls_mult()` fits the model by maximum
likelihood, its default estimator, and optimizes the profile likelihood with
the algorithm chosen by `method` (BFGS by default). The printed object gives a
compact view of the mean and dispersion coefficients together with the
convergence diagnostics.

```{r}
gls_fit <- gls_mult(fit)
gls_fit
```

```{r}
summary(gls_fit)
```

The summary reports fitted conditional variances in squared U.S. dollars, or
USD squared, because `expenditure` is measured in USD and $\exp(\eta_t)$ is a
variance. On the standard-deviation scale, the corresponding values range from
about USD `r formatC(min(sqrt(gls_fit$fitted_variances)), digits = 0, format = "f", big.mark = ",")`
to USD `r formatC(max(sqrt(gls_fit$fitted_variances)), digits = 0, format = "f", big.mark = ",")`
and are comparable to the OLS residual standard deviation of about USD
`r formatC(stats::sigma(fit), digits = 0, format = "f", big.mark = ",")`. The
million-scale variance entries are therefore expected for this response scale
and are not evidence of an exponentiation error. The same report gives the
mean and dispersion blocks, normal Wald tests, and the likelihood evaluated at
the accepted local solution together with AIC and BIC.

## Extracting the fit

The usual extractor functions work on the fitted object and, where relevant,
take a `model` argument to choose between the mean and dispersion blocks.

```{r}
coef(gls_fit)                        # mean coefficients
coef(gls_fit, model = "dispersion")  # log-variance coefficients
```

```{r}
vcov(gls_fit)                        # mean covariance
vcov(gls_fit, model = "dispersion")  # dispersion covariance
```

```{r}
tests(gls_fit)
confint(gls_fit)
```

```{r}
head(fitted(gls_fit))
head(gls_fit$fitted_variances)   # estimated conditional variances
head(gls_fit$weights)            # GLS weights exp(-eta)
```

By default `coef()` and `vcov()` return the mean block, while
`model = "dispersion"` selects the log-variance coefficients and their
covariance. The `tests()` and `confint()` methods report normal Wald inference
on the mean coefficients, using the model-based standard errors. The fitted
object also stores the conditional variances in `fitted_variances` and the GLS
weights $\exp(-\eta_t)$ in `weights`, so the estimated variance structure is
available for diagnostics or plotting. These weights have
inverse-squared-response units, and multiplying all of them by one positive
constant leaves the weighted least squares coefficients unchanged, so only
their relative scale affects those coefficients.

## Harvey two-step

Setting `estimator = "two_step"` fits the same model with Harvey's non-iterative
estimator.

```{r}
two_step <- gls_mult(fit, estimator = "two_step")
coef(two_step)
```

```{r}
data.frame(
  term = names(coef(gls_fit)),
  ml = coef(gls_fit),
  two_step = coef(two_step)
)
```

The two-step estimator uses one auxiliary regression, whereas maximum
likelihood locally optimises the profile likelihood and applies the stationarity
guard. In this example, the income slopes are similar, but the intercept and
`south` coefficients differ more noticeably, so the table demonstrates
estimator dependence rather than equivalence. Information criteria are defined
only for the accepted maximum likelihood fit because the two-step point does
not optimise the likelihood and therefore carries no comparable `logLik()`.

## Custom variance models

The `variance` argument decouples the dispersion model from the mean model, so
the analyst can let the variance depend on a different, and usually smaller, set
of regressors.

```{r}
gls_income <- gls_mult(fit, variance = ~ income_scaled)
coef(gls_income, model = "dispersion")
```

Here the variance is modelled as a function of income alone, while the mean
model still includes both income and the regional indicator. The default
`variance = NULL` instead reuses the full mean design as $Z$, which is convenient
but not always the most parsimonious choice for the variance.

## Model selection with AIC and BIC

Because each maximum likelihood fit carries a proper `logLik()` with $p + q$
degrees of freedom, where $p$ counts the mean coefficients and $q$ the
dispersion coefficients, the base `AIC()` and `BIC()` generics compare competing
variance specifications directly. The three fits below share the same mean model
but differ in their variance model: the full dispersion model on income and
region, a reduced model on income alone, and the special case `variance = ~ 1`,
which forces a constant variance and so reduces the fit to the homoskedastic
Gaussian model.

```{r}
full <- gls_mult(fit)                       # variance ~ income + south
income_only <- gls_mult(fit, variance = ~ income_scaled)
homoskedastic <- gls_mult(fit, variance = ~ 1)
AIC(full, income_only, homoskedastic)
BIC(full, income_only, homoskedastic)
```

The specification with the smallest criterion is preferred, so these tables let
the data adjudicate between richer and sparser variance models on the same
footing. Because they rest on the likelihood, the criteria are defined only for
maximum likelihood fits. Meaningful comparisons also require the competing fits
to have reached comparable likelihood solutions.

## Comparison with HC standard errors

It is instructive to place the feasible GLS standard errors next to the OLS and
HC standard errors for the same model.

```{r}
data.frame(
  term = names(coef(fit)),
  ols = sqrt(diag(vcov(fit))),
  fgls_ml = sqrt(diag(vcov(gls_fit))),
  hcbeta = sqrt(diag(vcov(hcinfer(fit, type = "hcbeta")))),
  hc3 = sqrt(diag(vcov(hcinfer(fit, type = "hc3"))))
)
```

The feasible GLS standard errors come from a re-weighted fit under an explicit
variance model, whereas the HC standard errors, here HC$\beta$ and HC3, keep the
OLS coefficients and only robustify their covariance without modelling the
variance. Agreement between the reported standard errors is descriptive and
does not by itself establish that the variance model is correctly specified.
Systematic differences can motivate inspection of the variance specification
and influential high-leverage observations, but they do not identify either
explanation on their own.

## Practical guidance

Several practical points help in routine use. Maximum likelihood is the
recommended default because it is efficient under the model and provides
information criteria, while the two-step estimator is a useful non-iterative
alternative. The mean covariance that `gls_mult()` reports is a model-based
plug-in, valid when the variance model is correctly specified rather than an HC
sandwich, so inference that should be agnostic about the variance form still
belongs to `hcinfer()`. Scaling the dispersion regressors sensibly can improve
numerical behaviour, because extreme scaling may require optimiser tuning
through `control`, such as a tighter `reltol`. The convergence code and
scale-invariant score check establish only a locally optimised stationary fit;
they do not prove that the accepted solution is a global maximum. For the
underlying theory of the HC estimators see
`vignette("hcinfer-methodology", package = "hcinfer")`, and for the full
argument list see `?gls_mult`.
