library(glmnet)
library(mvtnorm)

# Data generating process
nsim <- 1000 # numero simulazioni
n <- 50 # numero osservazioni
p <- 100 # numero varaibili indipendenti
significative <- 20 # numero variabili indipendenti significative
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
Sigma <- clusterGeneration::rcorrmatrix(p, alphad=1) # matrice di varianza e covarianza che rappresenta appunto la correlazione tra le p variabili indipendenti
mean_corr <- mean(Sigma[upper.tri(Sigma)]) # correlazione media delle variabili

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)
stime_ols <- fit.ols$coefficients

# Scelgo il lambda ottimale mediante CV
cv.ridge <- cv.glmnet(x=X,y=Y,family="gaussian", alpha=0)
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 - beta_vector)^2)
mse_ridge <- mean((stime_ridge - beta_vector)^2)


# Plots