14  Introduction

Statistical modeling applies statistical methods to real-world data to give empirical content to relationships. It aims to quantify phenomena and develop models and test hypotheses, making it a crucial field for economic research, policy analysis, and decision-making. The aim of the statistical modeling is to study the (unknown) mechanism that generates the data, i.e., the Data Generating Process (DGP). The statistical model is a function that approximates the DGP.

14.1 The matrix of data

Let’s consider $n$ realizations defining a sample for $i=1,2,\dots ,n$. Suppose we have $p$ dependent variables and $k$ explanatory variables (also known as predictors). The data matrix of the exogenous (regressors) ${\mathbf{X}}_{n,k}$ is defined as in Equation 13.1, while the matrix composed by the endogenous (dependent) variables $\mathbf{Y}$ reads $$\underset{n\times p}{\mathbf{Y}}=\left(\begin{array}{cccccc}{y}_{1,1}& y_{1,2}& \dots & y_{1,j}& \dots & y_{1,p}\\ y_{2,1}& y_{2,2}& \dots & y_{2,j}& \dots & y_{2,p}\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ y_{n,1}& y_{n,2}& \dots & y_{n,j}& \dots & y_{n,p}\end{array}\right)\text{.}$$ Hence, the complete matrix of data reads $$\begin{array}{}\text{(14.1)}& \underset{n\times (k+p)}{\mathbf{W}}=\left(\begin{array}{cc}\mathbf{Y}& \mathbf{X}\end{array}\right)=\left(\begin{array}{cccccc}{y}_{1,1}& \dots & y_{1,p}& x_{1,1}& \dots & x_{1,k}\\ ⋮& \ddots & ⋮& ⋮& \ddots & ⋮\\ y_{n,1}& \dots & y_{n,p}& x_{n,1}& \dots & x_{n,k}\end{array}\right)\text{.}\end{array}$$ In general, when $p=1$ then the model has only one equation to satisfy for $i=1,\dots ,n$, for example $$\begin{array}{}\text{(14.2)}& Y_i=b_0+b_1{x}_{i,1}+b_2{x}_{i,2}+\cdots +b_k{x}_{i,k}+u_i\text{.}\end{array}$$ Otherwise, when $p>1$ there are more than one dependent variable and the model is composed by $p$-equations for $i=1,\dots ,n$, i.e. the same linear model with $p$ equations reads: $$\begin{array}{}\text{(14.3)}& \{\begin{array}{l}{Y}_{i,1}=b_{0,1}+b_{1,1}{x}_{i,1}+b_{1,2}{x}_{i,2}+\cdots +b_{1,k}{x}_{i,k}+u_{i,1}\\ Y_{i,2}=b_{0,2}+b_{2,1}{x}_{i,1}+b_{2,2}{x}_{i,2}+\cdots +b_{2,k}{x}_{i,k}+u_{i,2}\\ ⋮\\ Y_{i,p}=b_{0,p}+b_{p,1}{x}_{i,1}+b_{p,2}{x}_{i,2}+\cdots +b_{p,k}{x}_{i,k}+u_{i,p}\end{array}\end{array}$$ Thus the matrix of the residuals components reads $$\underset{n\times p}{\mathbf{U}}=\left(\begin{array}{cccccc}{u}_{1,1}& u_{1,2}& \dots & u_{1,i}& \dots & u_{1,p}\\ u_{2,1}& u_{2,2}& \dots & u_{2,i}& \dots & u_{2,p}\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ u_{n,1}& u_{n,2}& \dots & u_{n,i}& \dots & u_{n,p}\end{array}\right)\text{,}$$

14.2 Joint, conditional and marginals

Let’s consider the bi-dimensional random vector in Equation 14.1 and let’s write the joint distribution of $\mathbf{X}$ and $\mathbf{Y}$, i.e.  $$\begin{array}{}\text{(14.4)}& \underset{\text{joint probability}}{\mathbb{P}({\mathbf{Y}}_1,\dots ,{\mathbf{Y}}_p\le {\mathbf{y}}_1,\dots ,{\mathbf{y}}_p,{\mathbf{X}}_1,\dots ,{\mathbf{X}}_p\le {\mathbf{x}}_1,\dots ,{\mathbf{x}}_k})=\underset{\text{distribution function}}{F_{\mathbf{Y},\mathbf{X}}({\mathbf{y}}_1,\dots ,{\mathbf{y}}_p,{\mathbf{x}}_1,\dots ,{\mathbf{x}}_k})\text{.}\end{array}$$ In the continuous case, there exists a joint density $f_{\mathbf{Y},\mathbf{X}}$ such that: $$\begin{array}{}\text{(14.5)}& F_{\mathbf{Y},\mathbf{X}}(\mathbf{y},\mathbf{x})={\int }_{-\mathrm{\infty }}^{\mathbf{x}}{\int }_{-\mathrm{\infty }}^{\mathbf{y}}{f}_{\mathbf{Y},\mathbf{X}}(\mathbf{y},\mathbf{x})d\mathbf{y}d\mathbf{x}\text{.}\end{array}$$ Moreover, from the joint distribution (Equation 14.4) it is possible to recover the marginals distributions, i.e.
$$\begin{array}{}\text{(14.6)}& \begin{array}{rl}& f_{\mathbf{Y}}(\mathbf{y})={\partial }_{\mathbf{y}}{F}_{\mathbf{Y},\mathbf{X}}(\mathbf{y},\mathbf{x})={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{f}_{\mathbf{Y},\mathbf{X}}(\mathbf{y},\mathbf{x})d\mathbf{x}\text{,}\\ & f_{\mathbf{X}}(\mathbf{x})={\partial }_{\mathbf{x}}{F}_{\mathbf{Y},\mathbf{X}}(\mathbf{y},\mathbf{x})={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{f}_{\mathbf{Y},\mathbf{X}}(\mathbf{y},\mathbf{x})d\mathbf{y}\text{.}\end{array}\end{array}$$

Then, given the marginals (Equation 14.6), it is possible to compute the unconditional moments, for example

• First moment: $\mathbb{E}\{\mathbf{Y}\}={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathbf{y}{f}_{\mathbf{Y}}(\mathbf{y})d\mathbf{y}$.
• Second moment: $\mathbb{E}\{{\mathbf{Y}}^2\}={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathbf{y}}^2{f}_{\mathbf{Y}}(\mathbf{y})d\mathbf{y}$.

Applying the Bayes theorem (Theorem 6.2), from the joint distribution (Equation 14.4) it is possible to recover the conditional distribution, i.e
$$\begin{array}{}\text{(14.7)}& f_{\mathbf{Y}\mid \mathbf{X}}(\mathbf{y}\mid \mathbf{x})=\frac{f_{\mathbf{Y},\mathbf{X}}(\mathbf{y},\mathbf{x})}{f_{\mathbf{X}}(\mathbf{x})}⟺\underset{\text{joint}}{f_{\mathbf{Y},\mathbf{X}}(\mathbf{y},\mathbf{x})}=\underset{\text{conditional}}{f_{\mathbf{Y}\mid \mathbf{X}}(\mathbf{y}\mid \mathbf{x})}\cdot \underset{\text{marginal}}{f_{\mathbf{X}}(\mathbf{x})}\text{.}\end{array}$$ Given the conditional distribution the conditional moments reads $$\begin{array}{rl}& \mathbb{E}\{\mathbf{Y}\mid \mathbf{X}\}={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathbf{y}{f}_{\mathbf{Y}\mid \mathbf{X}}(\mathbf{y}\mid \mathbf{x})d\mathbf{y}\text{,}\\ & \mathbb{E}\{{\mathbf{Y}}^2\mid \mathbf{X}\}={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\mathbf{y}}^2{f}_{\mathbf{Y}\mid \mathbf{X}}(\mathbf{y}\mid \mathbf{x})d\mathbf{y}\text{.}\end{array}$$

Example: Inference on a Multivariate Gaussian model

Example 14.1 Let’s consider a multivariate Gaussian setup, i.e.  $$\left(\begin{array}{c}\mathbf{Y}\\ \mathbf{X}\end{array}\right)\sim \text{MVN}(\left(\begin{array}{c}{\mathit{\mu }}_{\mathbf{Y}}\\ {\mathit{\mu }}_{\mathbf{X}}\end{array}\right),\left(\begin{array}{cc}{\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{Y}}& {\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{X}}\\ {\mathit{\mu }}_{\mathbf{X}\mathbf{Y}}& {\mathbf{\Sigma }}_{\mathbf{X}\mathbf{X}}\end{array}\right))\text{.}$$ If $(\mathbf{X}\mathbf{Y})$ are jointly normal, then the marginals are multivariate normal, i.e. $$\mathbf{Y}\sim \text{MVN}({\mathit{\mu }}_{\mathbf{Y}},{\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{Y}})\text{,}\mathbf{X}\sim \text{MVN}({\mathit{\mu }}_{\mathbf{X}},{\mathbf{\Sigma }}_{\mathbf{X}\mathbf{X}})\text{,}$$ and also the conditionals distributions, i.e.  $$\mathbf{Y}\mid \mathbf{X}\sim \text{MVN}({\mathit{\mu }}_{\mathbf{Y}\mid \mathbf{X}},{\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{Y}\mid \mathbf{X}})\text{,}\mathbf{X}\mid \mathbf{Y}\sim \text{MVN}({\mathit{\mu }}_{\mathbf{X}\mid \mathbf{Y}},{\mathbf{\Sigma }}_{\mathbf{X}\mathbf{X}\mid \mathbf{Y}})\text{.}$$ In such model’s setup the conditional expectation of $\mathbf{Y}$ given $\mathbf{X}$ reads $$\begin{array}{rl}\mathbb{E}\{\mathbf{Y}\mid \mathbf{X}\}& ={\mathit{\mu }}_{\mathbf{Y}\mid \mathbf{X}}=\\ & ={\mathit{\mu }}_{\mathbf{Y}}+{\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{X}}\cdot {\mathbf{\Sigma }}_{\mathbf{X}\mathbf{X}}^{-1}(\mathbf{X}-{\mathit{\mu }}_{\mathbf{X}})=\\ & ={\mathit{\mu }}_{\mathbf{Y}}-{\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{X}}\cdot {\mathbf{\Sigma }}_{\mathbf{X}\mathbf{X}}^{-1}{\mathit{\mu }}_{\mathbf{X}}+{\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{X}}\cdot {\mathbf{\Sigma }}_{\mathbf{X}\mathbf{X}}^{-1}\mathbf{X}=\\ & ={\mathit{\mu }}_{\mathbf{Y}}-{\mathbf{b}}_{\mathbf{Y}\mid \mathbf{X}}{\mathit{\mu }}_{\mathbf{X}}+{\mathbf{b}}_{\mathbf{Y}\mid \mathbf{X}}\mathbf{X}=\\ & ={\mathbf{a}}_{\mathbf{Y}\mid \mathbf{X}}+{\mathbf{b}}_{\mathbf{Y}\mid \mathbf{X}}\mathbf{X}\end{array}$$ and the conditional variance as $$\mathbb{V}\{\mathbf{Y}\mid \mathbf{X}\}={\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{Y}}-{\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{X}}\cdot {\mathbf{\Sigma }}_{\mathbf{X}\mathbf{X}}^{-1}{\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{X}}\text{.}$$ In this setup the parameters are:

• Joint distribution, $\mathit{\theta }=\{{\mathit{\mu }}_{\mathbf{Y}},{\mathit{\mu }}_{\mathbf{X}},{\mathbf{\Sigma }}_{\mathbf{X}\mathbf{X}},{\mathbf{\Sigma }}_{\mathbf{X}\mathbf{Y}},{\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{Y}}\}$.
• Conditional distribution, ${\mathit{\lambda }}_1=\{{\mathbf{a}}_{\mathbf{Y}\mid \mathbf{X}},{\mathbf{b}}_{\mathbf{Y}\mid \mathbf{X}},{\mathbf{\Sigma }}_{\mathbf{Y}\mathbf{Y}\mid \mathbf{X}}\}$.
• Marginal distribution, ${\mathit{\lambda }}_2=\{{\mathit{\mu }}_{\mathbf{X}},{\mathbf{\Sigma }}_{\mathbf{X}\mathbf{X}}\}$.

Noting that ${\mathit{\lambda }}_1$ is a function of $\mathit{\theta }$, i.e. $\tau =f({\mathit{\lambda }}_1)$ in the Gaussian case it is possible to prove that ${\mathit{\lambda }}_1$ and ${\mathit{\lambda }}_2$ are free to vary. Hence, imposing restrictions on ${\mathit{\lambda }}_1$ do not impose restrictions on ${\mathit{\lambda }}_2$. In general, if the parameters of interest are a function of the conditional distribution and ${\mathit{\lambda }}_1$ and ${\mathit{\lambda }}_2$ are free to vary, then the inference can be done without losing of information considering the conditional model. In this case we say that $\mathbf{X}$ is weakly exogenous for $\tau =f({\mathit{\lambda }}_1)$.

14.3 Conditional expectation model

Let’s consider a very general conditional expectation model with $p=1$, of which the linear models are a special case. In matrix notation it can be written as: $$\begin{array}{}\text{(14.8)}& \mathbf{y}=\mathbb{E}\{\mathbf{y}\mid \mathbf{X}\}+\mathbf{u}\text{,}\end{array}$$ where the conditional expectation errors are defined as: $$\begin{array}{}\text{(14.9)}& \mathbf{u}=\mathbf{y}-\mathbb{E}\{\mathbf{y}\mid \mathbf{X}\}\text{.}\end{array}$$

Proposition 14.1 In a conditional expectation model as in Equation 14.8, the residuals $\mathbf{u}$, defined as in Equation 14.9, have unconditional expectation and covariance with the regressors $\mathbf{X}$ equal to zero, i.e. $$\mathbb{E}\{\mathbf{u}\}=0,\mathbb{E}\{\mathbf{u}\mathbf{X}\}=0\text{.}$$ Moreover, the conditional expectation error is orthogonal to any transformation of the conditioning variables, i.e.  $$\begin{array}{}\text{(14.10)}& \mathbf{y}=\mathbb{E}\{\mathbf{y}\mid g(\mathbf{X})\}+\mathbf{u}⟹\mathbb{E}\{\mathbf{u}g(\mathbf{X})\}=0\text{.}\end{array}$$

Proof: Proposition 14.1

Proof. Let’s start the unconditional expectation of the residuals defined in Equation 14.8, i.e.  $$\begin{array}{rl}\mathbb{E}\{\mathbf{u}\}& =\mathbb{E}\{\mathbf{y}-\mathbb{E}\{\mathbf{y}\mid \mathbf{X}\}\}=\\ & =\mathbb{E}\{\mathbf{y}\}-\mathbb{E}\{\mathbb{E}\{\mathbf{y}\mid \mathbf{X}\}\}=\\ & =\mathbb{E}\{\mathbf{y}\}-\mathbb{E}\{\mathbf{y}\}=0\end{array}$$ Then, let’s compute the expected value of between the residuals and the regressors, i.e.  $$\mathbb{E}\{\mathbf{u}\}=0⟹\mathbb{C}v\{\mathbf{u},\mathbf{X}\}=\mathbb{E}\{\mathbf{u}\mathbf{X}\}\text{.}$$ For simplicity let’s assume that $\mathbf{X}$ can takes only values in $\{0,1\}$. Applying the tower property of conditional expectation one obtain: $$\begin{array}{rl}\mathbb{E}\{\mathbf{u}\mathbf{X}\}& =\mathbb{E}\{\mathbb{E}\{\mathbf{u}\mathbf{X}\mid \mathbf{X}\}\}=\\ & =\mathbb{E}\{\mathbf{u}\mathbf{X}\mid \mathbf{X}=0\}\mathbb{P}(\mathbf{X}=0)+\mathbb{E}\{\mathbf{u}\mathbf{X}\mid \mathbf{X}=1\}\mathbb{P}(\mathbf{X}=1)=\\ & =\mathbb{E}\{\mathbf{u}\mathbf{X}\mid \mathbf{X}=1\}\mathbb{P}(\mathbf{X}=1)\end{array}$$ Then, let’s substitute $\mathbf{u}$ from Equation 14.8 and $\mathbf{X}$ with $1$, i.e. $$\begin{array}{rl}\mathbb{E}\{\mathbf{u}\mathbf{X}\}& =\mathbb{E}\{(\mathbf{y}-\mathbb{E}\{\mathbf{y}\mid \mathbf{X}\})\mathbf{X}\mid \mathbf{X}=1\}\mathbb{P}(\mathbf{X}=1)=\\ & =\mathbb{E}\{\mathbf{y}\mid \mathbf{X}=1\}\mathbb{P}(\mathbf{X}=1)-\mathbb{E}\{\mathbb{E}\{\mathbf{y}\mid \mathbf{X}\}\mid \mathbf{X}=1\}\mathbb{P}(\mathbf{X}=1)=\\ & =\mathbb{E}\{\mathbf{y}\mid \mathbf{X}=1\}\mathbb{P}(\mathbf{X}=1)-\mathbb{E}\{\mathbf{y}\mid \mathbf{X}=1\}\mathbb{P}(\mathbf{X}=1)=0\end{array}$$ For a general transformation of the regressors as in Equation 14.10, the covariance is computed as: $$\begin{array}{rl}\mathbb{E}\{\mathbf{u}g(\mathbf{X})\}& =\mathbb{E}\{\mathbb{E}\{\mathbf{u}g(\mathbf{X})\mid \mathbf{X}\}\}=\\ & =\mathbb{E}\{g(\mathbf{X})\mathbb{E}\{\mathbf{u}\mid \mathbf{X}\}\}=\\ & =\mathbb{E}\{g(\mathbf{X})\mathbb{E}\{\mathbf{y}-\mathbb{E}\{\mathbf{y}\mid \mathbf{X}\}\mid \mathbf{X}\}\}=\\ & =\mathbb{E}\{g(\mathbf{X})[\mathbb{E}\{\mathbf{y}\mid \mathbf{X}\}-\mathbb{E}\{\mathbf{y}\mid \mathbf{X}\}]\}=0\end{array}$$

14.4 Uniequational linear models

Let’s consider an uni-equational linear model, i.e. with $p=1$ in (Equation 14.2), is expressed in compact matrix notation as: $$\begin{array}{}\text{(14.11)}& \underset{n\times 1}{\mathbf{y}}=\underset{n\times k}{\mathbf{X}}\underset{k\times 1}{\mathbf{b}}+\underset{n\times 1}{\mathbf{u}}\text{,}\end{array}$$ where $\mathbf{b}$ and $\mathbf{u}$ represent the true parameters and residuals in population. Let’s consider a sample of $n$-observations extracted from a population, then the matrix of the regressors $\mathbf{X}$ reads as in Equation 13.1, while the vectors of dependent variable and of the residuals reads $$\underset{n\times 1}{\mathbf{y}}=\left(\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right)\text{,}\underset{n\times 1}{\mathbf{u}}=\left(\begin{array}{c}{u}_1\\ ⋮\\ u_n\end{array}\right)\text{.}$$ Hence, the matrix of data $\mathbf{W}$ is composed by: $$\begin{array}{}\text{(14.12)}& \left(\begin{array}{cc}\mathbf{y}& \mathbf{X}\end{array}\right)=\left(\begin{array}{cccc}{y}_1& x_{1,1}& \dots & x_{1,k}\\ ⋮& ⋮& \ddots & ⋮\\ y_n& x_{n,1}& \dots & x_{n,k}\end{array}\right)\text{.}\end{array}$$

14.4.1 Estimators of b

In general, the true population parameter $\mathbf{b}$ is unknown and lives in the parameter space, i.e. $\mathbf{b}\in {\mathrm{\Theta }}_{\mathbf{b}}\subset {\mathbb{R}}^k$. In the following section we will define with $Q$ the estimator function while $\hat{\mathbf{b}}$ will denote one of its possible estimates.

Depending on the specification of the model, the function $Q$ takes as input the matrix of data and returns as output a vector of estimates, i.e.  $$Q:\mathbf{W}⟶{\mathrm{\Theta }}_{\mathbf{b}},\text{such that}Q(\mathbf{W})=\hat{\mathbf{b}}\in {\mathrm{\Theta }}_{\mathbf{b}}\text{.}$$ In this context, the fitted values $\hat{\mathbf{y}}$ can be seen as function of the estimate $\hat{\mathbf{b}}$, i.e.  $$\begin{array}{}\text{(14.13)}& \hat{\mathbf{y}}(\hat{\mathbf{b}})=\mathbf{X}\hat{\mathbf{b}}\text{.}\end{array}$$ Consequently, the fitted residuals $\hat{\mathbf{u}}$, which measure the discrepancies between the observed and the fitted values, are also a function of $\mathbf{b}$, i.e. $$\begin{array}{}\text{(14.14)}& \hat{\mathbf{u}}(\hat{\mathbf{b}})=\mathbf{y}-\hat{\mathbf{y}}(\hat{\mathbf{b}})=\mathbf{y}-\mathbf{X}\hat{\mathbf{b}}\text{.}\end{array}$$

Different optimal estimators of $\mathbf{b}$

As we will see in Chapter 15, the assumptions on the variance of the residuals determines the optimal estimator of $\mathbf{b}$. In general, the residuals could be

1. omoskedastic: the residuals are uncorrelated and their variance is equal for each observation.
2. heteroskedastic: the residuals are uncorrelated and their variance is difference for each observation.
3. autocorrelated: the residuals are correlated and their variance is equal and their variance is difference for each observation.

As shown in Figure 14.1, depending on the assumption (1,2 or 3) the optimal estimator of $\mathbf{b}$ is obtained with Ordinary Last Square for case 1, while Generalized Last Square for case 2 and 3.

Figure 14.1: Different classes of estimator for linear models.

For example, if the residuals are correlated, then their conditional variance-covariance matrix reads $$\mathbf{\Sigma }=\mathbb{V}\{\mathbf{u}{\mathbf{u}}^{\mathrm{\top }}\mid \mathbf{X}\}\text{,}$$ or more explicitly,
$$\begin{array}{}\text{(14.15)}& \underset{n\times n}{\mathbf{\Sigma }}=\left(\begin{array}{cccc}{u}_1^2& u_1{u}_2& \dots & u_1{u}_n\\ u_2{e}_1& u_2^2& \dots & u_2{u}_n\\ ⋮& ⋮& & ⋮\\ u_n{e}_1& u_n{u}_2& \dots & u_n^2\end{array}\right)=\left(\begin{array}{cccc}{\sigma }_1^2& {\sigma }_{1,2}& \dots & {\sigma }_{1,n}\\ {\sigma }_{2,1}& {\sigma }_2^2& \dots & {\sigma }_{2,n}\\ ⋮& ⋮& & ⋮\\ {\sigma }_{n,1}& {\sigma }_{n,2}& \dots & {\sigma }_n^2\end{array}\right)\text{.}\end{array}$$ Since the matrix $\mathbf{\Sigma }$ is symmetric the number of distinct elements above (or below) the diagonal reads $$(\genfrac{}{}{0ex}{}{n}{2})=\frac{n(n-1)}{2}\text{.}$$ Hence, given that the number of elements of $\mathbf{\Sigma }$ is $n\times n$, the number of unique values (free elements) are given by the $n$ variances and plus $\frac{n(n-1)}{2}$ covariances, i.e. $$\text{free elements}=n+\frac{n(n-1)}{2}=\frac{n(n+1)}{2}>n\text{.}$$