##A reproducible R session illustrating the packages and functions employed in the analyses featured in "Identifying the Relationship between Earnings and Prices" by C. Starica and J. Marton.
##The read.me file outlines the specific R version and key packages utilized in the analyses.

##Upload the needed libraries (supposed to be installed)
library(randomForest)
library(ranger)
library(gbm)
library(caret)
library(ggplot2)
library(pdp)
library(iml)

##Read the data
DAT = read.csv('Data file.csv', header = TRUE)
head(DAT)
##The data consists of 1000 firm-year observations of market value (dep.var), earnings(ib), price-to-book, leverage and the 48 Fama-French industry classification - see the article for the definition of the variables
dim(DAT)
names(DAT)

## Infer the conditional mean on the dataset. We can use different packages:
##with the package 'randomForest'
rf_model_rf <- randomForest(dep.var ~ ., data = DAT)

##with the package 'ranger'
rf_model_ranger <- ranger(dep.var ~., data = DAT)

##with the package 'gbm'
rf_model_gbm <- gbm(dep.var ~., distribution = 'gaussian', cv.folds = 5, interaction.depth = 4, data = DAT)

## Infer conditional quantiles on the dataset
##with the package 'gbm'
rf_model_q_gbm <- gbm(dep.var ~., distribution = list(name = "quantile", alpha = 0.95), cv.folds = 5, interaction.depth = 4, data = DAT)

#######################################
##Use package 'caret' to calculate Individual Conditional Expectations curve and to visualize the Partial Dependence function (the construct in equation (5) in the paper):

##A Predictor object is a container for the prediction model and the data. This ensures that the machine learning model can be analyzed in a robust way
mod <- Predictor$new(rf_model_rf, data = DAT)

## Compute the accumulated local effects for the first feature (earnings) with a partial dependence plot and ice curves
eff <- FeatureEffect$new(mod,
                         feature = "earnings", method = "pdp+ice",
                         grid.size = 100
                         )
##Visualization
plot(eff) +
    ggtitle("Partial dependence")

########################################
##Use package 'pdp' to visualize the ICE curves and the Partial Dependence function:
# Generate ICE curves for a feature (e.g., dep.var)
ice_data <- partial(rf_model_rf, pred.var = "earnings", ice = TRUE)

# Plot the ICE curves
plotPartial(ice_data)

# Plot ICE curves with color and PDP overlay
plotPartial(ice_data, center = FALSE, plot.pdp = TRUE, 
            pdp.col = "yellow", smooth = TRUE, ylab = 'Price', xlab = 'Earnings')

####################
##Two-dimensional Partial Dependence function (might take a while to run)
ice_data <- partial(rf_model_rf, pred.var = c("earnings", 'pb'))

# Plot the heat map of the PD surface
plotPartial(ice_data)