# Generate dataset individual level attribution

# Generate 50,000 observations, which will be our 50,000 consumers (of which ~50% is an existing customer)
n=50000
set.seed(1234)
consumers <- data.frame(c(1:n), rbinom(n, 1, 0.5))

names(consumers)[1] <- paste("Consumer_ID")
names(consumers)[2] <- paste("Existing_customer")

# Generate relationship length variable, which is positive (between 1 and 60 months) for existing customers
consumers$Relation_length <- ifelse(consumers$Existing_customer==1, sample(1:60, n, replace=TRUE), 0)

# Generate customer lifetime value variable, which is positive for existing customers and positively skewed
consumers$CLV <- ifelse(consumers$Existing_customer==1, rbeta(n,2,9)*10000, 0)

# Generate a variable indicating if the consumer is in the firm ad (treatment) group or the charity ad (control) group (~80/20 probability) 
consumers$Firm_banner <- rbinom(n, 1, 0.8)

# Generate a variable if the customer has signed up for the email newsletter, which is only the case for existing customers and more likely when the CLV is higher and the customer has been longer with the firm
consumers$Email_group <- ifelse(consumers$Existing_customer==1, 
                                ifelse(((log(consumers$CLV+1)+(consumers$Relation_length/30))-sample(4:12, n, replace=TRUE))>0,1,0),0)

summary(consumers)


# Generate customer journeys non-existing (i.e. non-current) customers  ----


# Simulate discrete Markov chains according to transition matrix P
run.mc.sim <- function( P, num.iters = 51 ) {
  
  # number of possible states
  num.states <- nrow(P)
  
  # stores the states X_t through time
  states     <- numeric(num.iters)
  
  # initialize variable for first state 
  states[1]    <- 1
  
  for(t in 2:num.iters) {
    
    # probability vector to simulate next state X_{t+1}
    p  <- P[states[t-1], ]
    
    ## draw from multinomial and determine state
    states[t] <-  which(rmultinom(1, 1, p) == 1)
  }
  return(states)
}

# Split up consumers dataframe for different groups of customers
consumers_new <- consumers[ which(consumers$Existing_customer==0), ]
consumers_exist_noemail <- consumers[ which(consumers$Existing_customer==1 & consumers$Email_group==0), ]
consumers_exist_email <- consumers[ which(consumers$Existing_customer==1 & consumers$Email_group==1), ]


# Setup transition matrix for non-existing (i.e. non-current) customers
P <- t(matrix(c( 0,0.5,0.02,0.13,0.05,0.1,0,0,0.2,0,
                 0,0.48,0.03,0.1,0.09,0.05,0,0,0.15,0.1,
                 0,0.25,0.1,0.1,0.1,0.07,0,0,0.28,0.1,
                 0,0.2,0.02,0.18,0.25,0.07,0,0,0.15,0.13,
                 0,0.15,0.01,0.05,0.25,0.02,0,0,0.32,0.2,
                 0,0.2,0.01,0.03,0.1,0.24,0,0,0.2,0.22,
                 0,0,0,0,0,0,0,0,0,0,
                 0,0,0,0,0,0,0,0,0,0,
                 0,0.2,0.01,0.02,0.05,0.02,0,0,0.5,0.2,
                 0,0,0,0,0,0,0,0,0,1), nrow=10, ncol=10))

num.chains     <- nrow(consumers_new)
num.iterations <- 51

# Each column stores the sequence of states for a single chains
set.seed(1234)
chain.states  <- matrix(NA, ncol=num.chains, nrow=num.iterations)

# Simulate chains
for(c in seq_len(num.chains)){
  chain.states[,c] <- run.mc.sim(P)
}


# Merge journey data 
consumers_new_journey <- as.data.frame(lapply(consumers_new, rep, num.iterations))
consumers_new_journey <- consumers_new_journey[order(consumers_new_journey$Consumer_ID),]
consumers_new_journey$Channel <- as.vector(chain.states)
consumers_new_journey$Channel_number <- rep(0:(num.iterations-1))
consumers_new_journey <- consumers_new_journey[which(consumers_new_journey$Channel>1 & consumers_new_journey$Channel<10),]


# Create variables which summarize the journeys

consumers_new$Banner_no_click <- colSums(chain.states==2)
consumers_new$Banner_click <- colSums(chain.states==3)
consumers_new$SEA_product_click <- colSums(chain.states==4)
consumers_new$SEA_brand_click <- colSums(chain.states==5)
consumers_new$Price_comp_click <- colSums(chain.states==6)
consumers_new$Email_no_click <- colSums(chain.states==7)
consumers_new$Email_click <- colSums(chain.states==8)
consumers_new$Direct_visit <- colSums(chain.states==9)
first_observation <- aggregate(consumers_new_journey, by=list(consumers_new_journey$Consumer_ID), FUN = function(x) {  first = head(x,1) } )
consumers_new$First_channel <- first_observation$Channel
last_observation <- aggregate(consumers_new_journey, by=list(consumers_new_journey$Consumer_ID), FUN = function(x) {  last = tail(x,1) } )
consumers_new$Last_channel <- last_observation$Channel
consumers_new$Amount_touchpoints <- last_observation$Channel_number



# Generate customer journeys existing customers w/o email -----------------


# Setup transition matrix for existing customers without e-mail subscription
P <- t(matrix(c( 0,0.2,0.01,0.1,0.15,0.12,0,0,0.42,0,
                 0,0.35,0.01,0.08,0.14,0.06,0,0,0.31,0.05,
                 0,0.3,0.08,0.04,0.12,0.07,0,0,0.29,0.1,
                 0,0.2,0.02,0.18,0.25,0.07,0,0,0.15,0.13,
                 0,0.15,0.01,0.05,0.25,0.02,0,0,0.32,0.2,
                 0,0.2,0.01,0.03,0.1,0.24,0,0,0.2,0.22,
                 0,0,0,0,0,0,0,0,0,0,
                 0,0,0,0,0,0,0,0,0,0,
                 0,0.2,0.01,0.02,0.05,0.02,0,0,0.5,0.2,
                 0,0,0,0,0,0,0,0,0,1), nrow=10, ncol=10))

num.chains     <- nrow(consumers_exist_noemail)
num.iterations <- 51

# Each column stores the sequence of states for a single chains
set.seed(1234)
chain.states  <- matrix(NA, ncol=num.chains, nrow=num.iterations)

# Simulate chains
for(c in seq_len(num.chains)){
  chain.states[,c] <- run.mc.sim(P)
}


# Merge journey data 
consumers_exist_noemail_journey <- as.data.frame(lapply(consumers_exist_noemail, rep, num.iterations))
consumers_exist_noemail_journey <- consumers_exist_noemail_journey[order(consumers_exist_noemail_journey$Consumer_ID),]
consumers_exist_noemail_journey$Channel <- as.vector(chain.states)
consumers_exist_noemail_journey$Channel_number <- rep(0:(num.iterations-1))
consumers_exist_noemail_journey <- consumers_exist_noemail_journey[which(consumers_exist_noemail_journey$Channel>1 & consumers_exist_noemail_journey$Channel<10),]


# Create variables which summarize the journeys

consumers_exist_noemail$Banner_no_click <- colSums(chain.states==2)
consumers_exist_noemail$Banner_click <- colSums(chain.states==3)
consumers_exist_noemail$SEA_product_click <- colSums(chain.states==4)
consumers_exist_noemail$SEA_brand_click <- colSums(chain.states==5)
consumers_exist_noemail$Price_comp_click <- colSums(chain.states==6)
consumers_exist_noemail$Email_no_click <- colSums(chain.states==7)
consumers_exist_noemail$Email_click <- colSums(chain.states==8)
consumers_exist_noemail$Direct_visit <- colSums(chain.states==9)
first_observation <- aggregate(consumers_exist_noemail_journey, by=list(consumers_exist_noemail_journey$Consumer_ID), FUN = function(x) {  first = head(x,1) } )
consumers_exist_noemail$First_channel <- first_observation$Channel
last_observation <- aggregate(consumers_exist_noemail_journey, by=list(consumers_exist_noemail_journey$Consumer_ID), FUN = function(x) {  last = tail(x,1) } )
consumers_exist_noemail$Last_channel <- last_observation$Channel
consumers_exist_noemail$Amount_touchpoints <- last_observation$Channel_number



# Generate customer journeys existing customers with email ----------------


# Setup transition matrix for existing customers with e-mail subscription
P <- t(matrix(c( 0,0.2,0.01,0.05,0.08,0.08,0.07,0.03,0.48,0,
                 0,0.31,0.01,0.06,0.14,0.06,0.07,0.03,0.27,0.05,
                 0,0.28,0.06,0.04,0.09,0.07,0.07,0.03,0.26,0.1,
                 0,0.2,0.02,0.18,0.16,0.07,0.07,0.03,0.15,0.12,
                 0,0.15,0.01,0.05,0.16,0.02,0.07,0.03,0.32,0.19,
                 0,0.16,0.01,0.03,0.08,0.2,0.07,0.03,0.2,0.22,
                 0,0.1,0.01,0.05,0.13,0.12,0.2,0.04,0.3,0.05,
                 0,0.08,0.01,0.03,0.15,0.09,0.1,0.15,0.3,0.09,
                 0,0.2,0.01,0.02,0.05,0.02,0,0,0.5,0.2,
                 0,0,0,0,0,0,0,0,0,1), nrow=10, ncol=10))

num.chains     <- nrow(consumers_exist_email)
num.iterations <- 51

# Each column stores the sequence of states for a single chains
set.seed(1234)
chain.states  <- matrix(NA, ncol=num.chains, nrow=num.iterations)

# Simulate chains
for(c in seq_len(num.chains)){
  chain.states[,c] <- run.mc.sim(P)
}


# Merge journey data 
consumers_exist_email_journey <- as.data.frame(lapply(consumers_exist_email, rep, num.iterations))
consumers_exist_email_journey <- consumers_exist_email_journey[order(consumers_exist_email_journey$Consumer_ID),]
consumers_exist_email_journey$Channel <- as.vector(chain.states)
consumers_exist_email_journey$Channel_number <- rep(0:(num.iterations-1))
consumers_exist_email_journey <- consumers_exist_email_journey[which(consumers_exist_email_journey$Channel>1 & consumers_exist_email_journey$Channel<10),]


# Create variables which summarize the journeys

consumers_exist_email$Banner_no_click <- colSums(chain.states==2)
consumers_exist_email$Banner_click <- colSums(chain.states==3)
consumers_exist_email$SEA_product_click <- colSums(chain.states==4)
consumers_exist_email$SEA_brand_click <- colSums(chain.states==5)
consumers_exist_email$Price_comp_click <- colSums(chain.states==6)
consumers_exist_email$Email_no_click <- colSums(chain.states==7)
consumers_exist_email$Email_click <- colSums(chain.states==8)
consumers_exist_email$Direct_visit <- colSums(chain.states==9)

first_observation <- aggregate(consumers_exist_email_journey, by=list(consumers_exist_email_journey$Consumer_ID), FUN = function(x) {  first = head(x,1) } )
consumers_exist_email$First_channel <- first_observation$Channel
last_observation <- aggregate(consumers_exist_email_journey, by=list(consumers_exist_email_journey$Consumer_ID), FUN = function(x) {  last = tail(x,1) } )
consumers_exist_email$Last_channel <- last_observation$Channel
consumers_exist_email$Amount_touchpoints <- last_observation$Channel_number



# Finalize the datasets ----------------------------------------------------

# Merge datasets

consumers <- rbind(consumers_new, consumers_exist_email, consumers_exist_noemail)
consumers <- consumers[order(consumers$Consumer_ID),]

consumers_journey <- rbind(consumers_new_journey, consumers_exist_email_journey, consumers_exist_noemail_journey)
consumers_journey <- consumers_journey[order(consumers_journey$Consumer_ID),]

# Remove unused material
rm(chain.states)
rm(consumers_exist_email)
rm(consumers_exist_email_journey)
rm(consumers_exist_noemail)
rm(consumers_exist_noemail_journey)
rm(consumers_new)
rm(consumers_new_journey)
rm(first_observation)
rm(last_observation)
rm(P)
rm(c)
rm(num.chains)
rm(num.iterations)
rm(run.mc.sim)


# Label Channels

consumers$First_channel <- as.character(consumers$First_channel)

consumers$First_channel[consumers$First_channel == "2"] <- "Banner_impression"
consumers$First_channel[consumers$First_channel == "3"] <- "Banner_click"
consumers$First_channel[consumers$First_channel == "4"] <- "SEA_product_click"
consumers$First_channel[consumers$First_channel == "5"] <- "SEA_brand_click"
consumers$First_channel[consumers$First_channel == "6"] <- "Price_comp_click"
consumers$First_channel[consumers$First_channel == "7"] <- "Email_received"
consumers$First_channel[consumers$First_channel == "8"] <- "Email_click"
consumers$First_channel[consumers$First_channel == "9"] <- "Direct_visit"


consumers$Last_channel <- as.character(consumers$Last_channel)

consumers$Last_channel[consumers$Last_channel == "2"] <- "Banner_impression"
consumers$Last_channel[consumers$Last_channel == "3"] <- "Banner_click"
consumers$Last_channel[consumers$Last_channel == "4"] <- "SEA_product_click"
consumers$Last_channel[consumers$Last_channel == "5"] <- "SEA_brand_click"
consumers$Last_channel[consumers$Last_channel == "6"] <- "Price_comp_click"
consumers$Last_channel[consumers$Last_channel == "7"] <- "Email_received"
consumers$Last_channel[consumers$Last_channel == "8"] <- "Email_click"
consumers$Last_channel[consumers$Last_channel == "9"] <- "Direct_visit"

# Create dummy indicating the consumer lives in a flyer region
consumers$Flyer_region <- rbinom(n, 1, 0.5)

rm(n)


# Create purchase variable (~25% of cases)
consumers$Purchase <- -2.90 + consumers$Existing_customer*.1 + log(consumers$Relation_length+1)*.06 + log(consumers$CLV+1)*.05 + consumers$Email_group*.08 + log(consumers$Banner_no_click+1)*.04 + log(consumers$Banner_no_click+1)*consumers$Firm_banner*.3 + log(consumers$Banner_no_click+1)*consumers$Firm_banner*consumers$Existing_customer*-.2 + log(consumers$Banner_click+1)*.04 + log(consumers$Banner_click+1)*consumers$Firm_banner*.3 + log(consumers$Banner_click+1)*consumers$Firm_banner*consumers$Existing_customer*-.2 + log(consumers$SEA_product_click+1)*.05 + consumers$Flyer_region*.2 + consumers$Flyer_region* log(consumers$Banner_no_click+1)*consumers$Firm_banner*.2 + log(consumers$SEA_brand_click+1)*.1 + log(consumers$Price_comp_click+1)*.12 + log(consumers$Email_no_click+1)*.02 + log(consumers$Email_click+1)*.05 + log(consumers$Direct_visit+1)*.14 + ifelse(consumers$Last_channel=="SEA_brand_click",1,0)*.1 + ifelse(consumers$Last_channel=="Price_comp_click",1,0)*.15 + ifelse(consumers$Last_channel=="Email_click",1,0)*.5 + ifelse(consumers$Last_channel=="Direct_visit",1,0)*.1 + runif(50000) + runif(50000) + runif(50000)

consumers$Purchase <- ifelse(consumers$Purchase>0,
                             ifelse(consumers$Last_channel=="Banner_impression",0,          ifelse(consumers$Last_channel=="Email_received",0,1)),0)

# Label Channels

consumers_journey$Channel_name <- as.character(consumers_journey$Channel)

consumers_journey$Channel_name [consumers_journey$Channel_name == "2"] <- "Banner_impression"
consumers_journey$Channel_name [consumers_journey$Channel_name == "3"] <- "Banner_click"
consumers_journey$Channel_name [consumers_journey$Channel_name == "4"] <- "SEA_product_click"
consumers_journey$Channel_name [consumers_journey$Channel_name == "5"] <- "SEA_brand_click"
consumers_journey$Channel_name [consumers_journey$Channel_name == "6"] <- "Price_comp_click"
consumers_journey$Channel_name [consumers_journey$Channel_name == "7"] <- "Email_received"
consumers_journey$Channel_name [consumers_journey$Channel_name == "8"] <- "Email_click"
consumers_journey$Channel_name [consumers_journey$Channel_name == "9"] <- "Direct_visit"


# Get descriptive statistics
library(psych)
describe(consumers)

