#--- Motecarlo ---#
library(glmnet)
library(mvtnorm)
library(Matrix)
# Functions

correlation_matrix <- function(p) {
  a = runif(1,0.1,0.9)
  b = runif(1,0.1,0.9)
  corr_min = min(a,b)
  corr_max = max(a,b)

  Sigma <- diag(1, p)
  Sigma[lower.tri(Sigma)] <- runif((p^2-p)/2, corr_min, corr_max)
  Sigma[upper.tri(Sigma)] <- t(Sigma)[upper.tri(Sigma)]
  Sigma
  Sigma <- nearPD(Sigma, corr = TRUE, keepDiag = TRUE) # trova la matrice definita positiva più vicina alla matrice di input utilizzando l’algoritmo di Higham
  Sigma$mat
  return(list(matrix = Sigma$mat, mean_correlation = mean(Sigma$mat[upper.tri(Sigma$mat)]))) # restituisce matrice e correlazione media delle variabili
}


# Data generating process
nsim <- 500 # numero simulazioni
n <- 200 # numero osservazioni
p <- 100 # numero varaibili indipendenti
significative <- 20 # numero variabili indipendenti significative
beta_vector <- runif(p,-5,5)
# beta_value <- 8 # valore del coeff beta
# beta_vector <- c(rep(beta_value, significative), rep(0,p-significative)) # vettore dei beta
sigma_X <- (1/n)^.5 # deviazione standard della distribuzione normale utilizzata per generare le variabili indipendenti
sigma_epsilon <- 1 # deviazione standard della distribuzione normale utilizzata per generare l'errore

df <- data.frame(mean_sigma_inv = numeric(0), mean_correlation = numeric(0), determinant = numeric(0), mse_ols = numeric(0), mse_ridge = numeric(0), optimal_lambda = numeric(0))

for (i in 1:nsim){
  set.seed(i+1)

  # Sigma <- clusterGeneration::rcorrmatrix(p, alphad=1) # matrice di varianza e covarianza che rappresenta appunto la correlazione tra le p variabili indipendenti
  res <- correlation_matrix(p)
  Sigma <- as.matrix(res[1]$matrix)
  mean_corr <- res[2]
  Sigma
  mean_corr

  # Multicollinearità
  determinant <- det(Sigma) # http://www.stat.unipg.it/~bart/metodi/lezione14.pdf Quando c’è una forte correlazione tra due o più variabili esplicative in un modello di regressione, il determinante della matrice di varianza e covarianza può diventare molto piccolo

  Sigma_inv <- solve(Sigma) # Matrice inversa
  mean_sigma_inv <- mean(Sigma_inv[upper.tri(Sigma_inv)])
  mean_sigma_inv

  epsilon <- rnorm(n,0,sigma_epsilon) # genero l'errore da una normale N(0,1)
  X <- mvtnorm::rmvnorm(n,sigma = Sigma) # genero la matrice di dati nxp, ogni riga rappresenta un'osservazione, ogni colonna una variabile indipendente
  Y <- X%*%beta_vector + epsilon # genero la variabile dipendente come combinazione lineare delle X e del vettore dei beta con l'aggiunta di un errore

  # OLS Estimator
  fit.ols <- lm(Y ~ X)
  summary(fit.ols)
  stime_ols <- fit.ols$coefficients
  stime_ols

  # Scelgo il lambda ottimale mediante CV
  cv.ridge <- cv.glmnet(x=X,y=Y,family="gaussian", alpha=0)
  # cv.ridge <- cv.glmnet(x[train,],y[train], nfolds = 10, type.measure = "mse",alpha=0)
  # plot(cv.ridge)
  optimal_lambda <- cv.ridge$lambda.min
  optimal_lambda

  # Ridge Estimator
  # Prova senza family
  fit.ridge <- glmnet(x=X,y=Y,family="gaussian", alpha=0, lambda = optimal_lambda)
  stime_ridge <- fit.ridge$beta

  # MSE
  mse_ols <- mean((stime_ols[-1] - beta_vector)^2) # il primo valore in stime_ols è l'intercetta (?)

  mse_ridge <- mean((stime_ridge - beta_vector)^2)

  new_row <- data.frame(mean_sigma_inv=mean_sigma_inv,mean_correlation=mean_corr, determinant = determinant, mse_ols = mse_ols, mse_ridge = mse_ridge, optimal_lambda=optimal_lambda)
  df <- rbind(df, new_row)
}


library(ggplot2)
plot(df$mean_correlation, df$mse_ols)
plot(df$mean_correlation, df$mse_ridge)

ggplot(df, aes(x = df$determinant)) +
  geom_point(aes(y = df$determinant), color = "red", alpha=0.7) +
  scale_y_log10()

ggplot(df, aes(x = df$mean_correlation)) +
  geom_point(aes(y = df$mse_ols), color = "red", alpha=0.7) +
  geom_point(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
  geom_smooth(aes(y = df$mse_ols), color = "red", alpha=0.7) +
  geom_smooth(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
  scale_y_log10()


ggplot(df, aes(x = df$mean_sigma_inv)) +
  geom_point(aes(y = df$mse_ols), color = "red", alpha=0.7) +
  geom_point(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
  geom_smooth(aes(y = df$mse_ols), color = "red", alpha=0.7) +
  geom_smooth(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
  scale_y_log10()

ggplot(df, aes(x = df$mean_sigma_inv)) +
  geom_point(aes(y = df$mse_ols), color = "red", alpha=0.7) +
  geom_point(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
  geom_smooth(aes(y = df$mse_ols), color = "red", alpha=0.7) +
  geom_smooth(aes(y = df$mse_ridge), color = "blue",alpha=0.3)

ggplot(df, aes(x = df$determinant)) +
  geom_point(aes(y = df$mse_ols), color = "red", alpha=0.7) +
  geom_point(aes(y = df$mse_ridge), color = "blue",alpha=0.3) +
  geom_smooth(aes(y = df$mse_ols), color = "red", alpha=0.7) +
  geom_smooth(aes(y = df$mse_ridge), color = "blue",alpha=0.3)

ggplot(df, aes(x = df$mean_correlation)) +
  geom_point(aes(y = df$optimal_lambda), color = "red", alpha=0.7) +
  geom_smooth(aes(y = df$optimal_lambda), color = "red",alpha=0.3) +
  scale_y_log10()

ggplot(df, aes(x = df$mean_sigma_inv)) +
  geom_point(aes(y = df$optimal_lambda), color = "red", alpha=0.7) +
  geom_smooth(aes(y = df$optimal_lambda), color = "red",alpha=0.3) +
  scale_y_log10()


# La multicollinearità va guardata sull'inversa
# Posso calcolare il determinante della matrie inversa per capire il livello di correlazione




# --- Scenario 2 --- #
library(glmnet)
library(mvtnorm)
library(Matrix) # Per trovare la semidefinita positiva piu vicina
# Functions

correlation_matrix <- function(p) {
  a = runif(1,0.1,0.9)
  b = runif(1,0.1,0.9)
  corr_min = min(a,b)
  corr_max = max(a,b)

  Sigma <- diag(1, p)
  Sigma[lower.tri(Sigma)] <- runif((p^2-p)/2, corr_min, corr_max)
  Sigma[upper.tri(Sigma)] <- t(Sigma)[upper.tri(Sigma)]
  Sigma
  Sigma <- nearPD(Sigma, corr = TRUE, keepDiag = TRUE) # trova la matrice definita positiva più vicina alla matrice di input utilizzando l’algoritmo di Higham
  Sigma$mat
  return(Sigma$mat) # restituisce matrice e correlazione media delle variabili
}

# Data generating process
nsim <- 500 # numero simulazioni
n <- 100 # numero osservazioni
p <- 200 # numero varaibili indipendenti
significative <- 20 # numero variabili indipendenti significative
beta_vector <- runif(p,-5,5)
# beta_value <- 8 # valore del coeff beta
# beta_vector <- c(rep(beta_value, significative), rep(0,p-significative)) # vettore dei beta
sigma_X <- (1/n)^.5 # deviazione standard della distribuzione normale utilizzata per generare le variabili indipendenti
sigma_epsilon <- 1 # deviazione standard della distribuzione normale utilizzata per generare l'errore

df <- data.frame(mse_ridge=numeric(0),mse_lasso=numeric(0))

for (i in 1:nsim){
  set.seed(i+1)

  # Sigma <- clusterGeneration::rcorrmatrix(p, alphad=1) # matrice di varianza e covarianza che rappresenta appunto la correlazione tra le p variabili indipendenti
  Sigma <- correlation_matrix(p)
  Sigma <- as.matrix(Sigma)
  Sigma

  epsilon <- rnorm(n,0,sigma_epsilon) # genero l'errore da una normale N(0,1)
  X <- mvtnorm::rmvnorm(n,sigma = Sigma) # genero la matrice di dati nxp, ogni riga rappresenta un'osservazione, ogni colonna una variabile indipendente
  Y <- X%*%beta_vector + epsilon # genero la variabile dipendente come combinazione lineare delle X e del vettore dei beta con l'aggiunta di un errore

  # Ridge Estimator
  fit.ridge <- glmnet(x=X,y=Y, alpha=0)
  stime_ridge <- fit.ridge$beta

  # Lasso Estimator
  fit.lasso <- glmnet(x=X,y=Y, alpha=1)
  stime_lasso <- fit.lasso$beta
  stime_lasso


  # MSE
  mse_ridge <- mean((stime_ridge - beta_vector)^2)
  mse_lasso <- mean((stime_lasso - beta_vector)^2)
  new_row <- data.frame(mse_ridge=mse_ridge,mse_lasso=mse_lasso)
  df <- rbind(df, new_row)
}

mean_mse_ridge <- mean(df$mse_ridge)
mean_mse_lasso <- mean(df$mse_lasso)
mean_mse_ridge
mean_mse_lasso