---
title: "Esercizio15052020"
author: "499EC"
date: "15/5/2020"
output:
  html_document:
    df_print: paged
  pdf_document: default
  word_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```

## Esercizio di analisi in componenti principali

Svolgiamo l'esercizio di *Cap8.R* inserendolo in un documento di R markdown.

### I dati

Prendiamo i dati della *Tabella 8.7* nel Biggeri *et al.*.

```{r import}
Tab8.7 <- read.table("Tab8.7.txt", sep="\t", h=T, dec=".")
Imprese <- Tab8.7[,-1]    # via l'indice della prima colonna
dimnames(Imprese)[[1]] <- c("Alfa", paste("U", 2:6, sep=""))
```

### Statistiche descrittive

Calcoliamo media e varianza

```{r pressure, echo=TRUE}

## calc. media e varianza
apply(Imprese, 2, mean)
apply(Imprese, 2, var) * (5/6) # aggiusto per var. campionaria  
```

### Analisi in component principali

Calc. le CP come da esempio del libro: (NB qui le variabili non vengono standardizzate!)

#### Scomposizione della matrice di covarianza

Scomponiamo la matrice di covarianza usando una eigendecomposition (v. eq.(8)):

```{r eigendecomp}
eigen.decomposition <- eigen(cov(Imprese))
autovalori <- eigen.decomposition$values
autovettori <- eigen.decomposition$vectors
```

I coefficienti di C1 e C2 saranno gli autovettori della matrice di covarianza delle variabili originarie: 

```{r }
## coefficienti di C1 in (13):
autovettori[,1]
## coefficienti di C2:
autovettori[,2]
```

### Interpretazione delle componenti principali

L'interpretazione delle CP (si veda p.358 Biggeri, in alto): C1 è positivamente correlata al current ratio, negativamente al rapporto di indebitamento; pertanto **è da considerarsi un indicatore di stabilità finanziaria** (più alto è, meglio è); C2 invece è residuale e non ha una chiara interpretazione.

#### Varianze delle CP

Calcoliamo le varianze delle 2 componenti principali C1, C2:


```{r }
autovalori
```

Controlliamo: la varianza totale deve restare invariata

```{r }
sum(autovalori)
sum(apply(Imprese, 2, var))
```

## Interpretazione grafica delle CP

Calcoliamo i punteggi delle componenti principali (scores) per ogni impresa (calcolati sulle variabili in deviazioni dalla propria media):

```{r }
Imprese.dev <- apply(Imprese, 2, FUN=function(x) x-mean(x))
Imprese.scores <- as.matrix(Imprese.dev) %*% autovettori
dimnames(Imprese.scores)[[2]] <- c("C1", "C2")
```

Con questi riproduciamo la *Tabella 8.8* del libro:

```{r }
round(Imprese.scores[,"C1"], 3)
```

Rappresentiamo ora le imprese nello spazio delle variabili originarie:



```{r }
plot(Imprese$BETA, Imprese$CR, xlim=c(0,3), ylim=c(0,2),
     pch=19)
text(Imprese$BETA, Imprese$CR + 0.1, row.names(Imprese))
```

e in quello definito dalle componenti principali:

```{r }
plot(Imprese.scores[,"C1"], Imprese.scores[,"C2"],
     xlim=c(-1.4,1.4), ylim=c(-1.4,1.4), pch=19)
text(Imprese.scores[,"C1"], Imprese.scores[,"C2"] + 0.15,
     row.names(Imprese))
abline(h=0, lty=2)
abline(v=0, lty=2)
## posizione delle imprese sulla prima CP:
points(Imprese.scores[,"C1"], rep(0, dim(Imprese)[[1]]),
       pch=19, col="orange")
```

## Illustrazione grafica della rotazione degli assi

Facciamo vedere le imprese nello spazio delle variabili originarie, assieme agli assi delle nuove variabili: le componenti principali (C1 in rosso e C2 in verde).

```{r }
## Figura 8.2 (a bis): imprese nello spazio scalato e
## assi delle componenti principali (vedi Fig. 8.1 (a))
plot(Imprese$BETA-mean(Imprese$BETA),
     Imprese$CR-mean(Imprese$CR),
     xlim=c(-1.5,1.5), ylim=c(-1,1), pch=19,
     main="Rotazione degli assi fattoriali")
text(Imprese$BETA-mean(Imprese$BETA),
     Imprese$CR-mean(Imprese$CR) + 0.1,
     row.names(Imprese))
abline(h=0, col=grey(0.8))
abline(v=0, col=grey(0.8))

## aggiungiamo gli assi delle CP:
## passano per l'origine (a = 0) e la pendenza delle rette
## è data da b = a21/a11 e, risp., b = a22/a12 in (13)
abline(c(0, autovettori[2,1]/autovettori[1,1]),
       lty=2, col="red")
abline(c(0, autovettori[2,2]/autovettori[1,2]),
       lty=2, col="green3")
text(1.4, -0.8, "C1", col="red")
text(0.95, 0.95, "C2", col="green3")

```