## import data
dati <- read.csv("DatiBoldrin.csv", header=TRUE, sep="\t")

## statistiche descrittive ROA

## tutta la distribuzione
summary(dati$ROA...2018)

## solo falliti
summary(dati[dati$Output == 1, "ROA...2018"])

## solo sopravvissuti
summary(dati[dati$Output == 0, "ROA...2018"])

## plot densità
plot(density(dati[dati$Output == 1, "ROA...2018"]), type="l",
     col="red", xlim=summary(dati$ROA...2018)[c(1,6)],
     ylim=c(0, 0.1))
lines(density(dati[dati$Output == 0, "ROA...2018"]), type="l",
     col="blue")

## se uno volesse aggiungere la distribuzione di tutti:
lines(density(dati[, "ROA...2018"]), lty=2,
      col="green3")

## PN/Attivo 

## plot dist
plot(density(dati[dati$Output == 1, "PN.Attivo...2018"]), type="l",
     col="red", xlim=c(-500,200),
     ylim=c(0, 0.03))
lines(density(dati[dati$Output == 0, "PN.Attivo...2018"]), type="l",
      col="blue")
abline(v=0, lty=2)

## se uno volesse aggiungere la distribuzione di tutti:
lines(density(dati[, "PN.Attivo...2018"]), lty=2,
      col="green3")

####### analisi discriminante

library(MASS)
altman <- lda(formula = Output ~ ROA...2018 + CCN.Attivo...2018 +
                  PN.Attivo...2018,
               data=dati)

## valori previsti (qui, in sample), prime dieci
predict(altman)$class[1:10]

## probabilità di fallimento previste (idem)
predict(altman)$posterior[1:10, ]

## per confronto: modello logit:

## logit
logit1 <- glm(formula = Output ~ ROA...2018 + CCN.Attivo...2018 +
                 PN.Attivo...2018,
              data=dati, family = binomial(link="logit"))
summary(logit1)

############### usando tutte le variabili:

## (modo automatico di costruire la formula)
fm <-as.formula(
   paste("Output ~",
         paste(dimnames(dati)[[2]][-c(1,13)], collapse=" + ")))

## Linear probability model (LPM)
lp.mod <- lm(fm, dati)
summary(lp.mod)

## logit
lgt.mod <- glm(fm, data=dati, family=binomial(link="logit"))
summary(lgt.mod)

## probit
pbt.mod <- glm(formula=fm, data=dati, family=binomial(link="probit"))

## mettiamoli a confronto:
library(texreg)
screenreg(list(LPM=lp.mod, Logit=lgt.mod, Probit=pbt.mod))

############ valutazione dell'efficacia previsiva:

## test and training datasets
training <- dati[c(1:100, 186:285), ]
test <- dati[-c(1:100, 186:285), ]

lgt.1 <- glm(fm, data=training, family=binomial(link="logit"))

## make test y.hat
y.hat.lgt <- predict(lgt.1, newdata=test, type="response")

y.hat.lgt.bin <- round(y.hat.lgt, 0)

## compare to reality:
table(data.frame(predicted=y.hat.lgt.bin, actual=test$Output))

####### LDA, same proc. #########
library(MASS)
lda.1 <- lda(fm, data=training)

## make test y.hat
y.hat.lda <- predict(lda.1, newdata=test)

## compare to reality:
table(data.frame(predicted=y.hat.lda$class, actual=test$Output))