-
Notifications
You must be signed in to change notification settings - Fork 0
/
IBMScenario2.R
353 lines (254 loc) · 15.3 KB
/
IBMScenario2.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
# IBM scenario 2: Licensing ####
# Author: Francesca Mancini
# Date created: 2018-06-25
# Date modified: 2018-08-14
library(doParallel)
library(foreach)
# initialise population of tourists and tour operators ####
set.seed(123)
# create dataframe of parameters
params <- data.frame(trends = rep(c(0.005, 0, -0.005), each = 15),
min_rating_shape1 = rep(c(0.8, 3, 1.9), each = 5),
min_rating_shape2 = rep(c(3, 0.8, 1.9), each = 5))
cl <- makeCluster(15)
registerDoParallel(cl)
results <- foreach(i = 1:45, .packages = c("dplyr", "RGeode")) %dopar% {
source("IBMfunctions.R")
# years and days
years <- 50
days <- 365
# fine
fine <- 1000
detection_prob <- 0.3
# create vector of behavioural phenotypes for tour operators
phenotypes <- c("trustful", "optimist", "pessimist", "envious", "undefined")
# the phenotypes will be sampled from a population with the following proportions of phenotypes
# Poncela-Casasnovas et al, 2016. Science Advances 2(8): e1600451-e1600451. doi:10.1126/sciadv.1600451
p_trustful <- 0.17
p_optimist <- 0.2
p_pessimist <- 0.21
p_envious <- 0.3
p_undefined <- 0.12
tour_ops <- data.frame(id = seq(1, 10, 1), price = rnorm(10, 30, 1), rating = rep(3, 10),
capacity = as.integer(runif(10, 10, 30)), bookings = rep(0, 10), bookings_year = rep(0, 10),
investment_infra = rep(0.001, 10), investment_ot = rep(0.001, 10),
time_with = rep(0, 10), time_with_year = rep(0, 10),
profit = rep(0, 10), profit_year = rep(0, 10),
phenotype = sample(phenotypes, 10, replace=TRUE,
prob=c(p_trustful, p_optimist, p_pessimist, p_envious, p_undefined)),
behaviour = character(10), tours = rep(365, 10), stringsAsFactors = FALSE)
tour_ops$behaviour <- factor(tour_ops$behaviour, levels = c("defect", "cooperate", "no choice"))
capacity_0 <- sum(tour_ops$capacity)
# generate daily time series of tourists
tourists_start_mean <- 200
days_tot <- years* days # number of days in total
#effect sizes
y_effect <- params[i, "trends"] # trends in demand
eff.season <- 50 # seasonal fluctuation
season.sd <- 15
# sampling noise
sampling.sd <-10
# creates a vector holding the day of the year for all of the days
season <- rep(1:days, length = days_tot)
# creates a vector holding the days
days_all <- 1:days_tot
# calculate the number of tourists as a linear function of the annual trend
# and a sinusoidal function of the day of year
# plus some noise
n_tourists <- round(rnorm(days_tot,
mean = tourists_start_mean + y_effect *days_all - rnorm(days_tot, eff.season, season.sd) * cos(2 * pi/days*season),
sd =sampling.sd))
n_tourists <- ifelse(n_tourists > 0, n_tourists, 0)
# create vectors to store wildlife and management parameters
effects <- rep(NA, years)
encounter_probs <- rep(NA, years)
encounter_probs[1] <- 0.5
max_times <- rep(NA, years)
max_times[1] <- 100000
profits <- vector("list", years)
investments <- vector("list", years)
bookings_year <- vector("list", years)
ratings_year <- vector("list", years)
prices <- vector("list", years)
withanimals <- vector("list", years)
behaviours_year <- vector("list", years)
for(y in 1:years){ # start year loop
# create year tourists population
tourists_pop <- data.frame(id = seq(1, 1000000, 1), price_max = c(rnorm(60000, 30, 1.5), rnorm(30000, 45, 3.5), rnorm(10000, 60, 3.5)),
rating_min = rbeta(1000000, params[i, "min_rating_shape1"], params[i, "min_rating_shape2"]) * 5, going = rep(NA, 1000000),
waiting = rep(0, 1000000), sample_p = rep(0.5, 1000000),
satisfaction = rep(NA, 1000000), satis_random = rep(NA, 1000000),
satis_animals = rep(NA, 1000000), satis_price = rep(NA, 1000000),
satis_invest = rep(NA, 1000000), satis_wait = rep(NA, 1000000),
stringsAsFactors=FALSE)
# calculates encounter probabilities
# and maximum time that can be sustainably be spent with animals
# according to effect of tourism in the previous year
max_times[y] <- ifelse(y == 1, max_times[y], time_with_animals(maxx = max_times[y-1], effect = effects[y-1]))
encounter_probs[y] <- ifelse(y == 1, encounter_probs[y], p_encounter(p_e = encounter_probs[y-1], effect = effects[y-1]))
# create a list to store daily ratings
ratings <- vector("list", days)
# create a list to store daily behaviours
behaviours <- vector("list", days)
for(d in 1:days){
day_of_sim <- d + (y - 1) * days
### daily tourists here
tourists <- tourists_pop[sample(nrow(tourists_pop), n_tourists[day_of_sim], replace = FALSE, prob = tourists_pop$sample_p), ]
if(nrow(tourists) == 0){
next
}
bookings <- booking(tourists, as.data.frame(tour_ops))
# extract dataframes from list
tourists <- bookings[[1]]
tour_ops <- bookings[[2]]
# avoid tour operators' profits being negative
# if booking * price - costs < 0 then tour operator does not run tour
tour_ops$bookings <- ifelse(((tour_ops$bookings * tour_ops$price) - (1.5 * 90)) < 0, 0, tour_ops$bookings)
tourists[which(tourists$going %in% which(tour_ops$bookings == 0)), c("going", "waiting")] <- c(NA, +1)
# add daily bookings to year bookings
tour_ops$bookings_year <- tour_ops$bookings_year + tour_ops$bookings
# calculate time spent with animals
tour_ops <- tour_ops %>%
mutate(time_with = ifelse(bookings > 0, encounter_time(n_ops = length(tour_ops$id), p_e = encounter_probs[y]), 0))
# choose behavioural strategy (cooperate or defect).
# If the time calculated above is more than what is allowed
# the tour operator makes a choice according to costs and benefits
tour_ops <- behaviour_choice(tour_ops = tour_ops, payoff_CC = CC(tour_ops$price, tour_ops$capacity, max_times[y], 90, 15, length(tour_ops$id)),
payoff_CD = CD(tour_ops$price, tour_ops$capacity, max_times[y], tour_ops$time_with, 90, 15, length(tour_ops$id)),
payoff_DC = DC(tour_ops$price, tour_ops$capacity, tour_ops$time_with, max_times[y], 90, 15, detection_prob, fine, length(tour_ops$id)),
payoff_DD = DD(tour_ops$price, tour_ops$capacity, tour_ops$time_with, 90, 15, detection_prob, fine), max_times[y])
tour_ops$behaviour <- factor(tour_ops$behaviour, levels = c("defect", "cooperate", "no choice"))
# if tour operator cooperates time_with is the maximum allowed
# if defects time_with is the time calculated above
tour_ops <- tour_ops %>%
mutate(time_with = case_when(bookings == 0 ~ 0,
behaviour == "defect" | behaviour == "no choice" ~ time_with,
behaviour == "cooperate" ~ max_times[y] / days / 10))
# tour operators update time spent with animals in the year and profits
tour_ops <- tour_ops %>%
mutate(time_with_year = time_with_year + time_with,
profit = case_when(bookings == 0 ~ 0,
TRUE ~ (bookings * price) - (1.5 * 90)),
profit_year = profit_year + profit)
# tourists calculate satisfaction
tourists <- tourists %>%
group_by(going) %>%
mutate(satis_animals = ifelse(is.na(going), as.integer(NA),
satisfaction_animals(tour_ops[which(tour_ops$id == unique(going)), "time_with"], 90, 15)),
satis_price = ifelse(is.na(going), as.integer(NA),
satisfaction_price(tour_ops[which(tour_ops$id == unique(going)), "price"], max(tour_ops$price),
tour_ops[which(tour_ops$id == unique(going)), "rating"], max(tour_ops$rating), 15, 0.7)),
satis_wait = ifelse(is.na(going), as.integer(NA), satisfaction_waiting(waiting, -60,0.1)),
satis_invest = ifelse(is.na(going), as.integer(NA),
satisfaction_other_investment(tour_ops[which(tour_ops$id == unique(going)), "investment_ot"],
max(tour_ops$investment_ot), 10, 0.3)),
satis_random = ifelse(is.na(going), as.integer(NA),
rbinom(1, 1, p = (satis_animals * satis_price * satis_wait * satis_invest)))) %>%
ungroup() %>%
rowwise() %>%
mutate(satisfaction = ifelse(is.na(going), as.integer(NA),
ifelse(satis_random == 1,
sum(satis_animals, satis_price, satis_wait, satis_invest, na.rm = TRUE) +
(satis_animals * satis_price * satis_wait * satis_invest),
sum(satis_animals, satis_price, satis_wait, satis_invest, na.rm = TRUE) -
(satis_animals * satis_price * satis_wait * satis_invest))))
# tour operators update rating accroding to tourist satisfaction and clean dataset
tour_ops <- tour_ops %>%
group_by(id) %>%
mutate(rating = ifelse(bookings == 0, rating,
median(c(colMeans(tourists[which(tourists$going == id), "satisfaction"], na.rm = T), rating), na.rm = T)),
bookings = 0,
time_with = 0,
profit = 0) %>%
ungroup()
# store daily ratings
ratings[[d]] <- data.frame(id=tour_ops$id, rating=tour_ops$rating)
# store daily behaviours
behaviours[[d]] <- data.frame(id = tour_ops$id, behaviour = tour_ops$behaviour)
# merge tourists back into tourists_pop to keep the waiting counter and sampling probability
tourists_pop$waiting[match(tourists$id, tourists_pop$id)] <- tourists$waiting
tourists_pop$sample_p[match(tourists$id, tourists_pop$id)] <- tourists$sample_p
}
# calculate and store median rating and standard deviation
ratings_year[[y]] <- data.frame(id = tour_ops$id, year=rep(y,length(tour_ops$id)),
rating = aggregate(rating ~ id, data = do.call(rbind, ratings), median)[,2],
SD = aggregate(rating ~ id, data = do.call(rbind, ratings), sd)[,2])
# calculate and store daily behaviours of tour operators
behaviours_year[[y]] <- data.frame(id = tour_ops$id, year=rep(y,length(tour_ops$id)),
defect = aggregate(behaviour ~ id, data = do.call(rbind, behaviours), table)$behaviour[,"defect"],
cooperate = aggregate(behaviour ~ id, data = do.call(rbind, behaviours), table)$behaviour[,"cooperate"])
# fine defecting tour operators
# identify defectors
tour_ops_defect <- data.frame(id = behaviours_year[[y]][which(behaviours_year[[y]]$defect > 0), "id"])
# fine
if(dim(tour_ops_defect)[1] != 0){
tour_ops_defect$fine <- replicate(length(tour_ops_defect$id), {fines(detection_prob, fine)})
tour_ops[which(tour_ops$id %in% tour_ops_defect$id), "profit_year"] <- tour_ops[which(tour_ops$id %in% tour_ops_defect$id), "profit_year"] - tour_ops_defect$fine
}
# calculate tourism effect in the past year
wide_time <- 0.00025/(max_times[y]/100000)
wide_capacity <- 0.2 / (max_times[y]/100000)
effects[y] <- tourism_effect(slope_time = wide_time, slope_capacity = wide_capacity,
init_capacity = capacity_0, new_capacity = sum(tour_ops$capacity),
withanimals = sum(tour_ops$time_with_year), maxx = max_times[y])
# decide on infrastruucture investment
tour_ops <- tour_ops %>%
mutate(investment_ot = invest_services(rating = rating, max_rating = max(rating), profit = profit_year),
investment_infra = invest_infrastructure(profit = profit_year, max_profit = capacity * price * days, capacity = capacity, ticket = price)) %>%
mutate(investment_infra = ifelse(profit_year - (investment_infra + investment_ot) > 0, investment_infra, 0.001),
capacity = ifelse(investment_infra > 0.001, capacity + as.integer(profit / (capacity * price * 14)), capacity))
# store tour operators investments in a list
investments[[y]] <- data.frame(id = tour_ops$id, year=rep(y,length(tour_ops$id)),
infrastructure = tour_ops$investment_infra,
services = tour_ops$investment_ot)
# store tour operators annual profits in a list
profits[[y]] <- data.frame(id=tour_ops$id, year=rep(y,length(tour_ops$id)), money=tour_ops$profit_year)
# remeber to use
# profits <- do.call("rbind", profits)
# to transform into data.frame
# store tour operators bookings in a list
bookings_year[[y]] <- data.frame(id=tour_ops$id, year=rep(y,length(tour_ops$id)), bookings = tour_ops$bookings_year)
# store prices
prices[[y]] <- data.frame(id=tour_ops$id, year=rep(y,length(tour_ops$id)), ticket_price = tour_ops$price)
# store time spent with animals
withanimals[[y]] <- data.frame(id=tour_ops$id, year=rep(y,length(tour_ops$id)), time = tour_ops$time_with_year)
# operators who have had no profits for the past 3 years retire
if(y > 2) {bankrupt <- as.numeric(names(which(table(do.call("rbind", lapply(profits[(y-3):y], subset, money == 0))$id) == 3)))
tour_ops <- subset(tour_ops, !(id %in% bankrupt))}
if(nrow(tour_ops) == 0){break}
# tour operators update prices according to costs
tour_ops <- tour_ops %>%
mutate(price = case_when(price_change(ticket = tour_ops$price, demand = sum(tour_ops$bookings_year, na.rm = T), supply = sum(tour_ops$capacity) * days) > (1.5 * 90) / tour_ops$capacity ~ price_change(ticket = price, demand = sum(tour_ops$bookings_year, na.rm = T), supply = sum(tour_ops$capacity) *days),
TRUE ~ price))
# new operator start? Every 6 years new TOs can start
# probability of new operators wanting to start given by demand / supply ratio
# probability is used in a binomial draw
if(y %in% seq(6, years, 6)) {
p_to <- sum(tour_ops$bookings_year, na.rm = T) / (sum(tour_ops$capacity) * days)
new_tour_ops <- rbinom(1, 1, p = ifelse(p_to > 1, 1, p_to))}
else{new_tour_ops <- 0}
# if binomial draw is one and the time spent with the animals
# in the past year is smaller than the maximum time allowed
# by an amount that is equal to how much a tour operator is allowed,
# one new TO starts in the next year
if(new_tour_ops == 1 &&
mean(max_times[(y-6):y]) - mean(unlist(lapply(withanimals[(y-6):y], function(x) sum(x)))) >= mean(unlist(lapply(withanimals[(y-6):y], function(x) sum(x)))) / dim(tour_ops)[1]) {
tour_ops <- rbind(tour_ops, data.frame(id = max(tour_ops$id) + 1, price = rnorm(1, mean(tour_ops$price), 1), rating = 3,
capacity = as.integer(runif(1, 10, 30)), bookings = 0, bookings_year = 0, investment_infra = 0.001,
investment_ot = 0.001,time_with = 0,
time_with_year = 0, profit = 0, profit_year = 0, phenotype = sample(phenotypes, 1, replace=TRUE,
prob=c(p_trustful, p_optimist, p_pessimist, p_envious, p_undefined)),
behaviour = character(1), tours = 365, stringsAsFactors = FALSE))}
# set profits bookings and time with animals back to 0
tour_ops$profit_year <- 0
tour_ops$time_with_year <- 0
tour_ops$bookings_year <- 0
# keep track of simulation
print(y)
}
list(effect = effects, e_probs = encounter_probs, time_max = max_times,
income = profits, invest = investments, bookings = bookings_year,
ratings = ratings_year, tickets = prices, time_with = withanimals, behaviours = behaviours_year)
}
stopCluster(cl)
saveRDS(results, "../SimResults/Scenario2.rds")