## Esercizio di regressione lineare su dati cross-section

## importiamo un dataset panel tratto dall'indagine "SHIW"
## della Banca d'Italia sui consumi delle famiglie (dal
## formato originale Stata); contiene i risultati dei 
## questionari per circa 7000 famiglie osservate negli
## anni 1991, 1993, 1995

library(foreign)
dati <- read.dta("panel.dta")

## estraiamo la cross-section relativa ai soli dati 1991
## (sovrascrivendo il data.frame originale)
dati <- dati[dati$anno==1991, ]

## siamo interessati a modellizzare il consumo (c) in funzione 
## del reddito (y)

## eliminiamo i record (le righe) dove il reddito è zero
dati <- dati[!(dati$y<=0), ]

## analisi esplorativa del legame tra c e y
plot(dati$y, dati$c)
## idem in semitrasparenza, per evidenziare dove si concentra
## la massa dei dati
library(scales)
plot(dati$y, dati$c, col=alpha("red", 0.1))

## una trasformazione logaritmica aiuta a "schiacciare"
## le osservazioni eccezionalmente grandi
plot(log(dati$y), log(dati$c), col=alpha("red", 0.4))
rug(log(dati$c), col="red")

## modello: c = a + b*y + u
## (scaliamo c e y in migliaia per migliorare la leggibilità)
mod <- lm(I(c/1000)~I(y/1000), dati)
summary(mod)

## modello in logaritmi: ln(c) = a + b*ln(y) + u
## NB: in questo caso la relazione originaria è assunta
## essere moltiplicativa! c = e^a * y^b * e^u
lmod <- lm(log(c) ~ log(y), dati)
abline(lmod, col="green3", lwd=3)

summary(lmod)

## scatterplot dei residui, ordinati per reddito!
## modello ai livelli 
plot(resid(mod)[order(dati$y)], col=alpha("red", 0.3))
abline(h=0)
## risultano evidenti eteroschedasticità e problemi alla
## forma funzionale

## modello ai logaritmi
plot(resid(lmod)[order(dati$y)], col=alpha("red", 0.3))
abline(h=0)

## ordinando per consumo anziché per reddito si evidenziano
## problemi più seri di non linearità
plot(resid(lmod)[order(dati$c)], col=alpha("red", 0.3))
abline(h=0)

## (il modello ai livelli sarebbe un disastro)
plot(resid(mod)[order(dati$c)], col=alpha("red", 0.3))
abline(h=0)