#####################################################################################################
##
##                                    Misspecification
##
## E. Flachaire, 16/12/2021
##
#####################################################################################################


###### Application 1: 

library(MASS);library(randomForest);library(gbm);library(glmnet)
data(Boston) ; nobs=nrow(Boston)
set.seed(12345) ; nfold=10 
Kfold=cut(seq(1,nobs),breaks=nfold,labels=FALSE)
mse.test=matrix(0,nfold,4)
# generate X^2 X^3 and pairwise interactions for the Lasso
Xcol=colnames(Boston)[-14]
Xsqr=paste0("I(",Xcol,"^2)", collapse="+")  # squared covariates
Xcub=paste0("I(",Xcol,"^3)", collapse="+")  # cubic covariates
fmla=paste0("medv~(.)^2+",Xsqr,"+",Xcub)
Z=model.matrix(as.formula(fmla), data=Boston)[,-1]
y=Boston[,14]
#
mysample=sample(1:nobs)   # random sampling (permutation)
for(i in 1:nfold){  # K-fold CV
  cat("K-fold loop: ", i, "\r")  
  test=mysample[which(Kfold==i)]
  train=mysample[which(Kfold!=i)]
  # OLS, Lasso, Random Forest, Boosting  
  fit.lm <- lm(medv~., data=Boston, subset=train)
  fit.la <- cv.glmnet(Z[train,], y[train], alpha=1)
  fit.rf <- randomForest(medv~., data=Boston, subset=train, mtry=6)
  fit.bo <- gbm(medv~., data=Boston[train,], distribution="gaussian", interaction.depth=6)
  # out-sample MSE
  mse.test[i,1]=mean((Boston$medv-predict(fit.lm,Boston))[-train]^2)
  mse.test[i,2]=mean((y-predict(fit.la,Z,s="lambda.min"))[-train]^2)
  mse.test[i,3]=mean((Boston$medv-predict(fit.rf,Boston))[-train]^2)
  mse.test[i,4]=mean((Boston$medv-predict(fit.bo,Boston))[-train]^2)
}
mse=colMeans(mse.test)  # test error
round(mse, digits=2)