#-----------------------------------------------------#
#                                                     #
#        LABORATORIO DI METODI STATISTICI             #
#                PER LA FINANZA                      #
#               PROF F. BARTOLUCCI                    #
#                                                     #
#                    LAB-05                          #
#-----------------------------------------------------#



##########################################################################
# Lezione 5
##########################################################################1
#                                                                           
#   Indice:                                                                 
#           
#   * Richiami di calcolo vettoriale e matriciale
#   * Regressione lineare semplice
#   * Regressione lineare semplice: forma matriciale
#   * Regressione lineare semplice: funzione lm()
#
##########################################################################


# -----------------------------------------------
# Modello di regressione lineare semplice
# -----------------------------------------------

data = read.table("regular.txt", header = TRUE) # regular.txt
head(data)

# covariata
x0 = data$reddito
# risposta
y = data$spesa 

# distribuzione di x ed y
summary(x0)
summary(y)

# scatterplot x vs y
plot(x0, y)



# stima dei parametri del modello di regressione semplice
# ---------------------------------------------------------
# beta
covXY = cov(x0,y)
varX = var(x0)
beta1hat = covXY / varX
beta1hat

# intercetta
meanY = mean(y)
meanY
meanX = mean(x0)
meanX
beta0hat = meanY - beta1hat*meanX
beta0hat

# stima dei valori di y
# -----------------------------
yhat = beta0hat + beta1hat * x0

# calcolo dei residui
# -------------------------
e = y - yhat

# stima della varianza dell'errore
# -----------------------------------
# residual sum of squares
e2 = e^2
RSS = sum(e2)

# varianza dell'errore
n = length(y)
s2 = RSS/(n-2)
s2

# varianza di beta1
# -----------------------
devX = sum((x0-meanX)^2)
varBeta1 = s2 /devX
varBeta1

# standard error di beta1
# ------------------------------
seBeta1 = sqrt(varBeta1)
seBeta1
# verifica della significativitą di beta1: H0: beta = 0, H1: beta != 0
# -------------------------------------------------------------------------
t = beta1hat / sqrt(varBeta1)
t
# calcoliamo il p-value
pValue = 2*(pt(abs(t), n-2, lower.tail = FALSE))
pValue
# R H0


# Calcolo r^2
# -------------------------------------
# Explained sum of squares
ESS = sum((yhat - meanY)^2)
# residual sum of squares
RSS = sum((y - yhat)^2)
# total sum of squares
TSS = RSS + ESS

# r^2
r2 = ESS / (ESS + RSS)
r2


# ------------------------------------------------------------
# Modello di regressione lineare semplice: forma matriciale
# ------------------------------------------------------------



# ----------------------------------------------------
# Richiami di calcolo vettoriale e matriciale
# ----------------------------------------------------

## Definizione di vettori (colonna)
a = c(-2, 3, 1, 0)
b = c(1, 3, 5, 2)

# trasposto (vettore riga)
t(b)

# estrazionedi alcuni elementi del vettore
a[2]
a[1:2]
a[c(1,3)]

# operazioni sugli elementi vettoriali
a-1
1-a
2*a 
a^2
log(b)

## Operazioni tra vettori
a+b
a-b
# prodotto '*' e divsione '/': moltiplico e divido ogni elemento di a per il corrispondente
# elemento di b; ottengo un vettore di dimensione 4
a*b
a/b
# prodotto scalare '%*%': moltiplico ogni elemento di a per il corrispondente
# elemento di b e sommo; ottengo uno scalare
a %*% b 
t(a) %*% b
# prodotto scalare '%*%': moltiplico ogni elemento di a
# per tutti gli elementi di t(b); ottengo una matrice 4x4
a %*% t(b)

# crezione di matrici
X = matrix(1:12, nrow = 4, ncol = 3)
X

# funzioni cbind e rbind
X = cbind(a, b)
X
X = rbind(a, b); X

## Operazioni con le matrici
X+1
X-1
2*X
X/10
X %*% a # Errore! Attenzione alla dimensione di X e di a
a %*% X # Ok
t(X) %*% a # Ok
Y = cbind(a, b)
X = matrix(rchisq(10, 2), nrow = 2, ncol = 5)
dim(Y)
dim(X)
YX = Y %*% X; YX

# Calcolare il determinante di una matrice
X = matrix(rnorm(9), 3, 3)
det(X)

# Calcolare l'inversa di una matrice
Xinv = solve(X)
Xinv
Xinv %*% X
round(Xinv %*% X, digits = 2)



# matrice del disegno associata al modello
x = cbind(rep(1,n),x0)

# Stima dei parametri
# ------------------------------------------------
betahat = solve(t(x) %*% x) %*% t(x) %*% y

# stima dei valori di y
# --------------------------
yhat = x %*% betahat

# stima della matrice di var. e cov. del vettore beta
# --------------------------------------------------------
V = s2 * solve(t(x)%*%x)
V
# stima degli errori standard di beta
se = sqrt(diag(V))

# statistica-test t per verificare la significativitą dei singoli parametri beta: H: beta = 0, H1: beta != 0
t = betahat/se

# calcoliamo il p-value
p_value = 2*(pt(abs(t), n-2, lower.tail = FALSE))

# creiamo una tabella riepilogativa delle stime
coeff = cbind(betahat, se, t, p_value)
coeff
# diamo un nome alle colonne
colnames(coeff) = c("Estimates", "Se", "t", "p-value")
coeff

# ---------------------------------------------------------------------
# Stima del modello di regressione lineare semplice: funzione lm()
# ---------------------------------------------------------------------
?lm

mod = lm(y ~ x0)
mod

# maggiori informazioni...
summary(mod)

# rappresentazione grafica
# --------------------------------
plot(x0, y)
abline(mod, col = "red")


## Estrazione delle informazioni dall'oggetto lm
# ------------------------------------------------
# coefficienti
coefficients(mod) # oppure:
mod$coefficients

# fitted values
fitted(mod) # oppure:
mod$fitted.values

# residui
residuals(mod) # oppure:
mod$residuals

# varianza dell'errore
summary(mod)$sigma^2

# r^2
summary(mod)$r.squared

# r^2 aggiustato che tiene conto del numero dei parametri del modello
summary(mod)$adj.r.squared


# salviamo e chiudiamo tutto
savehistory()