Laboratori in R
Ambiente e convenzioni dei laboratori
I laboratori utilizzano R e le funzioni dei pacchetti base e stats, senza richiedere archivi esterni. I dati sono piccoli esempi deterministici oppure simulazioni con processo generatore esplicito. I risultati teorici del processo non devono essere confusi con le stime casuali prodotte da una singola simulazione. La documentazione ufficiale di lm, glm, arima, optim e delle distribuzioni è raccolta nel riferimento R29. La registrazione di sessionInfo() permette di conservare le versioni dell’ambiente utilizzato.
I blocchi vanno eseguiti nell’ordine indicato quando richiamano una funzione definita in precedenza. Le funzioni di covarianza seguenti sono destinate a regressioni OLS non ponderate, a rango pieno e prive di osservazioni mancanti nel campione fornito. Non sono sostituti universali di implementazioni specialistiche per pesi, molti effetti assorbiti, pochi cluster o modelli non lineari. Le condizioni vengono rese esplicite per evitare che una formula corretta in un caso venga applicata a un altro modello soltanto perché accetta gli stessi oggetti.
Laboratorio 1 OLS, residui e previsioni
Si considerano cinque osservazioni, sufficienti per ricostruire a mano le formule e confrontarle con la regressione. L’intercetta stimata è 2,2, la pendenza è 0,6, RSS è 2,4 e R2 è 0,6. Le verifiche numeriche inserite nel codice interrompono l’esecuzione se questi risultati non sono riprodotti entro la tolleranza. Gli intervalli richiesti a predict sono distinti: confidence riguarda la media condizionata, prediction una nuova osservazione sotto le ipotesi classiche.
# R: regressione lineare elementare
x <- 1:5
y <- c(2, 4, 5, 4, 5)
dat <- data.frame(x = x, y = y)
fit <- lm(y ~ x, data = dat)
b <- coef(fit)
rss <- sum(residuals(fit)^2)
tss <- sum((y - mean(y))^2)
stopifnot(isTRUE(all.equal(unname(b), c(2.2, 0.6))))
stopifnot(abs(rss - 2.4) < 1e-10)
print(summary(fit))
print(confint(fit))
print(c(RSS = rss, R2 = 1 - rss / tss))
new <- data.frame(x = 3)
print(predict(fit, newdata = new, interval = "confidence"))
print(predict(fit, newdata = new, interval = "prediction"))
La previsione puntuale in x = 3 è quattro. L’errore standard della media è 0,4, mentre quello di una nuova osservazione è la radice di 0,96, circa 0,9798. La differenza deriva dal rumore individuale aggiuntivo. Il campione ha soltanto tre gradi di libertà residui, quindi gli intervalli classici usano una t con tre gradi di libertà e sono molto più ampi di quelli ottenuti sostituendo automaticamente 1,96.
Laboratorio 2 Covarianza HC3
La funzione costruisce la matrice sandwich HC3 usando la matrice del modello, i residui e il leverage. Controlla rango, gradi di libertà e denominatori prossimi a zero. Il successivo esempio genera una media lineare corretta ma una dispersione crescente con x. OLS e la versione con HC3 hanno gli stessi coefficienti: ciò che cambia è la misura dell’incertezza. La simulazione non garantisce che un errore standard sia sempre maggiore dell’altro in ogni realizzazione.
# R: covarianza HC3 per OLS non ponderati
vcov_hc3 <- function(fit) {
if (!inherits(fit, "lm") || inherits(fit, "glm"))
stop("Serve un modello OLS lm, non glm.")
if (!is.null(fit$weights)) stop("Pesi non supportati.")
X <- model.matrix(fit)
e <- fit$residuals
n <- nrow(X); k <- ncol(X)
if (n <= k || qr(X)$rank < k) stop("Rango o gradi di liberta insufficienti.")
h <- hatvalues(fit)
if (any(1 - h <= sqrt(.Machine$double.eps)))
stop("Leverage troppo vicino a uno.")
B <- solve(crossprod(X))
omega <- (e / (1 - h))^2
B %*% crossprod(X, X * as.numeric(omega)) %*% B
}
set.seed(2026)
n <- 400
x <- runif(n, 0, 4)
y <- 1 + 2 * x + (0.3 + x) * rnorm(n)
fit_h <- lm(y ~ x)
V_hc3 <- vcov_hc3(fit_h)
print(cbind(stima = coef(fit_h),
SE_classico = sqrt(diag(vcov(fit_h))),
SE_HC3 = sqrt(diag(V_hc3))))
Un intervallo asintotico normale può essere costruito con coef(fit_h) più o meno qnorm(0.975) volte gli errori standard robusti. Non sarebbe un intervallo t esatto nel modello eteroschedastico. La funzione non corregge una covariata endogena: se x fosse correlata con il termine non osservato, la stima della pendenza potrebbe convergere a un valore diverso da due anche con una matrice HC3 calcolata correttamente.
Laboratorio 3 IV e due stadi con covarianza coerente
Si genera x attraverso uno strumento z e una componente v. L’errore dell’esito contiene v, rendendo x endogena, mentre z rimane indipendente dall’errore per costruzione. Il parametro strutturale è due. La pendenza OLS converge a 2 + 0,7/1,64, circa 2,4268, sotto il processo specificato. IV usa soltanto la variazione di x collegata a z. La covarianza robusta viene costruita con i residui strutturali y − Xb, non con i residui della regressione sui valori previsti.
# R: 2SLS con uno strumento escluso
set.seed(2027)
n <- 500
z <- rnorm(n); v <- rnorm(n)
u <- 0.7 * v + rnorm(n)
x <- 0.8 * z + v
y <- 1 + 2 * x + u
X <- cbind(intercetta = 1, x = x)
Z <- cbind(intercetta = 1, z = z)
if (qr(Z)$rank < ncol(Z)) stop("Strumenti ridondanti.")
Xhat <- Z %*% qr.coef(qr(Z), X)
if (qr(Xhat)$rank < ncol(Xhat)) stop("Modello non identificato.")
A <- solve(crossprod(Xhat))
b_iv <- A %*% crossprod(Xhat, y)
e_iv <- as.numeric(y - X %*% b_iv)
V_iv <- A %*% crossprod(Xhat, Xhat * e_iv^2) %*% A
print(cbind(stima_IV = as.numeric(b_iv), SE_HC0 = sqrt(diag(V_iv))))
print(coef(lm(y ~ x)))
print(summary(lm(x ~ z)))
La forza del primo stadio deve essere valutata insieme alla numerosità e all’inferenza desiderata. Ridurre il coefficiente 0,8 dello strumento può rendere il rapporto IV instabile. La simulazione rende vera l’esclusione perché il processo è noto; in dati reali questa proprietà richiede una giustificazione economica. Il fatto che il codice produca una stima vicina a due in alcuni campioni non dimostra la validità di uno strumento osservazionale.
Laboratorio 4 Effetti fissi e clustering
La funzione cluster seguente consente dipendenza entro gruppi e applica una correzione CR1 elementare. Richiede che il vettore dei gruppi sia già allineato alle righe effettivamente utilizzate dal modello. In questo esempio non vi sono dati mancanti. L’inferenza con pochi gruppi richiederebbe procedure ulteriori; la funzione non le implementa. Il panel simulato contiene una componente individuale correlata con x e un coefficiente strutturale pari a 1,5.
# R: covarianza cluster per OLS non ponderati
vcov_cluster <- function(fit, cluster) {
if (!inherits(fit, "lm") || inherits(fit, "glm"))
stop("Serve un modello OLS lm.")
if (!is.null(fit$weights)) stop("Pesi non supportati.")
X <- model.matrix(fit); e <- fit$residuals
n <- nrow(X); k <- ncol(X)
if (length(cluster) != n || anyNA(cluster))
stop("Identificativi cluster non allineati.")
G <- length(unique(cluster))
if (G < 2 || n <= k || qr(X)$rank < k)
stop("Informazione insufficiente per la covarianza.")
S <- rowsum(X * as.numeric(e), group = cluster, reorder = FALSE)
B <- solve(crossprod(X))
(G / (G - 1)) * ((n - 1) / (n - k)) *
B %*% crossprod(S) %*% B
}
set.seed(2028)
N <- 80; TT <- 5
id <- rep(seq_len(N), each = TT)
time <- rep(seq_len(TT), times = N)
a <- rnorm(N)
x <- 0.6 * a[id] + rnorm(N * TT)
y <- 3 + 1.5 * x + a[id] + 0.2 * time + rnorm(N * TT)
panel <- data.frame(id, time, x, y)
stopifnot(!anyDuplicated(panel[c("id", "time")]))
fe <- lm(y ~ x + factor(id) + factor(time), data = panel)
V_fe <- vcov_cluster(fe, panel$id)
print(c(stima = coef(fe)["x"], SE_cluster = sqrt(V_fe["x", "x"])))
Le dummies individuali eliminano la componente a nel confronto della pendenza. L’inclusione di effetti temporali rappresenta differenze comuni fra date. Il clustering non modifica il coefficiente, ma consente che gli errori della stessa unità non siano indipendenti. Per grandi panel, l’assorbimento degli effetti è computazionalmente più efficiente della costruzione esplicita delle dummies; l’esempio utilizza la forma esplicita per rendere visibile il modello.
Laboratorio 5 Differenze nelle differenze
Il trattamento è assegnato a metà delle unità e inizia dal terzo periodo. Il processo simulato ha tendenze parallele per costruzione e un effetto omogeneo pari a due. La regressione contiene effetti individuali e temporali e usa la funzione cluster del laboratorio precedente. Il risultato campionario varia attorno all’effetto del processo; non si deve sostituire il numero due al coefficiente effettivamente stimato nell’output.
# R: DiD con panel e adozione comune
set.seed(2029)
N <- 100; TT <- 4
id <- rep(seq_len(N), each = TT)
time <- rep(seq_len(TT), N)
group <- as.integer(id <= N / 2)
post <- as.integer(time >= 3)
D <- group * post
a <- rnorm(N)
y <- a[id] + 0.5 * time + 2 * D + rnorm(N * TT)
did_data <- data.frame(id, time, group, post, D, y)
fit_did <- lm(y ~ D + factor(id) + factor(time), data = did_data)
V_did <- vcov_cluster(fit_did, did_data$id)
print(c(effetto = coef(fit_did)["D"],
SE_cluster = sqrt(V_did["D", "D"])))
print(aggregate(y ~ group + post, data = did_data, FUN = mean))
Il laboratorio non copre l’adozione scaglionata con effetti eterogenei: in quel contesto occorrono estimandi gruppo-tempo e confronti coerenti, non una sostituzione automatica della variabile D. La simulazione può essere modificata aggiungendo un trend specifico del gruppo trattato. In quel caso la regressione può attribuire al trattamento parte del cambiamento differenziale non causato da esso, mostrando concretamente il ruolo delle tendenze parallele.
Laboratorio 6 Logit, probabilità e poisson con offset
La prima simulazione genera un esito binario da un logit. Il coefficiente dell’indice non è una differenza di probabilità; l’effetto marginale medio di x, in questa specificazione senza trasformazioni o interazioni di x, è la media di p(1 − p) moltiplicata per il coefficiente. La seconda simulazione genera conteggi con esposizioni diverse e un tasso esponenziale. L’offset mantiene distinta la durata di esposizione dalla relazione delle covariate con il tasso.
# R: logit ed effetto marginale medio
set.seed(2030)
n <- 600
x <- rnorm(n)
p <- plogis(-0.4 + 0.7 * x)
y_bin <- rbinom(n, size = 1, prob = p)
fit_logit <- glm(y_bin ~ x, family = binomial())
p_hat <- fitted(fit_logit)
ame_x <- mean(p_hat * (1 - p_hat)) * coef(fit_logit)["x"]
print(c(coefficiente = coef(fit_logit)["x"],
odds_ratio = exp(coef(fit_logit)["x"]), AME = ame_x))
# R: Poisson con esposizione e covarianza sandwich HC0
exposure <- runif(n, 0.5, 3)
mu <- exposure * exp(0.2 + 0.3 * x)
y_count <- rpois(n, lambda = mu)
fit_p <- glm(y_count ~ x + offset(log(exposure)), family = poisson())
Xp <- model.matrix(fit_p)
mup <- fitted(fit_p)
ep <- y_count - mup
Ap <- solve(crossprod(Xp, Xp * as.numeric(mup)))
Vp <- Ap %*% crossprod(Xp, Xp * ep^2) %*% Ap
print(cbind(stima = coef(fit_p), SE_robusto = sqrt(diag(Vp))))
La matrice robusta Poisson non è la covarianza HC3 OLS: usa la curvatura del criterio Poisson e i suoi contributi di score. La simulazione soddisfa equidispersione, ma la covarianza sandwich consente di studiare anche distribuzioni diverse con media corretta e condizioni adeguate. Per probabilità di zero e intervalli distributivi serve un modello della distribuzione, non soltanto della media. L’effetto marginale logit calcolato è una stima puntuale; il suo errore standard richiederebbe metodo delta o un ricampionamento coerente.
Laboratorio 7 Previsioni a origine mobile
Si genera un AR(1) stazionario con media dieci e persistenza 0,7. Le ultime quaranta osservazioni vengono previste una alla volta utilizzando soltanto il passato disponibile a ciascuna origine. Il benchmark usa l’ultimo valore osservato. La divisione non è casuale, perché il problema riguarda il futuro. L’intercetta riportata da arima in questa parametrizzazione senza differenziazione è collegata alla media della serie: non va confusa con il termine costante c dell’equazione autoregressiva scritta in livelli.
# R: valutazione temporale senza informazione futura
set.seed(2031)
y_ts <- 10 + as.numeric(arima.sim(model = list(ar = 0.7), n = 240))
test_index <- 201:240
pred_ar <- pred_naive <- numeric(length(test_index))
for (j in seq_along(test_index)) {
target <- test_index[j]
train <- y_ts[seq_len(target - 1)]
mod <- arima(train, order = c(1, 0, 0),
include.mean = TRUE, method = "ML")
if (mod$code != 0) stop("Ottimizzazione AR non convergente.")
pred_ar[j] <- as.numeric(predict(mod, n.ahead = 1)$pred)
pred_naive[j] <- tail(train, 1)
}
actual <- y_ts[test_index]
metrics <- function(actual, predicted) {
e <- actual - predicted
c(MAE = mean(abs(e)), RMSE = sqrt(mean(e^2)))
}
print(rbind(AR1 = metrics(actual, pred_ar),
Naive = metrics(actual, pred_naive)))
La superiorità teorica della previsione condizionata sotto il processo noto non garantisce che il modello stimato vinca su ogni campione di quaranta errori. Stima dei parametri e variabilità della valutazione possono modificare l’ordinamento osservato. Per confrontare formalmente accuratezze servirebbe considerare i differenziali di perdita e la loro dipendenza. Il codice non usa il test per scegliere i ritardi; aggiungere una selezione richiederebbe di ripeterla all’interno di ogni campione di addestramento pertinente.
Laboratorio 8 Bootstrap delle coppie e randomizzazione esatta
Il primo esercizio ricampiona intere coppie di regressore ed esito. Il secondo enumera tutte le assegnazioni possibili di due trattamenti fra quattro unità. Sono due fonti di casualità differenti: il bootstrap approssima il campionamento dalla popolazione, mentre il test esatto usa il meccanismo di assegnazione sotto una nulla forte. Le due distribuzioni non devono essere interpretate come se rispondessero alla stessa domanda.
# R: bootstrap i.i.d. delle coppie
set.seed(2032)
n <- 150
boot_data <- data.frame(x = rnorm(n))
boot_data$y <- 1 + 0.8 * boot_data$x + rnorm(n)
B <- 999
boot_b <- replicate(B, {
idx <- sample.int(n, n, replace = TRUE)
db <- boot_data[idx, , drop = FALSE]
if (length(unique(db$x)) < 2) return(NA_real_)
unname(coef(lm(y ~ x, data = db))["x"])
})
if (anyNA(boot_b)) warning("Repliche degeneri escluse.")
print(c(stima = coef(lm(y ~ x, data = boot_data))["x"],
SE_bootstrap = sd(boot_b, na.rm = TRUE)))
print(quantile(boot_b, c(0.025, 0.975), na.rm = TRUE))
# R: distribuzione esatta sotto la nulla forte
outcomes <- c(1, 2, 5, 6)
assignments <- combn(seq_along(outcomes), 2)
stat <- apply(assignments, 2, function(treated) {
mean(outcomes[treated]) - mean(outcomes[-treated])
})
observed <- mean(outcomes[c(3, 4)]) - mean(outcomes[c(1, 2)])
p_exact <- mean(abs(stat) >= abs(observed))
stopifnot(abs(p_exact - 1 / 3) < 1e-12)
print(sort(stat))
print(p_exact)
L’intervallo percentile bootstrap non è un intervallo esatto in qualunque campione. Il test di randomizzazione ha invece un riferimento esatto per le sei assegnazioni equiprobabili e la nulla forte specificata. Il p-value un terzo non viene ridotto aumentando il numero di simulazioni: tutte le assegnazioni sono già state enumerate. Aumentare artificialmente le repliche non crea nuove assegnazioni né nuove unità sperimentali.
Laboratorio 9 GMM lineare sovraidentificato
Il modello contiene due parametri strutturali, intercetta e pendenza, e tre strumenti complessivi, costante e due strumenti esclusi. Esiste quindi una restrizione di sovraidentificazione. Il primo stadio GMM usa ponderazione identità; il secondo usa la covarianza stimata dei momenti. Le normalizzazioni sono mantenute in termini di medie, perciò la covarianza finale dei parametri contiene una divisione per n e la statistica J una moltiplicazione per n.
# R: GMM lineare a due stadi
set.seed(2033)
n <- 500
z1 <- rnorm(n); z2 <- rnorm(n); v <- rnorm(n)
x <- 0.8 * z1 + 0.4 * z2 + v
y <- 1 + 2 * x + 0.6 * v + rnorm(n)
X <- cbind(1, x)
Z <- cbind(1, z1, z2)
if (qr(Z)$rank < ncol(Z)) stop("Strumenti ridondanti.")
Q <- crossprod(Z, X) / n
q <- crossprod(Z, y) / n
if (qr(Q)$rank < ncol(X)) stop("Condizione di rango non soddisfatta.")
b1 <- solve(crossprod(Q), crossprod(Q, q))
e1 <- as.numeric(y - X %*% b1)
S1 <- crossprod(Z * e1) / n
W2 <- solve(S1)
A <- solve(t(Q) %*% W2 %*% Q)
b2 <- A %*% t(Q) %*% W2 %*% q
e2 <- as.numeric(y - X %*% b2)
S2 <- crossprod(Z * e2) / n
V2 <- A %*% t(Q) %*% W2 %*% S2 %*% W2 %*% Q %*% A / n
gbar <- colMeans(Z * e2)
J <- as.numeric(n * t(gbar) %*% W2 %*% gbar)
df_J <- ncol(Z) - ncol(X)
print(cbind(stima = as.numeric(b2), SE = sqrt(diag(V2))))
print(c(J = J, gradi_liberta = df_J,
p_value = pchisq(J, df = df_J, lower.tail = FALSE)))
Il riferimento chi-quadrato di J è asintotico e richiede validità dei momenti e regolarità. Un p-value elevato non dimostra che gli strumenti siano validi. Nella simulazione la validità deriva dalla costruzione indipendente degli strumenti, non dal test. Si può studiare una violazione aggiungendo direttamente z1 all’equazione dell’esito: il primo strumento diventerebbe invalido per il modello che lo esclude, ma il test potrebbe comunque non rifiutare in alcune realizzazioni.
Laboratorio 10 AIPW con cross-fitting
Si simula un trattamento selezionato attraverso x, con effetto costante pari a due e relazione dell’esito non lineare in x. Il propensity score e le funzioni dell’esito vengono stimati fuori dal fold nel quale si calcolano i contributi AIPW. Il codice mantiene le probabilità stimate, senza troncarle silenziosamente; valori estremi generano un avviso. La formula dell’errore standard è quella basata sui contributi individuali nel quadro i.i.d. regolare, non una correzione per clustering o per confondimento non osservato.
# R: AIPW con cinque fold
set.seed(2034)
n <- 600; K <- 5
x <- rnorm(n)
e_true <- plogis(0.2 + 0.4 * x)
D <- rbinom(n, 1, e_true)
y <- 1 + 2 * D + x + 0.5 * x^2 + rnorm(n)
d <- data.frame(y, D, x)
fold <- sample(rep(seq_len(K), length.out = n))
e_hat <- m0_hat <- m1_hat <- rep(NA_real_, n)
for (k in seq_len(K)) {
te <- which(fold == k)
tr <- d[fold != k, , drop = FALSE]
if (length(unique(tr$D)) != 2) stop("Manca un gruppo nel training.")
ps <- glm(D ~ x, data = tr, family = binomial())
m0 <- lm(y ~ x + I(x^2), data = tr[tr$D == 0, , drop = FALSE])
m1 <- lm(y ~ x + I(x^2), data = tr[tr$D == 1, , drop = FALSE])
e_hat[te] <- predict(ps, newdata = d[te, ], type = "response")
m0_hat[te] <- predict(m0, newdata = d[te, ])
m1_hat[te] <- predict(m1, newdata = d[te, ])
}
if (anyNA(e_hat) || any(e_hat <= 0 | e_hat >= 1))
stop("Propensity score non utilizzabile.")
if (min(e_hat) < 0.01 || max(e_hat) > 0.99)
warning("Probabilita estreme: verificare sovrapposizione e pesi.")
score <- m1_hat - m0_hat + D * (y - m1_hat) / e_hat -
(1 - D) * (y - m0_hat) / (1 - e_hat)
ate <- mean(score)
se_ate <- sd(score) / sqrt(n)
print(c(ATE = ate, SE = se_ate,
inferiore = ate - qnorm(0.975) * se_ate,
superiore = ate + qnorm(0.975) * se_ate))
Le specificazioni delle funzioni ausiliarie coincidono qui con il processo generatore, ma una singola stima non deve essere esattamente due. Modificare una delle due componenti permette di esplorare la doppia robustezza in simulazioni ripetute. Modificarle entrambe o omettere un confondente che influenza trattamento ed esito può produrre distorsione. Cross-fitting e flessibilità statistica non sostituiscono l’ignorabilità e la positività richieste per identificare l’ATE.
Laboratorio 11 Ricorsione e quasi-verosimiglianza GARCH
Il laboratorio considera una serie di innovazioni a media zero, evitando di aggiungere una stima della media per concentrare l’attenzione sulla varianza. La parametrizzazione impone ω positivo, α e β positivi e α + β inferiore a uno. Questa restrizione rappresenta il modello con secondo momento finito adottato nell’esercizio; non comprende tutte le possibili strutture GARCH. La varianza iniziale viene posta al livello incondizionato del modello candidato.
# R: simulazione e stima GARCH(1,1) a media zero
set.seed(2035)
n <- 1000
z <- rnorm(n)
h <- eps <- numeric(n)
h[1] <- 1
for (t in seq_len(n)) {
if (t > 1) h[t] <- 0.1 + 0.1 * eps[t - 1]^2 + 0.8 * h[t - 1]
eps[t] <- sqrt(h[t]) * z[t]
}
unpack <- function(par) {
w <- exp(par[1])
a <- c(0, par[2], par[3])
pr <- exp(a - max(a)); pr <- pr / sum(pr)
c(omega = w, alpha = pr[2], beta = pr[3])
}
nll <- function(par, y) {
p <- unpack(par)
if (any(!is.finite(p)) || sum(p[2:3]) >= 1) return(1e100)
ht <- numeric(length(y))
ht[1] <- p[1] / (1 - p[2] - p[3])
for (t in 2:length(y))
ht[t] <- p[1] + p[2] * y[t - 1]^2 + p[3] * ht[t - 1]
if (any(!is.finite(ht)) || any(ht <= 0)) return(1e100)
ans <- 0.5 * sum(log(2 * pi) + log(ht) + y^2 / ht)
if (!is.finite(ans)) return(1e100)
ans
}
res <- optim(c(log(0.1), 0, log(8)), nll, y = eps,
method = "BFGS", control = list(maxit = 2000))
if (res$convergence != 0) stop("Stima GARCH non convergente.")
print(unpack(res$par))
Il codice produce stime puntuali, non una matrice di covarianza robusta dei parametri. Per costruire intervalli occorrerebbe propagare la trasformazione e utilizzare un’inferenza coerente con la distribuzione delle innovazioni e con la procedura di stima. Anche una convergenza numerica riuscita richiede controlli dei residui standardizzati, sensibilità alle inizializzazioni e verifica della plausibilità del modello. La restrizione α + β < 1 non deve essere imposta soltanto per nascondere una persistenza empirica prossima al confine.
Laboratorio 12 Posterior coniugate e previsione
Le distribuzioni coniugate permettono di separare nettamente il calcolo esatto del modello dall’approssimazione numerica di una simulazione. Nel caso beta-binomiale, i quantili della posterior sono calcolati direttamente. Nel caso normale, si confrontano intervallo della media e intervallo predittivo di una nuova osservazione. I numeri sono condizionati alle prior e alle distribuzioni adottate, non probabilità prive di ipotesi.
# R: beta-binomiale
prior_a <- 2; prior_b <- 2
successi <- 7; prove <- 10
post_a <- prior_a + successi
post_b <- prior_b + prove - successi
print(c(media = post_a / (post_a + post_b),
moda = (post_a - 1) / (post_a + post_b - 2)))
print(qbeta(c(0.025, 0.975), post_a, post_b))
# R: media normale con varianza nota
m0 <- 0; s0_sq <- 4
n <- 9; sigma_sq <- 9; ybar <- 3
sn_sq <- 1 / (1 / s0_sq + n / sigma_sq)
mn <- sn_sq * (m0 / s0_sq + n * ybar / sigma_sq)
zcrit <- qnorm(0.975)
print(c(media_posterior = mn, varianza_posterior = sn_sq))
print(mn + c(-1, 1) * zcrit * sqrt(sn_sq))
print(mn + c(-1, 1) * zcrit * sqrt(sigma_sq + sn_sq))
# Conservare la configurazione dell'ambiente insieme ai risultati.
sessionInfo()
L’ultima istruzione registra l’ambiente software, ma non rende automaticamente riproducibili fonti esterne che possono cambiare. In un progetto reale vanno conservate anche versioni dei dati, definizioni del campione e istruzioni di costruzione delle variabili. Questi laboratori sono autosufficienti proprio perché utilizzano dati deterministici o processi simulati espliciti. Il passaggio a dati osservazionali richiede una valutazione sostanziale delle ipotesi, non soltanto la sostituzione del nome del dataset.
Esercizi con soluzione
Esercizio 1 Media della popolazione e composizione
Una popolazione comprende trenta unità del gruppo A e settanta del gruppo B. Le medie dell’esito sono rispettivamente dieci e quattordici. Si chiede la media complessiva e si vuole valutare la proposta di fare la media semplice delle due medie.
Soluzione. La media complessiva è (30 × 10 + 70 × 14)/100 = 12,8. La media semplice delle medie sarebbe dodici e attribuirebbe implicitamente lo stesso peso ai due gruppi, come se rappresentassero metà della popolazione ciascuno. I due numeri rispondono quindi a popolazioni di riferimento diverse. Se il campione contenesse deliberatamente uguali numerosità per gruppo, ma l’obiettivo restasse la popolazione trenta-settanta, occorrerebbe ripristinare i pesi pertinenti. Un calcolo formalmente semplice può essere sbagliato quando il denominatore non rappresenta l’estimando.
Esercizio 2 Legge della varianza totale
Nella popolazione dell’esercizio precedente, le varianze interne ai gruppi A e B sono quattro e nove. Determinare la varianza complessiva e distinguere dispersione interna e differenze fra gruppi.
Soluzione. La legge della varianza totale dà Var(Y) = E[Var(Y|G)] + Var[E(Y|G)]. La componente interna è 0,3 × 4 + 0,7 × 9 = 7,5. La componente fra gruppi è 0,3(10 − 12,8)2 + 0,7(14 − 12,8)2 = 3,36. La varianza complessiva è 10,86. Fare soltanto la media ponderata delle varianze produrrebbe 7,5 e ometterebbe la dispersione dovuta alle diverse medie. Questo principio spiega perché la variabilità aggregata possa essere elevata anche quando ogni sottogruppo è relativamente omogeneo. Le formule riguardano qui momenti di popolazione; le correzioni dei divisori delle varianze campionarie richiedono attenzione separata.
Esercizio 3 Ricostruire OLS da statistiche sufficienti
Un campione di sei osservazioni ha x̄ = 4, ȳ = 10, Sxx = 20, Sxy = 30 e TSS = 60. Calcolare intercetta, pendenza, RSS, R2 e stima classica della varianza dell’errore.
Soluzione. La pendenza è 30/20 = 1,5 e l’intercetta è 10 − 1,5 × 4 = 4. La somma spiegata è β̂1Sxy = 45, quindi RSS = 60 − 45 = 15 e R2 = 45/60 = 0,75. Con due parametri stimati, i gradi di libertà sono quattro e s2 = 15/4 = 3,75. Sotto omoschedasticità, Var stimata della pendenza è 3,75/20 = 0,1875 e il suo errore standard circa 0,4330. Il risultato OLS è una descrizione del campione; per attribuire alla pendenza proprietà inferenziali o causali occorrono ipotesi aggiuntive. R2 pari a 0,75 non le dimostra.
Esercizio 4 Cambiare unità senza cambiare il modello
Il salario, espresso in migliaia di euro, è previsto da ŷ = 8 + 0,6x, dove x misura anni di una determinata formazione. Riscrivere il modello con salario in euro e durata in mesi.
Soluzione. Indicando con Y€ = 1000Y e con M = 12X, si ottiene Ŷ€ = 8000 + 600X = 8000 + 50M. L’intercetta diventa ottomila euro e la pendenza cinquanta euro per mese. Le previsioni economiche sono le stesse: dodici mesi aggiuntivi corrispondono a seicento euro nella componente prevista. Errori standard e intervalli dei coefficienti devono essere riscalati nello stesso modo, mentre statistiche t per ipotesi equivalenti e R2 non cambiano. Il cambiamento di unità non modifica l’identificazione: una relazione associativa non diventa causale quando il coefficiente è espresso in una scala più intuitiva.
Esercizio 5 Ipotesi nulla e intervallo
Una pendenza stimata è 0,12, con errore standard classico 0,04 e ottanta gradi di libertà. Calcolare la statistica per H0: β = 0 e quella per H0: β = 0,10. Usare 1,990 come valore critico bilaterale approssimato al 5%.
Soluzione. Per la prima nulla, t = 0,12/0,04 = 3, quindi si rifiuta al livello indicato nel quadro classico appropriato. Per la seconda, t = (0,12 − 0,10)/0,04 = 0,5 e non si rifiuta. L’intervallo è 0,12 ± 1,990 × 0,04, circa [0,0404; 0,1996]. Zero è escluso, 0,10 è incluso. La significatività rispetto a zero non dimostra che il coefficiente sia maggiore di 0,10. La domanda economica deve essere tradotta nella nulla pertinente, non sostituita automaticamente dal test predefinito del software.
Esercizio 6 Test f parziale
Un modello non ristretto con cinque parametri, compresa l’intercetta, usa cento osservazioni e ha RSS = 180. Il modello ristretto impone due vincoli lineari e ha RSS = 200. Calcolare la statistica F classica.
Soluzione. La perdita di adattamento per restrizione è (200 − 180)/2 = 10. La varianza residua stimata nel modello non ristretto è 180/(100 − 5) = 180/95. La statistica è quindi 10/(180/95) = 95/18 ≈ 5,2778, da confrontare con una F con due e novantacinque gradi di libertà nelle condizioni classiche. Il test riguarda i due vincoli congiuntamente. Non si ottiene la stessa domanda guardando separatamente due test t. Con eteroschedasticità, la formula basata soltanto sulle RSS non è la procedura robusta generale: occorre un test di Wald con covarianza appropriata.
Esercizio 7 Varianza di un effetto combinato
Due coefficienti stimati sono β̂1 = 0,4 e β̂2 = 0,1. Le varianze sono 0,01 e 0,0025, mentre la covarianza è −0,002. Si vuole stimare θ = β1 + 2β2.
Soluzione. La stima è θ̂ = 0,4 + 2 × 0,1 = 0,6. La varianza è Var(β̂1) + 4Var(β̂2) + 4Cov(β̂1,β̂2) = 0,01 + 0,01 − 0,008 = 0,012, con errore standard circa 0,1095. Ignorare la covarianza darebbe varianza 0,02 e una misura diversa dell’incertezza. La covarianza negativa riduce qui la variabilità della combinazione, ma non lo farebbe necessariamente per una differenza anziché una somma. Ogni effetto marginale, moltiplicatore o confronto fra gruppi deve essere scritto come funzione dei coefficienti prima di costruirne l’inferenza.
Esercizio 8 Distorsione da variabile omessa
Il modello pertinente è Y = 5 + 1,5X − 2Z + u. La pendenza della proiezione di Z su X e costante è 0,4. Si stima soltanto Y su X e costante, assumendo valide le condizioni di ortogonalità del modello completo.
Soluzione. La pendenza della regressione corta è 1,5 + (−2) × 0,4 = 0,7. La distorsione rispetto al coefficiente di X a Z fissata è −0,8. Z è positivamente associata a X ma riduce Y, quindi il suo contributo omesso attenua l’associazione attribuita a X. Non si tratta di errore campionario eliminabile aumentando n. L’errore standard può diventare molto piccolo attorno a 0,7, lasciando intatta la differenza rispetto a 1,5. La formula presuppone il modello specificato: se esistono altri meccanismi di endogeneità, il confronto richiede un’analisi più ampia.
Esercizio 9 Logaritmi e variazioni finite
In un modello log-log, il coefficiente di X è due. Calcolare la variazione della componente esponenziale di Y quando X aumenta del 10%. Confrontarla con l’approssimazione elastica.
Soluzione. La variazione del logaritmo è 2ln(1,10). Il rapporto della componente esponenziale è exp[2ln(1,10)] = 1,102 = 1,21, quindi l’incremento esatto è 21%. L’approssimazione basata sull’elasticità produce 2 × 10% = 20%, adeguata soltanto come approssimazione per variazioni piccole. Per interpretare il rapporto come variazione della media aritmetica di Y servono anche condizioni sul fattore di ritrasformazione dell’errore. La precisione dell’algebra non sostituisce queste condizioni. Un coefficiente correttamente interpretato nella scala logaritmica può essere usato male quando viene trasferito senza cautela alla media in livelli.
Esercizio 10 Interazioni e differenze fra gruppi
Il modello stimato è ŷ = 30 + 4x − 6D + 3Dx. Determinare le due pendenze e la differenza fra gruppi in x = 0 e in x = 4.
Soluzione. Per D = 0, la previsione è 30 + 4x e la pendenza quattro. Per D = 1, è 24 + 7x e la pendenza sette. La differenza fra gruppi è −6 + 3x: vale −6 in x = 0 e sei in x = 4. Il coefficiente negativo della dummy non implica quindi che il gruppo D = 1 abbia sempre un esito inferiore. Il punto di incrocio è x = 2, che deve essere valutato rispetto al supporto osservato. Per verificare la differenza in x = 4 si testa −6 + 4 × 3 attraverso la combinazione dei coefficienti e la sua covarianza, non soltanto il coefficiente della dummy.
Esercizio 11 VIF e precisione
La regressione di una covariata sugli altri regressori produce Rj2 = 0,96. Calcolare il VIF e il fattore di aumento dell’errore standard nella formula classica, a parità degli altri elementi.
Soluzione. VIF = 1/(1 − 0,96) = 25. Poiché l’errore standard è la radice della varianza, il fattore corrispondente è cinque. La covariata conserva soltanto il 4% della propria variazione lineare non spiegata dagli altri regressori. Questo segnala difficoltà nel separarne il contributo, non prova che sia inutile. Eliminarla potrebbe ridurre la varianza ma introdurre distorsione se fosse un confondente necessario. Un VIF è una diagnostica dell’informazione lineare e non una regola causale di selezione. Inoltre, il confronto a parità degli altri elementi è teorico: rimuovere regressori può cambiare anche la varianza residua e l’estimando.
Esercizio 12 Ponderazione per precisione
Due misure indipendenti dello stesso parametro valgono quattro e dieci e hanno varianze note uno e nove. Calcolare la media ponderata per l’inverso della varianza e la sua varianza.
Soluzione. I pesi grezzi sono uno e un nono; normalizzandoli si ottengono 0,9 e 0,1. La stima è 0,9 × 4 + 0,1 × 10 = 4,6. La sua varianza è 0,92 × 1 + 0,12 × 9 = 0,9. La media semplice sarebbe sette e avrebbe varianza (1 + 9)/4 = 2,5. Il vantaggio nasce dal fatto che la seconda misura è molto meno precisa. Se le varianze fossero errate o i due valori misurassero parametri differenti, la giustificazione cambierebbe. I pesi di precisione non devono essere confusi con pesi di campionamento o frequenza.
Esercizio 13 Cluster e numerosità effettiva
Un campione contiene seicento osservazioni in gruppi di venticinque, con correlazione intraclasse 0,08. Usando la formula semplificata per una media con gruppi uguali, calcolare l’effetto del disegno e la numerosità effettiva approssimata.
Soluzione. DEFF = 1 + (25 − 1) × 0,08 = 2,92. La numerosità effettiva è circa 600/2,92 = 205,48 e il fattore dell’errore standard è √2,92 ≈ 1,7088 rispetto al calcolo indipendente. I gruppi sono ventiquattro, informazione importante per l’approssimazione cluster. La formula descrive un caso omogeneo e non sostituisce una covarianza appropriata per una regressione. Mostra però perché seicento righe non equivalgano sempre a seicento informazioni indipendenti. Aggiungere unità all’interno degli stessi gruppi e aggiungere nuovi gruppi possono produrre guadagni di precisione molto diversi.
Esercizio 14 IV, primo stadio e LATE
Un’offerta assegnata casualmente aumenta la partecipazione dal 25% al 65% e l’esito medio da cinquanta a cinquantasei. Calcolare ITT, primo stadio e rapporto di Wald. Precisare l’interpretazione.
Soluzione. ITT = 56 − 50 = 6. Il primo stadio è 0,65 − 0,25 = 0,40 e il rapporto è 6/0,40 = 15. L’effetto dell’offerta è sei; sotto esclusione, monotonicità, rilevanza e le altre condizioni pertinenti, quindici è l’effetto del trattamento sui compliers. Non è automaticamente l’ATE di tutta la popolazione. Se l’offerta producesse un effetto diretto sull’esito pari a uno, il rapporto incorporerebbe una componente aggiuntiva di 1/0,40 = 2,5 nel semplice modello additivo. La validità dell’esclusione è quindi decisiva anche quando la randomizzazione dell’offerta è credibile.
Esercizio 15 Sovraidentificazione GMM
Un modello stima tre parametri con cinque condizioni di momento. La statistica J efficiente vale sei. Individuare i gradi di libertà e calcolare il p-value usando il riferimento asintotico chi-quadrato.
Soluzione. I gradi di libertà sono cinque meno tre, quindi due. Per una chi-quadrato con due gradi di libertà, la probabilità di superare sei è exp(−6/2) = exp(−3) ≈ 0,0498. Il risultato si trova vicino alla soglia convenzionale del 5%. Non identifica quale momento sia eventualmente invalido e dipende dalle regolarità del test. Un numero elevato di strumenti o identificazione debole possono compromettere l’approssimazione. La soglia non deve essere interpretata come una separazione netta fra strumenti certamente validi e certamente invalidi. L’argomentazione economica sulle condizioni di momento rimane necessaria.
Esercizio 16 Odds ratio e probabilità
Un logit implica un odds ratio pari a tre per il passaggio di una dummy da zero a uno. La probabilità iniziale di un profilo è 0,20. Determinare la nuova probabilità e la variazione in punti percentuali.
Soluzione. Le odds iniziali sono 0,20/0,80 = 0,25. Moltiplicandole per tre diventano 0,75. La nuova probabilità è 0,75/(1 + 0,75) = 3/7 ≈ 0,4286. L’aumento è circa 22,86 punti percentuali, non quaranta punti e non un aumento della probabilità al 60%. Il rapporto delle probabilità è circa 2,1429, distinto dall’odds ratio tre. Con un altro livello iniziale, lo stesso coefficiente produrrebbe una variazione assoluta diversa. La comunicazione dell’effetto richiede quindi un profilo specificato o una media sulla distribuzione delle covariate.
Esercizio 17 Conteggi ed esposizione
Un modello esponenziale prevede un tasso di 1,5 eventi per mese. Un’unità è osservata per quattro mesi. Una covariata binaria ha coefficiente ln(1,2). Calcolare i conteggi attesi nelle due condizioni.
Soluzione. Senza la caratteristica, la media del conteggio è 4 × 1,5 = 6. Con la caratteristica, il tasso viene moltiplicato per 1,2 e la media diventa 7,2. Il cambiamento è del 20% nella media e di 1,2 eventi nell’intero periodo. Un’esposizione diversa cambierebbe la differenza assoluta ma non il rapporto nella specificazione adottata. Queste conclusioni sulla media valgono sotto il modello della media esponenziale. Per affermare che la probabilità di zero sia exp(−6) oppure exp(−7,2) occorre assumere anche una distribuzione Poisson, non soltanto usare PPML con covarianza robusta.
Esercizio 18 Effetto marginale tobit
In un Tobit con censura inferiore a zero, μ = σ e il coefficiente di una covariata è due. Determinare l’effetto marginale sulla media osservata, usando Φ(1) ≈ 0,8413.
Soluzione. Il rapporto z = μ/σ è uno. La probabilità di osservare un valore positivo è circa 0,8413 e l’effetto sulla media complessiva è βjΦ(z) = 2 × 0,8413 = 1,6826. Due è invece il coefficiente della variabile latente. L’effetto sulla probabilità di positività sarebbe φ(1) × 2/σ e richiederebbe il valore di σ. L’effetto sulla media dei soli positivi utilizza ancora una formula diversa. La distinzione evita di attribuire lo stesso numero a partecipazione, intensità condizionata e media totale. Le formule dipendono dalle assunzioni distributive e dalla varianza costante del Tobit classico.
Esercizio 19 Effetti fissi in due periodi
Tre unità hanno cambiamenti di x pari a uno, due e tre e cambiamenti di y pari a due, tre e cinque. Si stima un modello a effetti individuali in due periodi, senza effetto temporale comune aggiuntivo. Determinare la pendenza.
Soluzione. In due periodi il within è equivalente alla regressione in prime differenze senza intercetta per questa specificazione. La pendenza è ∑ΔxΔy/∑(Δx)2 = (1 × 2 + 2 × 3 + 3 × 5)/(1 + 4 + 9) = 23/14 ≈ 1,6429. Aggiungere un’intercetta alla regressione differenziata introdurrebbe invece un cambiamento temporale comune e produrrebbe un modello diverso. Le componenti individuali costanti sono eliminate, ma eventuali shock variabili correlati con Δx restano nel residuo. L’equivalenza algebrica non dimostra quindi esogeneità causale dei cambiamenti.
Esercizio 20 Differenze nelle differenze
La media del gruppo trattato passa da quaranta a quarantanove; quella del controllo da trenta a trentaquattro. Costruire il controfattuale del trattato e la DiD. Spiegare quale shock invaliderebbe l’interpretazione causale.
Soluzione. Il controllo cresce di quattro. Sotto tendenze parallele, la media controfattuale del trattato sarebbe 40 + 4 = 44. L’effetto stimato è 49 − 44 = 5, equivalente a (49 − 40) − (34 − 30). Se dopo l’intervento il gruppo trattato ricevesse anche un beneficio distinto che aumentasse l’esito di tre unità, la DiD potrebbe attribuirlo al trattamento studiato. Il parallelismo riguarda l’evoluzione senza trattamento, non l’uguaglianza dei livelli iniziali. La differenza iniziale di dieci non invalida da sola il disegno; uno shock differenziale concomitante può farlo.
Esercizio 21 Regressione discontinua fuzzy
Alla soglia di assegnazione, la probabilità di partecipazione salta da 0,20 a 0,70 e la media dell’esito presenta un salto di tre. Calcolare l’effetto locale IV e discutere il ruolo della banda.
Soluzione. Il salto del trattamento è 0,50 e il rapporto fuzzy è 3/0,50 = 6. Sotto le condizioni locali di continuità, esclusione e monotonicità, sei è un effetto medio per le unità il cui trattamento cambia a causa della soglia. La banda determina quante osservazioni vicine contribuiscono alla stima dei limiti: una banda piccola riduce alcune componenti di bias ma può aumentare molto la varianza. Il rapporto non è automaticamente valido per unità lontane dalla soglia. Inoltre, usare intervalli convenzionali senza considerare bias locale e incertezza del primo stadio può produrre precisione apparente non giustificata.
Esercizio 22 Durata e rischio cumulato
Un modello esponenziale ha rischio costante 0,2 per mese. Calcolare durata media, sopravvivenza a cinque mesi e probabilità di evento entro cinque mesi.
Soluzione. La durata media è 1/0,2 = 5 mesi. Il rischio cumulato a cinque mesi è 0,2 × 5 = 1. La sopravvivenza è exp(−1) ≈ 0,3679 e la probabilità di evento è 1 − exp(−1) ≈ 0,6321. Il fatto che la durata media sia cinque mesi non implica che metà delle unità abbia avuto l’evento a cinque mesi: la distribuzione è asimmetrica. La mediana è ln(2)/0,2, circa 3,4657 mesi. Il rischio 0,2 non è una probabilità mensile esatta del 20%; quella probabilità è 1 − exp(−0,2), circa 18,13%.
Esercizio 23 Previsione AR(1)
Si considera Yt = 1 + 0,5Yt-1 + εt, con innovazioni di varianza quattro e processo stazionario. Il valore corrente è sei. Determinare media incondizionata, previsione a due passi e varianza dell’errore di previsione con parametri noti.
Soluzione. La media è 1/(1 − 0,5) = 2. La previsione a due passi è 2 + 0,52(6 − 2) = 3. La varianza dell’errore è 4(1 + 0,52) = 5 e la deviazione standard √5 ≈ 2,2361. La varianza incondizionata della serie è invece 4/(1 − 0,25) = 16/3. L’incertezza della previsione cresce verso quel limite con l’orizzonte. Se i parametri fossero stimati, la formula a parametri noti non incorporerebbe tutta la loro incertezza. La previsione puntuale tre non implica che il valore futuro sarà vicino a tre con certezza.
Esercizio 24 Ardl e lungo periodo
Il modello è yt = c + 0,4yt-1 + 0,3xt + 0,2xt-1 + ut. Calcolare effetto immediato, moltiplicatore di lungo periodo e coefficiente del termine di correzione.
Soluzione. L’effetto immediato è 0,3. Il moltiplicatore di lungo periodo è (0,3 + 0,2)/(1 − 0,4) = 5/6 ≈ 0,8333. Il coefficiente del termine di correzione è 0,4 − 1 = −0,6. La relazione di equilibrio include una costante c/0,6 e pendenza 5/6. Il modello può essere riscritto con Δy, Δx e lo scostamento di livello ritardato. La stabilità e l’eventuale cointegrazione devono essere giustificate prima di interpretare economicamente il lungo periodo. Un rapporto numerico definito non dimostra da solo una relazione stabile fra livelli.
Esercizio 25 Cointegrazione
X segue un random walk e Y = 7 − 0,5X + v, con v stazionario. Stabilire se Y e X possono essere cointegrate e spiegare perché un’analisi sulle sole differenze può perdere informazione.
Soluzione. La combinazione Y − 7 + 0,5X coincide con v ed è stazionaria. Se entrambe le serie in livelli sono I(1), esiste quindi una relazione cointegrante nella struttura indicata. Una regressione soltanto delle differenze può ignorare lo scostamento dall’equilibrio precedente, che può essere utile per rappresentare l’aggiustamento. Questo non significa che ogni regressione in livelli abbia inferenza convenzionale valida: test di cointegrazione e stima del lungo periodo richiedono procedure appropriate. La cointegrazione descrive una relazione statistica stabile e non assegna automaticamente una direzione causale da X a Y.
Esercizio 26 GARCH e persistenza
Un GARCH(1,1) ha ω = 0,05, α = 0,15 e β = 0,80. La varianza corrente è due e l’innovazione corrente al quadrato è sei. Calcolare varianza incondizionata e previsione successiva.
Soluzione. La persistenza è 0,95 e, sotto le condizioni di secondo momento, la varianza incondizionata è 0,05/(1 − 0,95) = 1. La previsione successiva è 0,05 + 0,15 × 6 + 0,80 × 2 = 2,55. La previsione a due passi, usando soltanto l’informazione corrente, è 1 + 0,95(2,55 − 1) = 2,4725. Il processo può dunque avere un livello medio di varianza uno e trovarsi in una fase con varianza molto maggiore. La deviazione standard successiva è √2,55, non 2,55. La condizione di varianza finita non garantisce automaticamente l’esistenza di tutti i momenti superiori.
Esercizio 27 Media, mediana e perdita
Cinque esiti sono uno, due, tre, quattro e venti. Calcolare la previsione costante che minimizza la perdita quadratica e quella che minimizza la perdita assoluta.
Soluzione. La media è trenta diviso cinque, quindi sei, e minimizza la somma dei quadrati. La mediana è tre e minimizza la somma dei valori assoluti. L’osservazione venti sposta molto la media ma non la posizione centrale ordinata. Non esiste una previsione universalmente migliore senza specificare la perdita: sei è appropriata per costi quadratici simmetrici, tre per costi assoluti simmetrici. La robustezza della mediana rispetto a un grande valore verticale non implica robustezza a ogni forma di selezione o endogeneità. Nel modello condizionato, la stessa distinzione motiva regressione della media e regressione quantilica.
Esercizio 28 Ridge e LASSO ortogonali
In un problema standardizzato con X′X/n = I, il vettore z = X′y/n è (0,9; −0,4; 0,1). Si usa λ = 0,2 nella normalizzazione RSS/(2n). Calcolare ridge e lasso.
Soluzione. Ridge divide ciascun elemento per 1 + λ = 1,2 e produce (0,75; −0,3333; 0,0833). Lasso applica la soglia morbida e produce (0,7; −0,2; 0). La terza covariata viene eliminata da lasso, mentre ridge la mantiene ridotta. Questa differenza deriva dalla geometria delle penalità, non da un test che abbia dimostrato l’assenza di un effetto della terza variabile. Se le colonne non fossero ortogonali, non sarebbe corretto applicare indipendentemente queste formule ai coefficienti OLS. Il parametro di penalizzazione dovrebbe essere scelto con una procedura di validazione coerente con l’obiettivo.
Esercizio 29 Due modelli e due metriche
Su cinque osservazioni di test, il modello A ha errori assoluti sempre pari a due. Il modello B ha errori zero, zero, zero, zero e sei. Confrontare MAE e RMSE.
Soluzione. Per A, MAE = 2 e RMSE = 2. Per B, MAE = 6/5 = 1,2 e RMSE = √(36/5) ≈ 2,6833. B è preferibile secondo la perdita assoluta media, A secondo la perdita quadratica. Il disaccordo non è un errore: le metriche penalizzano diversamente l’unico errore grande. Occorre scegliere il criterio che rappresenta il costo della decisione. Cinque osservazioni sono inoltre poche per generalizzare l’ordinamento. La valutazione deve essere indipendente dalle scelte del modello e, nei dati temporali, rispettare l’ordine e la dipendenza degli errori.
Esercizio 30 Aggiornamento bayesiano
Una probabilità p ha prior Beta(1,1). Si osservano otto successi in dodici prove indipendenti con parametro comune. Determinare posterior, media e moda e interpretare la previsione della prova successiva.
Soluzione. La posterior è Beta(9,5). La media è 9/14 ≈ 0,6429 e la moda è 8/12 = 2/3, poiché entrambi i parametri sono maggiori di uno. La probabilità predittiva di successo nella prova successiva è la media posterior, sotto lo stesso modello di scambiabilità delle prove. La frequenza osservata è 2/3 e non coincide con la media posterior perché la prior contribuisce all’aggiornamento. Una prior diversa produrrebbe una posterior diversa, soprattutto con pochi dati. L’intervallo credibile si calcola dai quantili della Beta e si interpreta condizionatamente a dati, prior e modello, non come certezza indipendente dalle ipotesi.