\chapter{Degrees of freedom correction}
\label{chap:df}


\section{Introduction}

This chapter gives a brief account of the issue of correction for
degrees of freedom in the context of econometric modeling, leading up
to a discussion of the policies adopted in gretl in this regard.  We
also explain how to supplement the results produced automatically by
gretl if you want to apply such a correction where gretl does not, or
vice versa.

The first few sections are quite basic; experts are invited to skip to
section~\ref{sec:df-grey}.

\section{Back to basics}
\label{sec:df-basics}

It's well known that given a sample, $\{x_i\}$, of size $n$ from a
normally distributed population, the Maximum Likelihood (ML) estimator
of the population variance, $\sigma^2$, is
%
\begin{equation}
\label{eq:sigma-hat}
\hat{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \bar{x})^2
\end{equation}
%
where $\bar{x}$ is the sample mean, $n^{-1} \sum_{i=1}^n x_i$.  It's
also well known that $\hat{\sigma}^2$, while it is a consistent
estimator, is biased, and it is commonly replaced by the ``sample
variance'', namely,
%
\begin{equation}
\label{eq:sample-variance}
s^2 = \frac{1}{n-1} \sum_{i=1}^n (x_i - \bar{x})^2
\end{equation}

The intuition behind the bias in (\ref{eq:sigma-hat}) is
straightforward.  First, the quantity we seek to estimate is defined
as
%
\[
\sigma^2 = E\left[(x_i - \mu)^2\right]
\]
%
where $\mu = E(x)$. It is clear that if $\mu$ were observable, a
perfectly good estimator would be
%
\[
\tilde{\sigma}^2 = \frac{1}{n} \sum_{i=1}^n (x_i - \mu)^2.
\]
%
But this is not a practical option: $\mu$ is generally
unobservable. We therefore substitute $\bar{x}$ for the unknown
$\mu$. It is easily shown that $\bar{x}$ is the least-squares
estimator of $\mu$, and also (assuming normality) the ML estimator.
It is unbiased, but is of course subject to sampling error; in any
given sample it is highly unlikely that $\bar{x} = \mu$.  Given that
$\bar{x}$ is the least-squares estimator, the sum of squared
deviations of the $x_i$ from \textit{any} value other than $\bar{x}$
must be greater than the summation in (\ref{eq:sigma-hat}).  But since
$\mu$ is almost certainly not equal to $\bar{x}$, the sum of squared
deviations of the $x_i$ from $\mu$ will surely be greater than the sum
of squared deviations in (\ref{eq:sigma-hat}). It follows that the
expected value of $\hat{\sigma}^2$ falls short of the population
variance.

The proof that $s^2$ is indeed the unbiased estimator can be found in
any good statistics textbook, where we also learn that the magnitude
$n-1$ in (\ref{eq:sample-variance}) can be brought under a general
description as the ``degrees of freedom'' of the calculation at
hand. (Given $\bar{x}$, the $n$ sample values provide only $n-1$ items
of information since the $n^{\rm th}$ value can always be deduced via the
formula for $\bar{x}$.)

\section{Application to OLS regression}
\label{sec:df-ols}

The argument above carries over into the usual calculation of standard
errors in the context of OLS regression as applied to the linear model
$y = X\beta + u$.  If the disturbances, $u$, are assumed to be
independently and identically distributed (IID), then the variance of
the OLS estimator, $\hat\beta$, is given by
%
\[
\mbox{Var}\left(\hat\beta\right) = \sigma^2 (X'X)^{-1}
\]
%
where $\sigma^2$ is the variance of the error term and $X$ is an
$n\times k$ matrix of regressors.  But how should the unknown
$\sigma^2$ be estimated?  The ML estimator is
%
\begin{equation}
\label{eq:ols-sigma2}
\hat\sigma^2 = \frac{1}{n} \sum_{i=1}^n \hat{u}^2_i
\end{equation}
%
where the $\hat{u}^2_i$ are squared residuals, $u_i = y_i - X_i\beta$.
But this estimator is biased and we typically use the unbiased
counterpart
%
\begin{equation}
\label{eq:ols-s2}
s^2 = \frac{1}{n-k} \sum_{i=1}^n \hat{u}_i^2
\end{equation}
%
in which $n - k$ is the number of degrees of freedom given $n$
residuals from a regression where $k$ parameters are estimated.

The standard estimator of the variance of $\hat\beta$ in the context
of OLS is then $V = s^2 (X'X)^{-1}$.  And the standard errors of the
individual parameter estimates, $s_{\hat{\beta}_i}$, being the square
roots of the diagonal elements of $V$, inherit a degrees of freedom
correction from the estimator $s^2$.

Going one step further, consider hypothesis testing in the context of
OLS.  Since the variance of $\hat\beta$ is unknown and must itself be
estimated, the sampling distribution of the OLS coefficients is not,
strictly speaking, normal.  But if the \textit{disturbances} are
normally distributed (besides being IID) then, even in small samples,
the parameter estimates will follow a distribution that can be
specified exactly, namely the Student $t$ distribution with degrees
of freedom equal to the value given above, $\nu = n-k$.

That is, besides using a df correction in computing the standard
errors of the OLS coefficients, one uses the same $\nu$ in selecting
the particular distribution to which the ``$t$-ratio'',
$(\hat{\beta}_i-\beta^0)/s_{\hat{\beta}_i}$, should be referred in
order to determine the marginal significance level or $p$-value for
the null hypothesis that $\beta_i = \beta^0$.  This is the payoff to
df correction: we get test statistics that follow a known distribution
in small samples.  (In big enough samples the point is moot, since the
quantitative distinction between $\hat\sigma^2$ and $s^2$ vanishes.)

So far, so good.  Everyone expects df correction in plain OLS standard
errors just as we expect division by $n-1$ in the sample variance.
And users of econometric software expect that the $p$-values reported
for OLS coefficients will be based on the $t(\nu)$
distribution---although they are not always sufficiently aware that
the validity of such statistics in small samples depends on the
assumption of normally distributed errors.

\section{Beyond OLS}
\label{sec:df-beyond}

The situation is different when we move beyond estimation of the
classical linear model via OLS.  We may wish to estimate nonlinear
models (sometimes by least squares), and many models of interest to
econometricians are commonly estimated via maximization of a
likelihood function, or via the generalized method of moments (GMM).

In such cases we do not, in general, have exact small-sample results
to rely upon; in particular, we cannot assume that coefficient
estimates follow the $t$ distribution.  Rather, we typically appeal to
asymptotic results in statistical theory.  We seek \textit{consistent}
estimators which, although they may be biased, nonetheless converge in
probability to the corresponding parameter values as the sample size
goes to infinity.  Under the right conditions, laws of large numbers
and central limit theorems entitle us to expect that test statistics
will converge to the normal distribution, or the $\chi^2$ distribution
for multivariate tests, given big enough samples.

\subsection{To ``correct'' or not?}

The question arises, should we or should we not apply a df
``correction'' in reporting variance estimates and standard errors
for models that depart from the classical linear specification?

The argument against applying df adjustment is that it lacks a
theoretical basis: it does not produce test statistics that follow any
known distribution in small samples.  In addition, if parameter
estimates are obtained via ML, it makes sense to report ML estimates
of variances even if these are biased, since it is the ML quantities
that are used in computing the criterion function and in forming
likelihood-ratio tests.

On the other hand, pragmatic arguments for doing df adjustment are (a)
that it makes for closer comparability between regular OLS estimates
and nonlinear ones, and (b) that it provides a ``pinch of salt'' in
relation to small-sample results---that is, it inflates standard
errors, confidence intervals and $p$-values somewhat---even if it
lacks rigorous justification.

Note that even for fairly small samples, the difference between the
biased and unbiased estimators in equations (\ref{eq:sigma-hat}) and
(\ref{eq:sample-variance}) above will be small.  For example, if
$n=30$ then $s^2 = \frac{30}{29} \hat\sigma^2$.  In econometric
modelling proper, however, the difference can be quite substantial.
If $n=50$ and $k=10$, the $s^2$ defined in (\ref{eq:ols-s2}) will be
$50/40 = 1.25$ as large as the $\hat\sigma^2$ in
(\ref{eq:ols-sigma2}), and standard errors will be about 12 percent
larger.\footnote{A fairly typical situation in time-series
  macroeconometrics would be have between 100 and 200 quarterly
  observations, and to be estimating up to maybe 30 parameters
  including lags.  In this case df correction would make a difference
  to standard errors on the order of 10 percent.}  One can make a case
for inflating the standard errors obtained via nonlinear estimators as
a precaution against taking results to be ``more precise than they
really are''.

In rejoinder to the last point, one might equally say that savvy
econometricians should know to apply a discount factor (albeit an
imprecise one) to small-sample estimates outside of the classical,
normal linear model---or even that they should distrust such results
and insist on large samples before making inferences. This line of
thinking suggests that test statistics such as
$z = \hat{\beta}_i/\hat\sigma_{\hat{\beta}_i}$ should be referred to
the distribution to which they conform asymptotically---in this case
$N(0,1)$ for $H_0: \beta_i = 0$---if and only if the conditions for
appealing to asymptotic results can be considered as met.  From this
point of view df adjustment may be seen as providing a false sense of
security.

\section{Consistency and awkward cases}
\label{sec:df-grey}

Consistency (in the ordinary sense of uniformity of treatment) is a
bugbear when dealing with this issue.  To give a simple example,
suppose an econometrics program follows the policy of applying df
correction for OLS estimation but not for ML estimation.  One is, of
course, free to estimate the classical, normal linear model via ML, in
which case $\hat\beta$ should be numerically identical to that
obtained via OLS.  But the user of the software will obtain two
different sets of standard errors depending on the estimation method.
Admittedly, this example is not very troublesome; presumably one would
apply ML to the classical linear model only to make a pedagogical
point.

Here is a more awkward case.  An unrestricted vector autoregression
(VAR) is a system of equations, but the ML estimate of this system,
given normal errors, is equivalent to equation-by-equation OLS.
Should df correction be applied to VARs? Consistency with OLS argues
Yes. However, a popular extension of the VAR methodology is the vector
error-correction model (VECM).  VECMs are closely related to VARs and
one might well be interested in making comparisons across the two, but
a VECM is a nonlinear system and the cointegrating vectors that lie at
the heart of this model must be estimated via Maximum Likelihood.  So
perhaps VAR results should \textit{not} be df adjusted, for
comparability with VECMs.

Another ``grey area'' is the class of Feasible Generalized Least
Squares (FGLS) estimators---for example, weighted least squares
following the estimation of a skedastic function, or estimators
designed to handle first-order autocorrelation, such as
Cochrane--Orcutt.  These depart from the classical linear model, and
the theoreretical basis for inference in such models is asymptotic,
yet according to econometric tradition standard errors are generally
df adjusted.

Yet another awkward case: ``robust'' (heteroskedasticity- and/or
autocorrelation-consistent) standard errors in the context of OLS.
Such estimators are justified by asymptotic arguments and in general
we cannot determine their small-sample distributions. That would argue
for referring the associated test statistics to the normal
distribution. But comparability with classical standard errors pulls
in the other direction. Suppose in a particular case a robust
estimator produces a standard error that is numerically
indistinguishable from the classical one: if the former is referred to
the normal distribution and the latter to the $t$ distribution,
switching to robust standard errors will give a smaller $p$-value for
the coefficient in question, making it appear ``more significant,''
and arguably this is misleading.


\section{What gretl does}
\label{sec:df-policy}

First of all, the third column in gretl model output---following
``coefficient'' and ``std.\ error''---is labeled either ``t-ratio'' or
``z.'' This is your signal: ``t-ratio'' indicates that the estimated
standard error employs a degrees of freedom adjustment and the
reported $p$-value is obtained from the Student $t$ distribution,
while ``z'' indicates that no such adjustment is applied and the
$p$-value comes from the standard normal distribution.

If you see that gretl is applying a df adjustment but you don't want
this, the first point to check is whether you can switch to the
asymptotic variant by using an option flag or other command.
\begin{itemize}
\item The \texttt{ols} and \texttt{tsls} commands support a
  \option{no-df-corr} option to suppress degrees of freedom
  adjustment. In the case of Two-Stage Least Squares it's certainly
  arguable that df correction should not be performed by default,
  however gretl does this, largely for comparability with other
  software (for example \textsf{Stata}'s \texttt{ivreg} command).  But
  you can override the default if you wish.
\item The \texttt{estimate} command, for systems of equations, also
  supports the \option{no-df-corr} option when the specified
  estimation method is OLS or TSLS. (For other estimators supported by
  gretl's \texttt{system} command no df adjustment is applied by
  default.)
\item By default gretl uses the $t$ distribution for statistics based
  on robust standard errors under OLS. However, users can specify that
  $p$-values be calculated using the standard normal distribution
  whenever the \option{robust} option is passed to an estimation
  command, by means of the following ``set'' command
%
\begin{code}
set robust_z on
\end{code}
\end{itemize}

If these possibilities do not apply, it is fairly straightforward to
``purge'' regression results of df correction, as illustrated in the
following script fragment. We assume that a model has just been
estimated, so that the model-related accessors
(\verb|$stderr|, \verb|$coeff| and so on) are available.
%
\begin{code}
matrix se = $stderr * sqrt($df/$T)
matrix zscore = $coeff ./ se
matrix pv = 2 * pvalue(z, abs(zscore))
matrix M = $coeff ~ se ~ zscore ~ pv
cnameset(M, "coeff stderr z p-value")
print M
\end{code}

This will print the original coefficient estimates along with
asymptotic standard errors and the associated
$z$-scores and (two-sided) normal $p$-values. The converse case is
left as an exercise for the reader.

\subsection{VARs}

As mentioned above, Vector Autoregressions constitute a particularly
awkward case, with considerations of consistency of treatment pulling
in two opposite directions. For that reason gretl has adopted an
``agnostic'' policy in relation to such systems. We do not offer a
\dollar{vcv} accessor, but instead accessors named \dollar{xtxinv}
(the matrix $X'X^{-1}$ for the system as a whole) and \dollar{sigma}
(an estimate of the cross-equation variance--covariance matrix,
$\Sigma$). It's then up to the user to build an estimate of the
variance matrix of the parameter estimates---call it $V$---should that
be required.

Note that \dollar{sigma} gives the Maximum Likelihood Estimator
(without a degrees of freedom adjustment) so if you do
%
\begin{code}
matrix Vml = $sigma ** $xtxinv
\end{code}
%
(where ``\texttt{**}'' represents Kronecker product) you obtain the
MLE of the variance matrix of the parameter estimates. But if you want
the unbiased estimator you can do
%
\begin{code}
matrix S = $sigma * $T/($T-$ncoeff)
matrix Vu = S ** $xtxinv
\end{code}
%
to employ a suitably inflated variant of the
$\Sigma$ estimate. (For VARs, and also VECMs, \dollar{ncoeff} gives
the number of coefficients per equation.)

The second variant above is such that the vector of standard errors
produced by
%
\begin{code}
matrix SE = sqrt(diag(Vu))
\end{code}
%
agrees with the standard errors printed as part of the per-equation
VAR output.

A fuller example of usage of the \dollar{xtxinv} accessor is given in
Listing~\ref{GCtest}: this shows how one can replicate the $F$-tests
for Granger causality that are displayed by default by the \cmd{var}
command, with the refinement that, depending on the setting of the
\verb|USE_F| flag, these tests can be done using a small sample
correction as in gretl's output or in asymptotic ($\chi^2$) form.

\begin{script}[htbp]
  \scriptcaption{Computing statistics to test for Granger causality}
  \label{GCtest}
\begin{scode}
open denmark.gdt
list LST = LRM LRY IBO IDE
scalar p = 2 # lags in VAR
scalar USE_F = 1 # small sample correction?

var p LST --quiet

k = nelem(LST)
matrix theta = vec($coeff)
matrix V = $sigma ** $xtxinv
if USE_F
  scalar df = $T - $ncoeff
  V *= $T/df
endif

matrix GC = zeros(k, k)
cnameset(GC, LST)
rnameset(GC, LST)

matrix idx = seq(1,p) + 1
loop i = 1..k
  loop j = 1..k
    GC[i,j] = qform(theta[idx]', invpd(V[idx,idx]))
    idx += (j==k)? p+1 : p
  endloop
endloop

if USE_F
  GC /= p
  matrix pvals = pvalue(F, p, df, GC)
else
  matrix pvals = pvalue(X, p, GC)
endif

cnameset(pvals, LST)
rnameset(pvals, LST)
print GC pvals
\end{scode}
\end{script}

Vector Error Correction Models are more complex than VARs in this
respect, since we employ Johansen's variance estimator for the
``$\beta$'' terms. 
This means for example that the \dollar{xtxinv} accessor treats each estimated 
error correction (EC) term as one regressor on its own, such that the sampling uncertainty
of the loading coefficients is thereby addressed (after Kronecker-multiplying with 
\dollar{sigma} as before) . The ``internals'' of the EC terms
are of course made up of the integrated (levels) variables, and the special
\dollar{jvbeta} accessor is responsible for the variance of the cointegration coefficients,
where degrees-of-freedom corrections are not available. 

But as soon as the loading coefficients attached to the EC terms are restricted, there is 
no common set of regressors with freely varying coefficients in the VECM system anymore,
and therefore in these cases the formulas above are misleading. The \dollar{xtxinv} 
accessor can still be retrieved (because it does not involve the coefficients), but in 
the restricted $\alpha$ case it should no longer be used as shown above. 
The notion of system degrees of freedom then also becomes fuzzier since the number of 
regressors can vary across equations.

%%% Local Variables:
%%% mode: latex
%%% TeX-master: "gretl-guide"
%%% End:
