VERSION FINAL PERIODT PERIODT
library(deSolve)
library(ggplot2)
library(doParallel)
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
# Define the time vector for simulation
times <- seq(0, 40, by = 1)
# Set the initial conditions of the state variable
initial_conditions <- c(
children = 4000000,
sexually_active_adults = 12000000,
supply_of_prostitutes = 28254
)
# Define the model parameters
parameters <- c(
age_of_consent = 18, # years
sexual_lifetime = 45, # years
birth_rate = 17 / 1000, # 17 children per 1000 sexually active adult
normal_fraction_of_young_adults_willing_to_consider_a_job_in_prostitutuion = 1 / 1000, # 1 in 1000 youngsters
average_lifetime_in_prostitution = 10, # years
delivery_time = 1, # year
naked_costs_per_sex_service = 20, # pounds
clients_per_prostitute = 50, # clients
normal_percentage_of_johns = 12 / 100, # 12%
average_visit_frequency = 24 # times per year
)
# Define the function to calculate average_price_of_sex_services
calc_average_price_of_sex_services <- function(t, naked_costs_per_sex_service, supply_demand_effect_on_the_price_of_sex_services) {
if (t == 0) {
return(2 * naked_costs_per_sex_service)
} else {
return(2 * naked_costs_per_sex_service * supply_demand_effect_on_the_price_of_sex_services)
}
}
# Define the model function
modelo <- function(t, state, parameters) {
with(as.list(c(state, parameters)), {
# Flow variables
prostitution_outflow = supply_of_prostitutes / average_lifetime_in_prostitution
# Auxiliary or endogenous variables
societal_acceptance_of_prostitution = (1 - supply_of_prostitutes / 100000) / 3 # avg delay time of 3 years
effect_of_societal_acceptance_on_the_percentage_of_johns = (0.5 + societal_acceptance_of_prostitution / 1.5)
price_effect_on_the_demand_for_prostitution = approx(
c(0, 50, 100, 200, 400, 800, 1200), # x = avg_price_of_sex_services
c(4, 1.5, 1, 0.5, 0.2, 0.1, 0.1), # y = price_effect_on_the_demand_for_prostitution
xout = 2 * naked_costs_per_sex_service)$y # placeholder for xout
demand_for_prostitutes = ((sexually_active_adults / 2) * normal_percentage_of_johns * effect_of_societal_acceptance_on_the_percentage_of_johns * price_effect_on_the_demand_for_prostitution) / clients_per_prostitute
supply_demand_ratio = supply_of_prostitutes / demand_for_prostitutes
supply_demand_effect_on_the_price_of_sex_services = approx(
c(0, 0.05, 0.25, 1, 2, 3), # x = supply_demand_ratio
c(500 / 100, 300 / 100, 200 / 100, 100 / 100, 60 / 100, 50 / 100), # y = supply_demand_effect_on_the_price_of_sex_services
xout = supply_demand_ratio)$y # $y para que solo arroje y
average_price_of_sex_services = calc_average_price_of_sex_services(t, naked_costs_per_sex_service, supply_demand_effect_on_the_price_of_sex_services)
price_effect_on_the_demand_for_prostitution = approx(
c(0, 50, 100, 200, 400, 800, 1200), # x = avg_price_of_sex_services
c(4, 1.5, 1, 0.5, 0.2, 0.1, 0.1), # y = price_effect_on_the_demand_for_prostitution
xout = average_price_of_sex_services)$y # $y para que solo arroje y
demand_for_prostitutes = ((sexually_active_adults / 2) * normal_percentage_of_johns * effect_of_societal_acceptance_on_the_percentage_of_johns * price_effect_on_the_demand_for_prostitution) / clients_per_prostitute
# Auxiliary or endogenous variables
profitability_of_prostitution = (average_price_of_sex_services - naked_costs_per_sex_service) / naked_costs_per_sex_service
expected_gap_between_demand_and_supply = demand_for_prostitutes + prostitution_outflow - supply_of_prostitutes
# Flow variables
births = sexually_active_adults * birth_rate
from_minor_to_adult = children / age_of_consent
loss_of_sexual_activity = sexually_active_adults / sexual_lifetime
new_dutch_prostitutes = (from_minor_to_adult / 2) * normal_fraction_of_young_adults_willing_to_consider_a_job_in_prostitutuion * profitability_of_prostitution * societal_acceptance_of_prostitution
new_foreign_prostitutes = (expected_gap_between_demand_and_supply * profitability_of_prostitution) / delivery_time
total_annual_prostitution_revenues = supply_of_prostitutes * average_visit_frequency * clients_per_prostitute * average_price_of_sex_services
# State variables
dchildren <- births - from_minor_to_adult
dsexually_active_adults <- from_minor_to_adult - loss_of_sexual_activity
dsupply_of_prostitutes <- new_dutch_prostitutes + new_foreign_prostitutes - prostitution_outflow
# Return the results of the state variables
return(list(
c(dchildren, dsexually_active_adults, dsupply_of_prostitutes),
births = births,
from_minor_to_adult = from_minor_to_adult,
loss_of_sexual_activity = loss_of_sexual_activity,
new_foreign_prostitutes = new_foreign_prostitutes,
new_dutch_prostitutes = new_dutch_prostitutes,
prostitution_outflow = prostitution_outflow,
demand_for_prostitutes = demand_for_prostitutes,
societal_acceptance_of_prostitution = societal_acceptance_of_prostitution,
average_price_of_sex_services = average_price_of_sex_services # Add this to output
))
})
}
# Integration method
intg.method <- "rk4"
# Perform the simulation using the 'ode' function from the deSolve package
out <- ode(
y = initial_conditions, # Initial conditions
times = times, # Simulation time
func = modelo, # Model function
parms = parameters, # Parameters
method = intg.method # Integration method
)
# Convert the output to a data frame
out <- as.data.frame(out)
# Plot the simulation results
ggplot(out, aes(x = time, y = supply_of_prostitutes)) + geom_line() + ggtitle("Supply of Prostitutes Over Time")
## EXPLICACIÓN
ggplot(out, aes(x = time, y = children)) + geom_line() + ggtitle("Number of Children Over Time")
## EXPLICACIÓN
ggplot(out, aes(x = time, y = sexually_active_adults)) + geom_line() + ggtitle("Number of Sexually Active Adults Over Time")
## EXPLICACIÓN
ggplot(out, aes(x = time, y = prostitution_outflow)) + geom_line() + ggtitle("Prostitution Outflow Over Time")
## EXPLICACIÓN
ggplot(out, aes(x = time, y = average_price_of_sex_services)) + geom_line() + ggtitle("Average Price of Sex Services Over Time")
## EXPLICACIÓN
# Reduce the range and steps of the uncertainties
x1 <- seq(10/100, 20/100, 5/100) # clients per prostitute
x2 <- seq(10/100, 20/100, 5/100) # naked costs per sex service
x3 <- seq(20/100, 30/100, 5/100) # visit frequency
x4 <- seq(20/100, 30/100, 5/100) # birth rate
x5 <- seq(20/100, 30/100, 5/100) # normal fraction of young dutch adults willing to consider a job in prostitution
Xs <- expand.grid(list(clients_of_prostitutes = x1, naked_costs_per_sex_service = x2, visit_frequency = x3, birth_rate = x4, normal_fraction_of_young_dutch_adults_willing_to_consider_a_job_in_prostitution = x5))
Xs$Run.ID <- 1:nrow(Xs)
# Set up parallel backend
numCores <- parallel::detectCores() - 1
cl <- parallel::makeCluster(numCores)
registerDoParallel(cl)
# Perform the simulations in parallel
out_all <- foreach(i = 1:nrow(Xs), .combine = rbind, .packages = c("deSolve")) %dopar% {
parameters.Xs <- parameters
parameters.Xs["clients_of_prostitutes"] <- Xs$clients_of_prostitutes[i]
parameters.Xs["naked_costs_per_sex_service"] <- Xs$naked_costs_per_sex_service[i]
parameters.Xs["visit_frequency"] <- Xs$visit_frequency[i]
parameters.Xs["birth_rate"] <- Xs$birth_rate[i]
parameters.Xs["normal_fraction_of_young_dutch_adults_willing_to_consider_a_job_in_prostitution"] <- Xs$normal_fraction_of_young_dutch_adults_willing_to_consider_a_job_in_prostitution[i]
out <- ode(y = initial_conditions, times = times, func = modelo, parms = parameters.Xs, method = intg.method)
out <- data.frame(out)
out$Run.ID <- Xs$Run.ID[i]
out$clients_of_prostitutes <- Xs$clients_of_prostitutes[i]
out$naked_costs_per_sex_service <- Xs$naked_costs_per_sex_service[i]
out$visit_frequency <- Xs$visit_frequency[i]
out$birth_rate <- Xs$birth_rate[i]
out$normal_fraction_of_young_dutch_adults_willing_to_consider_a_job_in_prostitution <- Xs$normal_fraction_of_young_dutch_adults_willing_to_consider_a_job_in_prostitution[i]
return(out)
}
# Stop the parallel backend
stopCluster(cl)
# Convert to data frame and ensure it contains the expected columns
out_all <- as.data.frame(out_all)
# Check for NA values
na_counts <- colSums(is.na(out_all))
print(na_counts)
## time
## 0
## children
## 0
## sexually_active_adults
## 0
## supply_of_prostitutes
## 216
## births
## 0
## from_minor_to_adult
## 0
## loss_of_sexual_activity
## 0
## new_foreign_prostitutes
## 270
## new_dutch_prostitutes
## 270
## prostitution_outflow
## 216
## demand_for_prostitutes
## 270
## societal_acceptance_of_prostitution
## 216
## average_price_of_sex_services
## 270
## Run.ID
## 0
## clients_of_prostitutes
## 0
## naked_costs_per_sex_service
## 0
## visit_frequency
## 0
## birth_rate
## 0
## normal_fraction_of_young_dutch_adults_willing_to_consider_a_job_in_prostitution
## 0
# Remove rows with NA values if necessary
out_all <- na.omit(out_all)
# Verify that the columns are present
print(colnames(out_all))
## [1] "time"
## [2] "children"
## [3] "sexually_active_adults"
## [4] "supply_of_prostitutes"
## [5] "births"
## [6] "from_minor_to_adult"
## [7] "loss_of_sexual_activity"
## [8] "new_foreign_prostitutes"
## [9] "new_dutch_prostitutes"
## [10] "prostitution_outflow"
## [11] "demand_for_prostitutes"
## [12] "societal_acceptance_of_prostitution"
## [13] "average_price_of_sex_services"
## [14] "Run.ID"
## [15] "clients_of_prostitutes"
## [16] "naked_costs_per_sex_service"
## [17] "visit_frequency"
## [18] "birth_rate"
## [19] "normal_fraction_of_young_dutch_adults_willing_to_consider_a_job_in_prostitution"
# Plot results
ggplot(out_all, aes(x = time, y = supply_of_prostitutes, group = Run.ID, colour = birth_rate)) +
geom_line() +
scale_color_gradient(low = "blue", high = "orange")
## Se pueden mostrar diferentes escenarios donde se observa una fuerte correlación entre las variables "birth_rate" y "supply_of_prostitutes". Dado que las variables que definen "supply_of_prostitutes" están determinadas por "birth_rate", se puede esperar que muestren una correlación positiva en cada escenario diferente.
ggplot(out_all, aes(x = time, y = supply_of_prostitutes, group = Run.ID, colour = clients_of_prostitutes)) +
geom_line() +
scale_color_gradient(low = "blue", high = "orange")
## Se puede observar en esta gráfica cómo influye la cantidad de clientes de las prostitutas en la oferta de las mismas. Esto tiene sentido al tener en cuenta que la variable “clients_of_prostitutes”, en los diferentes escenarios de incertidumbre, sigue la tendencia de correlación positiva al constituir las variables que eventualmente forman “supply_of_prostitutes”.
ggplot(out_all, aes(x = time, y = sexually_active_adults, group = Run.ID, colour = visit_frequency)) +
geom_line() +
scale_color_gradient(low = "blue", high = "orange")
## En el gráfico puede observarse como la variable “visit_frequency” y “sexually_active_adults” están directamente relacionadas. Sexually active adults, especially decir la cantidad de hombres sexualmente activos afecta directamente la cantidad de visitas que reciben las prostitutas.
ggplot(out_all, aes(x = time, y = sexually_active_adults, group = Run.ID, colour = naked_costs_per_sex_service)) +
geom_line() +
scale_color_gradient(low = "blue", high = "orange")
## En el gráfico se puede observar la variable “naked_costs_per_sex_services” con una correlación positiva con “sexually_active_adults”. Esto es congruente, puesto que sexually_active_adults está directamente ligado a la demanda de prostitutas, y por ley de oferta y demanda, tendería a aumentar el precio de un servicio para cumplir con el abasto.
ggplot(out_all, aes(x = time, y = children, group = Run.ID, colour = new_foreign_prostitutes)) +
geom_line() +
scale_color_gradient(low = "blue", high = "orange")
## En el gráfico se muestra una tendencia entre la variable “children” y la variable “new_foreign_prostitutes”. Dicha tendencia se debe a que si hay mas niños, habrán mas adultos sexualmente activos, lo que hará que se aumente la demanda y por ende la oferta de prostitutas.
ggplot(out_all, aes(x = time, y = children, group = Run.ID, colour = new_dutch_prostitutes)) +
geom_line() +
scale_color_gradient(low = "blue", high = "orange")
## En el gráfico se muestra una tendencia entre la variable “children” y la variable “new_dutch_prostitutes”. Dicha tendencia se debe a que si hay mas niños, habrán mas adultos sexualmente activos, lo que hará que se aumente la demanda y por ende la oferta de prostitutas.