SNR <- c(1, 2, 4)
gam <- 2^seq(0.5, 5, len=19)

mMSE <- sMSE <- array(NA, dim=c(N, length(SNR), length(gam)))
pb <- txtProgressBar(1, N, style=3)
for (i in 1:N) {
  for (j in 1:length(SNR)) {
    Data <- genData(n, p, J0=p0, SNR=SNR[j])
    pData <- genData(n, p, J0=p0, SNR=SNR[j])
    for (k in 1:length(gam)) {
      fit <- with(Data, ncvreg(X, y, gamma=gam[k]))
      pred <- predict(fit, pData$X)
      lam <- fit$lambda[which.min(apply(pred-pData$y, 2, crossprod))]
      mMSE[i,j,k] <- crossprod(Data$beta - coef(fit, lambda=lam)[-1])

      fit <- with(Data, ncvreg(X, y, gamma=gam[k]+1, penalty='SCAD'))
      pred <- predict(fit, pData$X)
      lam <- fit$lambda[which.min(apply(pred-pData$y, 2, crossprod))]
      sMSE[i,j,k] <- crossprod(Data$beta - coef(fit, lambda=lam)[-1])
    }
  }
  setTxtProgressBar(pb, i)
}
save(gam, SNR, mMSE, sMSE, file='sim2.RData')