#--- 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 <- 100 # numero simulazioni
n <- 50 # numero osservazioni
p <- 20 # 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

# mse_ols_list <- list()
# mse_ridge_list <- list()
# lambda_list <- list()
# correlation_list <- list()
df <- data.frame(correlation = numeric(0), mse_ols = numeric(0), mse_ridge = numeric(0), optimal_lambda = numeric(0))

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

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

  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
  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(correlation=mean_corr, 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$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_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()