## Esercizio: costruire le matrici di transizione tra scaglioni di reddito
## per le famiglie italiane tra il 1991 e il 1993 sulla base dell'indagine
## 'SHIW' della Banca d'Italia
## (grazie a Erich Battistin per aver predisposto questo dataset 15 anni fa) 

## Lettura dei dati dal formato proprietario di Stata verso R
library(foreign)
panel <- read.dta("panel.dta") # aggiungere 'path' appropriato se il file
                               # non è nella directory di lavoro

## Per ogni famiglia (nquest), l'unità di osservazione, vengono misurate
## varie caratteristiche. A noi interessa 'y': il reddito totale.

## Un po' di statistiche descrittive (ma attenzione che tre anni sono
## mischiati assieme)

summary(panel$y) 
tapply(panel$y, panel$anno, summary)

plot(density(panel[panel$anno == 1991, "y"]), col="red",
     main = "Distribuzione del reddito familiare, 1991-1995")
lines(density(panel[panel$anno == 1993, "y"]), col="blue")
lines(density(panel[panel$anno == 1995, "y"]), col="green3")

## Ora assegniamo a ciascuno il quintile di reddito di appartenenza, sulla
## base della distribuzione totale (l'intervallo è codificato con 1:5)
## Chiamiamo 'qy' la nuova variabile che restituisce il quintile di
## appartenenza di ogni famiglia

panel$qy <- unclass(cut(panel$y, breaks=c(-Inf, quantile(panel$y, 0:5/5)[-1])))

## verifichiamo che ogni quintile contenga lo stesso numero di 'casi'
table(panel$qy)

## Per migliore comprensione: questi sono gli estremi di ogni intervallo
## (in migliaia di lire)
summary(cut(panel$y, breaks=c(-Inf, quantile(panel$y, 0:5/5)[-1])))

## Ora qualche magheggio per arrivare a due variabili che, per ogni famiglia,
## ci dicono in che quintile (stato) cadeva nel 1991 (t-1), e nel 1993 (t)

## separa per anno
panel91 <- panel[panel$anno == 1991, ]
panel93 <- panel[panel$anno == 1993, ]

## merge per nquest (famiglia), per ottenere qy91 e qy93 per ogni nquest
mypan <- merge(panel91, panel93, by="nquest")

## adesso qy.x è lo scaglione 1991, qy.y è lo scaglione 1993
## rinominiamoli per chiarezza:
dimnames(mypan)[[2]][which(dimnames(mypan)[[2]]=="qy.x")] <- "qy91"
dimnames(mypan)[[2]][which(dimnames(mypan)[[2]]=="qy.y")] <- "qy93"

## costruiamo la tabella di transizione
transtab <- table(mypan[,c("qy91","qy93")])
transtab

## tassi di permanenza e transizione: dividiamo transtab per i
## totali di riga (ovvero: i totali per stato in (t-1))

apply(transtab, 1, sum)
permtrans <- transtab/apply(transtab, 1, sum)
round(permtrans, 2)

## controllo: i totali di riga devono essere 1
apply(permtrans, 1, sum)