Introduzione all’Econometria

L’econometria è la scienza che si occupa di studiare relazioni economiche in modo quantitativo. L’enfasi può essere più sull’analisi e sviluppo di metodi econometrici, oppure più sulla loro applicazione: in quest’ultimo caso parliamo di econometria applicata. Questa richiede un bagaglio di conoscenze statistiche e matematiche, come l’analisi funzionale elementare (derivate, continuità, ottimizzazione ecc.), il calcolo vettoriale e matriciale, le nozioni di variabile casuale, valore atteso, funzione di densità ecc.
Queste conoscenze vengono ulteriormente sviluppate e combinate con la teoria economica per formulare modelli empirici della realtà economica. Infine, è fondamentale nell’econometria applicata la conoscenza di un appropriato software (es. Stata / Eviews / … ) o di un appropriato linguaggio (es. R / Python / …) per analizzare i dati economici.

L’obiettivo di ECONOMETRICS - PARTE 1 è quello di sintetizzare la teoria del modello di regressione multivariata nel caso di dati cross-section (\(n\) individui osservati tutti nello stesso periodo temporale): ipotesi del modello teorico, metodi di stima, proprietà degli stimatori, verifica del modello, previsioni e risultati asintotici.

L’econometria, in generale, è utile per testare la validità di teorie economiche (che specificano quantomeno la direzione della relazione tra due variabili) o per ricavare delle regolarità empiriche che forniscano spunti per sviluppare successivamente un’appropriata teoria economica. Se una teoria economica esiste già, allora dalla significatività statistica dei coefficienti di regressione è possibile anche dedurre un effetto causale. Se non esiste una teoria economica di fondamento, non è possibile inferire relazioni casuali, ma parleremo solo di correlazione tra i due fenomeni.

Iniziamo l’introduzione dal caso più semplice: modello di regressione univariata, che cerca di spiegare il comportamento di una variabile quantitativa continua \(Y\) attraverso una sola variabile esplicativa \(X\). Se \(X\) è deterministica, allora parliamo di “Fixed-X Setting”, mentre se \(X\) è stocastica parleremo di “Random-X Setting”. Qualunque sia il setting, la maggior parte dei risultati rimane valido (grazie alla legge del valore atteso reiterato e alla decomposizione della varianza).

Supponiamo che tra \(Y\) ed \(X\) esista una relazione funzionale di tipo lineare, ma che a causa della natura aleatoria di \(Y\) tale relazione non valga in modo esatto: ogni realizzazione di \(Y\) sarà spiegata dalla rispettiva realizzazione di \(X\) a meno di un termine di errore \(\epsilon\).

\[ Y = \beta_1 + \beta_{2}X + \epsilon\]

Dato che \(\epsilon = Y - \beta_1 + \beta_{2}X\), anche \(\epsilon\) è a sua volta una variabile aleatoria. Sia che si veda \(\epsilon\) come componente non spiegata di \(Y\), o come ulteriore variabile esplicativa, il trattamento del modello di regressione non cambia. La relazione di cui sopra prende il nome di modello statistico lineare. Se ad esso aggiungiamo un metodo di campionamento casuale, allora parliamo di modello di regressione lineare.

Ipotesi del modello

La relazione di cui sopra è valida all’interno della popolazione ed i valori \(\beta_1 , \beta_2\) sono scalari reali ma ignoti, che vanno stimati con opportune tecniche statistiche. Se le seguenti ipotesi (dette OLS deboli) sono valide all’interno della popolazione, allora tale ipotesi varranno anche per il modello di regressione lineare:

  1. \(Y\) è lineare e stabile nei parametri: \(\mathbb{E}(Y|X) = \beta_1 + \beta_2X\)
  2. Nessun regressore può essere ricavato come combinazione lineare degli altri
  3. \(Y\) ha varianza costante (omoschedasticità): \(\mathbb{Var}(Y_i)=\sigma^2\) per ogni \(i\)
  4. \(Y_i\) e \(Y_j\) sono non correlate: \(\mathbb{Cov}(Y_i, Y_j) = 0\) per ogni \(i \ne j\)

Queste ipotesi possono essere riformulate in termini dell’errore \(\epsilon\):

  1. L’errore \(\epsilon\) ha valore atteso condizionato nullo: \(\mathbb{E}(\epsilon|X) = \mathbb{E}(\epsilon)=0\)
  2. Nessun regressore può essere ricavato come combinazione lineare degli altri
  3. \(\epsilon\) ha varianza costante (omeschedasticità): \(\mathbb{Var}(\epsilon_i)=\sigma^2\) per ogni \(i\)
  4. \(\epsilon_i\) e \(\epsilon_j\) sono non correlati: \(\mathbb{Cov}(\epsilon_i, \epsilon_j) = 0\) per ogni \(i \ne j\)

Se le ipotesi 1-4 sono vere, allora i parametri \(\beta_1, \beta_2\) del modello possono essere stimati con OLS (ordinary least squares) e lo stimatore OLS dei parametri \(\hat\beta_{OLS}\) è BLUE (best unbiased linear estimator), dove l’aggettivo “best” indica che è lo stimatore a varianza minima all’interno della classe dei stimatori lineari corretti (non distorti).

Per stimarli con OLS, basta risolvere rispetto al vettore \(\beta = (\beta_1, \beta_2)\) il seguente problema di minimizzazione: \[\min_{\beta \in \mathbb{R}^2} \sum_{i=1}^{n}( y_i - \beta_1 - \beta_{2}x_i ) = (y - X\beta)^\mathsf{T}(y - X\beta) = L(\beta)\]

dove stavolta \(X\) indica la matrice dei dati, cioè la matrice di dimensione \(n \times p\) dove \(n\) è l’ampiezza campionaria e \(p\) è il numero di variabili esplicative, intercetta inclusa. \(y\) invece è il vettore delle realizzazioni della v.c. \(Y\) per ciascuna unità campionaria. Nel caso di regressione semplice univariata, abbiamo \(p=2\) e la matrice dei dati \(X\) di dimensione \(n \times 2\) è fatta nel seguente modo:

\[\underset{n\times 2}{\mathbf{X}} = \left[ \begin{array}{cc} 1 & x_{11}\\ 1 & x_{21}\\ \cdots & \cdots\\ 1 & x_{n1}\\ \end{array}\right]\]

Mentre più in generale, nel caso di \(p\) variabili esplicative abbiamo: \[\underset{n\times 1}{\mathbf{y}} = \left[ \begin{array}{c} y_1 \\ \cdots\\ y_i \\ \cdots\\ y_n \\ \end{array}\right] \qquad \underset{n\times p}{\mathbf{X}} = \left[ \begin{array}{cccccc} x_{1}^\mathsf{T} \\ x_{2}^\mathsf{T} \\ \cdots \\ x_{i}^\mathsf{T} \\ \cdots\\ x_{n}^\mathsf{T}\\ \end{array}\right] = \left[ \begin{array}{cccccc} x_{11} & x_{12} & \cdots & x_{1j} & \cdots & x_{1p} \\ x_{21} & x_{22} & \cdots & x_{2j} & \cdots & x_{2p} \\ \cdots & \cdots & \cdots & \cdots & \cdots \\ x_{i1} & x_{i2} & \cdots & x_{ij}& \cdots & x_{ip} \\ \cdots & \cdots & \cdots & \cdots & \cdots\\ x_{n1} & x_{n2} & \cdots & x_{nj} & \cdots & x_{np}\\ \end{array}\right]\]

con \(x_{i}^\mathsf{T}\) vettore riga delle realizzazioni delle \(p\)-variabili esplicative (intercetta inclusa) per l’\(i\)-esimo individuo / unità statistica.

Metodi di stima

Vediamo il metodo OLS in azione. L’applicazione del metodo OLS consiste nella risoluzione del problema di minimizzazione precedentemente formulato: vogliamo cercare i valori costanti per \(\beta_1\) e per \(\beta_2\) in modo che la somma dei residui al quadrato sia minima. Per risolvere il problema di minimizzazione, calcoliamo le derivate parziali della funzione obiettivo rispetto ai parametri e risolviamo rispetto al vettore di parametri le \(p\) equazioni di derivate parziali uguali a zero. Infine, si valutano le condizioni del secondo ordine: se le derivate parziali seconde sono positive quando calcolate in corrispondenza del vettore soluzione, allora tale vettore soluzione è davvero un punto di minimo.

\[\frac{\partial L(\beta)}{\partial \beta^\mathsf{T}} = \frac{\partial}{\partial\beta^\mathsf{T}}(y^\mathsf{T}y + \beta^\mathsf{T}X^\mathsf{T}X\beta - 2\beta^\mathsf{T}X^\mathsf{T}y) = 2X^\mathsf{T}X\beta - 2X^\mathsf{T}y=0\]

Che risolta per il vettore \(\beta\) dà come soluzioni:

\(\hat\beta_{OLS} = (X^\mathsf{T}X)^{-1}X^\mathsf{T}y\)

Nel caso univariato, abbiamo visto che la matrice \(X\) ha dimensione \(n \times 2\) ed il prodotto matriciale che descrive la soluzione vettoriale di \(\beta\) ha dimensione \(2 \times 1\), i.e. \(\hat\beta = (\hat\beta_0, \hat\beta_1)\).

\[\underset{n \times 2}{\mathbf{X}} = \left[ \begin{array}{cc} 1 & x_{11} \\ 1 & x_{21} \\ \cdots & \cdots \\ 1 & x_{n1} \\ \end{array}\right] = \left[ \begin{array}{cc} 1 & x \\ \end{array}\right]\]

\[ \underset{2\times n}{\mathbf{X}}^\mathsf{T} = \left[ \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{11} & x_{21} & \cdots & x_{n1} \\ \end{array}\right]\]

\[\underset{2\times 2}{\mathbf{X^\mathsf{T}\mathbf{X}}} = \left[ \begin{array}{cc} 1^\mathsf{T}1 & 1^\mathsf{T}x \\ x^\mathsf{T}1 & x^\mathsf{T}x \\ \end{array}\right] = \left[ \begin{array}{cc} n & \sum_{i=1}^{n}x_i \\ \sum_{i=1}^{n}x_i & \sum_{i=1}^{n}x_i^2 \end{array}\right]\]

\[\underset{2\times 2}{(\mathbf{X^\mathsf{T}\mathbf{X}})^{-1}} = \frac{1}{n\sum_{i=1}^{n}x_i^2 - (\sum_{i=1}^{n}x_i)^2} \times \left[ \begin{array}{cc} \sum_{i=1}^{n}x_i^2 & -\sum_{i=1}^{n}x_i \\ -\sum_{i=1}^{n}x_i & n \end{array}\right] \]

\[\underset{2\times 1}{\mathbf{X^\mathsf{T}}y} = \left[ \begin{array}{c} 1^{\mathsf{T}}y \\ x^{\mathsf{T}}y \\ \end{array}\right] = \left[ \begin{array}{c} \sum_{i=1}^{n}y_i \\ \sum_{i=1}^{n}x_{i}y_{i} \end{array}\right]\]

Infine, moltiplicando nell’ordine corretto le matrici, otteniamo gli stimatori OLS per \(\beta_0\) e \(\beta_1\). \

\[\hat\beta_{OLS} = \frac{1}{n\sum_{i=1}^{n}x_i^2 - (\sum_{i=1}^{n}x_i)^2} \times \left[ \begin{array}{cc} \sum_{i=1}^{n}x_i^2 & -\sum_{i=1}^{n}x_i \\ -\sum_{i=1}^{n}x_i & n \end{array}\right] \times \left[ \begin{array}{c} \sum_{i=1}^{n}y_i \\ \sum_{i=1}^{n}x_{i}y_{i} \end{array}\right] = \left[ \begin{array}{c} \hat\beta_1 \\ \hat\beta_2 \\ \end{array}\right]\]

dove il parametro stimato \(\beta_2\) misura la pendenza della retta di regressione, i.e. l’effetto marginale su \(Y\) di un incremento unitario di \(X\), mentre il parametro \(\beta_1\) rappresenta il valor medio di \(Y\) quando \(X=0\).

\[\hat\beta_2 = \frac{n\sum_{i=1}^{n}x_{i}y_{i} - \sum_{i=1}^{n}x_{i}\sum_{i=1}^{n}y_i}{n\sum_{i=1}^{n}x_i^2 - (\sum_{i=1}^{n}x_i)^2} \qquad \hat\beta_1 = n^{-1}(\sum_{i=1}^{n}y_i - \hat\beta_2\sum_{i=1}^{n}x_i)=\bar{y} - \hat\beta_2\bar{x}\] Un altro metodo di stima consiste nel metodo della massima verosimiglianza, ma richiede l’ipotesi OLS forte secondo cui è specificata la distribuzione probabilistica della variabile dipendente \(Y\) come una Normale di media \(X\beta\) e varianza \(\sigma^2\). In tal caso, nota la distribuzione di Y condizionata ad X, possiamo calcolare la likelihood function come prodotto delle n-marginali (una per ciascuna osservazione) e massimizzare tale funzione rispetto al vettore di parametri \(\beta\). Si dimostra che \(\beta_{ols}\) e \(\beta_{ML}\) sono equivalenti, mentre la varianza \(\sigma^2\) stimata con ML è distorta.

Infine, un ulteriore metodo è quello dei momenti (stima MM), dove si uguaglia il momento teorico della distribuzione di Y|X al rispettivo momento campionario. In questo caso, si uguaglia la media teorica \(x_i^\mathsf{T}\beta\) (funzione del vettore \(\beta\)) con la media campionaria condizionata \(\mathbb{E}(y_i|X=x_i)\) e si risolve l’equazione (o sistema di equazioni) rispetto al parametro \(\beta\). Se \(p=2\) servono almeno due realizzazioni di \(X\) con valori diversi, diciamo \(x_i\) e \(x_j\), in modo da calcolare: (1 x_i)*(beta1 beta2) = beta1 + beta2 x_i == mean(y | X = x_i) ed in modo analogo (1 x_j)(beta1 beta2) = beta1 + beta2 x_j == mean(y|X=x_j) e risolvere le due equazioni rispetto alle due incognite beta1, beta2. Purtroppo è un metodo abbastanza arbitrario: non è garantita la consistenza dello stimatore e, in alcuni casi, potrebbe non esistere finito il momento teorico di ordine r>1.

# Dimostra che le formule risolutive per beta1 e beta2 sono equivalenti al calcolo del vettore beta_ols
set.seed(321)
n = 100 
x = rnorm(n, mean=10, sd=1)
epsilon = rnorm(n)
y = 5 + 2.5*x + epsilon # true beta = c(5, 2.5); true sigma = 1

beta_2 = (n*sum(x*y) - sum(x)*sum(y))/(n*sum(x^2) - (sum(x))^2)
beta_1 = mean(y) - beta_2*mean(x)

X = as.matrix(cbind(rep(1,n), x))
beta_ols = solve(t(X)%*%X)%*%t(X)%*%y
round(beta_ols,2) == round(c(beta_1, beta_2),2)
##   [,1]
##   TRUE
## x TRUE
# Dimostra che beta_ols e beta_ML (max likelihood) coincidono quando y ~ N(X*beta, sigma), 
# In modo numerico: inizializza con start = beta_ols per avere uguaglianza esatta

start = beta_ols

nlogL <- function(param, x, y) {
0.5*sum((y - param[1] - param[2]*x)^2)
}

Vectorize(nlogL, vectorize.args='param')
## function (param, x, y) 
## {
##     args <- lapply(as.list(match.call())[-1L], eval, parent.frame())
##     names <- if (is.null(names(args))) 
##         character(length(args))
##     else names(args)
##     dovec <- names %in% vectorize.args
##     do.call("mapply", c(FUN = FUN, args[dovec], MoreArgs = list(args[!dovec]), 
##         SIMPLIFY = SIMPLIFY, USE.NAMES = USE.NAMES))
## }
## <environment: 0x000000001493e250>
ML <- optim(start, fn=nlogL, y=y, x=x, hessian=T)
stima <- ML$par

data.frame(row.names=c("beta_1", "beta_2"), ML=stima, OLS=beta_ols, True=c(5, 2.5))
##              ML      OLS True
## beta_1 4.100646 4.100646  5.0
## beta_2 2.589699 2.589699  2.5
# In modo analitico: risolvo min(nlogL) wrt beta
beta_ml2 = sum((x - mean(x))*(y - mean(y)))/sum((x - mean(x))^2)
beta_ml1 = mean(y) - mean(x)*beta_ml2
data.frame(row.names = c("beta_1", "beta_2"), ML=c(beta_ml1, beta_ml2), OLS=beta_ols, True=c(5, 2.5))
##              ML      OLS True
## beta_1 4.100646 4.100646  5.0
## beta_2 2.589699 2.589699  2.5
# Dimostra che il coefficiente di determinazione R2 coincide con il coefficiente di correlazione al quadrato nel caso univariato
r = cor(x,y)
mod = lm(y ~ x)
R2 = summary(mod)$r.squared
c(R2, r^2)
## [1] 0.8806638 0.8806638
# Dimostra che la correlazione tra y e yhat al quadrato è sempre pari al coefficiente di determinazione R2

yhat = predict(mod, as.data.frame(X))
round(cor(y, yhat)^2,10) == round(R2,10)
## [1] TRUE