Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #--- Motecarlo ---#
- library(glmnet)
- library(mvtnorm)
- library(Matrix)
- # Functions
- high_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 <- 5 # 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
- # 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 <- high_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%*%beta_vector + epsilon)
- 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)
- 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)
- # Lists update
- # mse_ols_list[[i]] <- mse_ols
- # mse_ridge_list[[i]] <- mse_ridge
- # lambda_list[[i]] <- optimal_lambda
- # correlation_list[[i]] <- mean_corr
- 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)
- ggplot(df, aes(x = df$mean_correlation)) +
- geom_line(aes(y = mse_ols), color = "red") +
- geom_line(aes(y = mse_ridge), color = "blue")
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement