```{r, message = FALSE, warning = FALSE}
library(xgboost)
library(tidyverse)
load("data_mat_rk.RData")

tick <- levels(data_mat_rk$CompanyRelatedID)        # IDs of stocks
dates <- levels(as.factor(data_mat_rk$DataDate))    # List of dates
# N <- length(tick)                   

oos <- "2007-12-06"
dates_oos <- dates[dates>oos]                       # Out-of-sample dates
Tt <- length(dates_oos)                             # Number of dates
nb_port <- 6                                        # Number of portfolios
q_port <- 1 / nb_port                     
q_port <- seq(0, 1, by  = q_port)                   # Quantiles for portfolio construction
r_port <- matrix(0, ncol = nb_port, nrow = Tt)      # Portfolio returns
offset_date <- 0                                    # 0 => is equivalent to 5 years
max_depth <- 5                                      # Maximum nb of levels for the trees
q <- 0.2                                            # Filtering intensity: be careful, inverses for bulk filter

ptm <- proc.time()                                  # Initiate CPU time
for(t in 1:Tt){
    if(t%%12==0){print(dates_oos[t])}
    x <- filter(data_mat_rk, 
                DataDate < as.Date(dates_oos[t],origin="1970-01-01")-365,  # No forward-looking bias!
                DataDate > dates[t + offset_date])                         # Rolling window
    # The line below determines the filtering intensity. Omit it to remove the filter.
    x <- filter(x, R12M_F >= quantile(x$R12M_F,1-q) | R12M_F<=quantile(x$R12M_F,q))      # Filtered data: extreme
    # x <- filter(x, R12M_F >= quantile(x$R12M_F,q) & R12M_F <= quantile(x$R12M_F,1-q))  # Filtered data: bulk
    train_data <- dplyr::select(x,-CompanyRelatedID,-DataDate,-R1M_F,-R12M_F) %>% data.matrix()
    train_label <- dplyr::select(x, R12M_F) %>% data.matrix()                            # Training label
    train_matrix <- xgb.DMatrix(data = train_data, label = train_label)                  # Full training data
    fit <- xgb.train(eta = 0.7,                    # Learning rate, 0.7 for light, 0.5 for intermed. 0.3 for deep
                     max_depth = max_depth,        # Maximum depth
                     gamma = 0,                    # Penalization
                     lambda = 1,                   # Penalization
                     objective = "reg:linear",     # Objective function
                     data = train_matrix,          # Data source
                     nrounds = 25)                 # Nb of trees, 25 for light, 50 for intermed., 100 for deep
    y <- filter(data_mat_rk, DataDate == dates_oos[t]) # Current values
    test_data <- dplyr::select(y,-CompanyRelatedID,-DataDate,-R1M_F,-R12M_F) %>% data.matrix()
    test_label <- dplyr::select(y, R1M_F) %>% data.matrix()          
    test_matrix <- xgb.DMatrix(data = test_data, label = test_label) # Full testing data
    pred <- predict(fit, test_matrix)                                # Prediction
    z <- filter(data_mat_rk, DataDate == dates_oos[t])               # Current data
    for(j in 2:length(q_port)){ # Loop creating the portfolios (could use quantcut)
        ind <- which(pred <= quantile(pred, q_port[j]) & pred >= quantile(pred,q_port[j-1]))
        r_port[t,j-1] <- mean(z$R1M_F[ind]) # Aggregating stock returns into portfolios
    }
}
proc.time() - ptm                           # CPU time difference

colMeans(r_port)                            # Average of portfolio returns
t.test(r_port[,6]-r_port[,1])               # t-test on difference of extreme portfolios

ew <- data_mat_rk %>%                       # Benchmark: EW portfolio
    group_by(DataDate) %>%
    summarise(Ret = mean(R1M_F))
mean(ew$Ret)                                # Average return of benchmark
```


