Tutorials and Examples
Below are several examples that showcase the functionality of the package and how it can be used.
MC Simulation for the SD Model
A Monte Carlo simulation requires generating data for a model and estimating it. We now do this for the Search and Discovery Model.
First, we load the package. We also load the Distributions package because it specifies the Distribution type we're using here.
using StructuralSearchModels
using DistributionsNow, we create the model with defined parameter values.
# Define model
m = SD(
β = [-0.05, 3.0], # preference parameters
Ξ = 3.5, # discovery value
ρ = [-0.1], # parameter governing decrease in discovery value across positions
ξ = 2.5, # search value
dE = Normal(), # Distribution of εᵢⱼ (shocks revealed from search)
dV = Normal(), # Distribution of νᵢⱼ (shocks revealed on list)
dU0 = Uniform(), # distribution of outside option shock η
zdfun = "log" # functional form specification
)SD{Float64}
β: Array{Float64}((2,)) [-0.05, 3.0]
Ξ: Float64 3.5
ρ: Array{Float64}((1,)) [-0.1]
ξ: Float64 2.5
dE: Distributions.Normal{Float64}
dV: Distributions.Normal{Float64}
dU0: Distributions.Uniform{Float64}
zdfun: String "log"
information_structure: InformationStructureSpecification{Float64}
heterogeneity: HeterogeneitySpecification{Float64}
cs: Nothing nothing
cd: Nothing nothing
To generate sample data, we specify for how many consumers we want to generate data, and for how many sessions per consumer. We also specify a seed, which is not required but recommended. As a default, generate_data generates products with product characteristics drawn from a multivariate normal distribution and a last product characteristic which is a dummy for the outside option. In this case, this means that the last element of $\beta$ is the preference parameter for the outside option. Generating Products provides an example of how products with more characteristics can be generated.
n_consumers = 1000
n_sessions_per_consumer = 1
seed = 1
d = generate_data(m, n_consumers, n_sessions_per_consumer; seed)DataSD{Float64}
consumer_ids: Array{Int64}((1000,)) [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 991, 992, 993, 994, 995, 996, 997, 998, 999, 1000]
product_ids: Array{Vector{Int64}}((1000,))
product_characteristics: Array{Matrix{Float64}}((1000,))
positions: Array{Vector{Int64}}((1000,))
session_characteristics: Nothing nothing
consideration_sets: Array{Vector{Bool}}((1000,))
purchase_indices: Array{Int64}((1000,)) [1, 1, 1, 1, 8, 26, 1, 1, 1, 1 … 1, 1, 13, 1, 8, 1, 1, 1, 1, 1]
min_discover_indices: Array{Int64}((1000,)) [31, 2, 2, 2, 20, 26, 2, 2, 2, 2 … 3, 2, 20, 27, 27, 2, 2, 2, 2, 2]
search_paths: Array{Vector{Int64}}((1000,))
stop_indices: Array{Int64}((1000,)) [31, 2, 2, 2, 23, 26, 2, 2, 2, 2 … 4, 4, 20, 31, 31, 2, 2, 2, 2, 2]
The data now is all in the data object d, which is a DataSD type. By default, the data contains all fields, including where consumers stopped scrolling and the search order. These fields can be dropped for the simulation by setting them to nothing. This is useful if you want to simulate data that does not include this information.
d.stop_indices = nothingTo estimate the model, we need to specify an estimator. The package currently includes one estimator: SMLE. It can be set up with the following command. See the Estimation section for more details on available options. By default, the SMLE estimator uses 100 simulation draws.
e = SMLE(200) # Simulated MLE estimator with 200 simulation drawsSMLE
options_optimization: @NamedTuple{algorithm::Optim.LBFGS{Nothing, LineSearches.InitialStatic{Float64}, LineSearches.BackTracking{Float64, Int64}, Returns{Nothing}}, differentiation::ADTypes.AutoForwardDiff{nothing, Nothing}}
options_problem: Tuple{} ()
options_solver: @NamedTuple{iterations::Int64}
numerical_integration_method: DefaultNI
numerical_integration_method_heterogeneity: QMC
conditional_on_search: Bool false
parameter_rescaling: Nothing nothing
We now estimate the model using the estimate function. The function takes the model, the estimator, and the data as arguments. We also specify a seed for the RNG. The function returns the estimated model, the estimates, the likelihood at the estimates, the result of the optimization solver, and standard errors. Asymptotic standard errors in the SMLE are computed using the inverse Hessian.
seed = 1
m_hat, estimates, likelihood_at_estimates, result_solver,
std_errors = estimate(m, e, d; seed)
m_hatSD{Float64}
β: Array{Float64}((2,)) [-0.08035830614876124, 2.9269376685586357]
Ξ: Float64 3.41390257079509
ρ: Array{Float64}((1,)) [-0.08696878774523484]
ξ: Float64 2.42992507688633
dE: Distributions.Normal{Float64}
dV: Distributions.Normal{Float64}
dU0: Distributions.Uniform{Float64}
zdfun: String "log"
information_structure: InformationStructureSpecification{Float64}
heterogeneity: HeterogeneitySpecification{Float64}
cs: Nothing nothing
cd: Nothing nothing
To compare the estimates with the true parameter values, we can get the parameter values of the model using the vectorize_parameters function.
true_params = vectorize_parameters(m)
est_params = vectorize_parameters(m_hat)
hcat(est_params, true_params, std_errors)5×3 Matrix{Float64}:
-0.0803583 -0.05 0.0160144
2.92694 3.0 0.0472033
3.4139 3.5 0.0529486
-0.0869688 -0.1 0.00735926
2.42993 2.5 0.0491273If we want to compute standard errors after estimation, the calculate_standard_errors function can be used. When doing so, it is important to set the seed argument to the same value as used for the estimation.
std_errors = calculate_standard_errors(m_hat, e, d; seed)5-element Vector{Float64}:
0.016014439838011416
0.04720331518388762
0.0529485514873095
0.00735925965313398
0.049127325353165864MC Simulation for the WM Model
We can do the same for the WM model. The only difference is that we specify the model differently. The other functions are the same.
using StructuralSearchModels
using Distributions
# Define model
m = WM(
β = [-0.1, 3.0], # preference parameters
ρ = [-0.1], # parameter governing decrease in discovery value across positions
ξ = 1.5, # search value
dE = Normal(), # Distribution of εᵢⱼ (shocks revealed from search)
dV = Normal(), # Distribution of νᵢⱼ (shocks revealed on list)
dU0 = Uniform(), # distribution of outside option shock η
zsfun = "linear" # functional form specification
)
n_consumers = 1000
n_sessions_per_consumer = 1
seed = 1
d = generate_data(m, n_consumers, n_sessions_per_consumer; seed)
e = SMLE(200) # Simulated MLE estimator with 200 simulation draws
seed = 1
m_hat, estimates, likelihood_at_estimates,
result_solver, std_errors = estimate(m, e, d; seed) ;
m_hatWM{Float64}
β: Array{Float64}((2,)) [-0.06351499678944007, 3.029058517962158]
ξ: Float64 1.495311841795068
ρ: Array{Float64}((1,)) [-0.08756232802224653]
dE: Distributions.Normal{Float64}
dV: Distributions.Normal{Float64}
dU0: Distributions.Uniform{Float64}
zsfun: String "linear"
information_structure: InformationStructureSpecification{Float64}
heterogeneity: HeterogeneitySpecification{Float64}
cs: Nothing nothing
Adding Search Over Vertical Attributes
To specify that vertical attributes are revealed when searching alternative, we additionally specify the search_value_specification. Below is an example where the first attribute in the data.product_characteristics is revealed on the list. This attribute enters utility through $x_j^l\beta$ and also informs consumers beliefs about what will be revealed when searching through $x_j^l\gamma$. The second characteristic is revealed on search and enters only $x_j^d\kappa$. The following code sets up this model. Note that the other elements in the three parameter vectors are set to zero. Moreover, note that the last element in $\beta$ is for the outside option.
m = SD(
β = [-0.05, 0.0, 0.8],
Ξ = 1.0,
ρ = [-0.1],
ξ = 1.0,
dE = Normal(),
dV = Normal(),
dU0 = Normal(),
zdfun = "linear",
information_structure = InformationStructureSpecification(
γ = [0.02, 0.0, 0.0],
κ = [0.0, 0.1, 0.0],
indices_characteristics_β_union = 1:1,
indices_characteristics_γ_union= 1:0,
indices_characteristics_κ_union = 2:2,
)
)SD{Float64}
β: Array{Float64}((3,)) [-0.05, 0.0, 0.8]
Ξ: Float64 1.0
ρ: Array{Float64}((1,)) [-0.1]
ξ: Float64 1.0
dE: Distributions.Normal{Float64}
dV: Distributions.Normal{Float64}
dU0: Distributions.Normal{Float64}
zdfun: String "linear"
information_structure: InformationStructureSpecification{Float64}
heterogeneity: HeterogeneitySpecification{Float64}
cs: Nothing nothing
cd: Nothing nothing
Computing Consumer Welfare and Revenues
To compute consumer welfare and revenues for an estimated model m_hat and data d, we can use the calculate_welfare and calculate_revenues functions.
First, we need to compute the costs using the calculate_costs! function. This function computes the search and discovery costs from the estimates for $\xi$ and $\Xi$ and adds them to the model in-place. Cost computation is done using n_draws_cd draws of the attribute distribution.
Once the costs are computed, we can compute welfare, which gives a Dict with welfare values. The dictionary reports the total welfare, the paid discovery costs, and the effective value of the chosen alternative. This is done as an average across all sessions, only those that had at least one search, and only those resulted in a purchase.
Going back to the example for the SD model, we can compute the costs and welfare as follows:
# Compute costs
n_draws_cd = 100_000
calculate_costs!(m_hat, d, n_draws_cd; seed)
# Compute welfare
n_draws_welfare = 1000
calculate_welfare(m_hat, d, n_draws_welfare; seed)Dict{Any, Any} with 3 entries:
:average => Dict(:discovery_costs=>0.033262, :utility=>3.4782…
:conditional_on_purchase => Dict(:discovery_costs=>0.0692072, :utility=>3.604…
:conditional_on_search => Dict(:discovery_costs=>0.0747511, :utility=>3.307…To compute revenues, we can use the calculate_revenues function. It requires specifying which product characteristic is the price attribute and the number of draws to use for the demand computation of each product. Note, this computation is done using the simulation procedure described in Greminger (2025).
index_price_attribute = 1
n_draws_demand = 50
total_revenues, revenues_by_session = calculate_revenues(m_hat, d, index_price_attribute, n_draws_demand; seed)Dict{Symbol, Any} with 4 entries:
:demand_individual => [0.1, 0.12, 0.04, 0.06, 0.08, 0.16, 0.1, 0.12, 0.08, …
:revenues_individual => [0.0802777, -0.0360898, -0.0531847, -0.0482984, -0.02…
:demand => 97.38
:revenues => -15.8848Computing Fit Measures
To evaluate how well an estimated model fits the data, we can use the calculate_fit_measures function. This function computes the fit measures for the estimated model and the data. It requires specifying the number of draws to use for the simulation. It returns a dictionary with the fit measures.
For the previously estimated SD model, we can compute the fit measures as follows:
n_draws = 1000
fit_measures = calculate_fit_measures(m_hat, d, n_draws; seed)Dict{Symbol, Any} with 7 entries:
:bounds => [([0.0, 0.148, 0.082, 0.077, 0.072, 0.067, 0.064,…
:click_stats_sim => Dict{Symbol, Any}(:mean_position_clicked=>11.4798…
:stop_probabilities_sim => [0.0, 0.577432, 0.038329, 0.028247, 0.02273, 0.01…
:click_stats_data => Dict{Symbol, Any}(:mean_position_clicked=>11.3608…
:purchase_stats_sim => Dict{Symbol, Any}(:purchase_probability=>0.098759…
:stop_probabilities_data => [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0…
:purchase_stats_data => Dict{Symbol, Any}(:purchase_probability=>0.089, :…These fit measures then can be accessed in the usual way. For example, if we want to compare the position-specific click probabilities, we can do this as follows.
prob_sim = fit_measures[:click_stats_sim][:click_probability_per_position]
prob_data = fit_measures[:click_stats_data][:click_probability_per_position]
hcat(prob_sim, prob_data) # compare simulated and data click probabilities31×2 Matrix{Float64}:
0.0 0.0
0.166502 0.161
0.09736 0.101
0.091599 0.094
0.086973 0.083
0.081229 0.086
0.077343 0.082
0.074299 0.089
0.071032 0.067
0.068552 0.061
⋮
0.045369 0.041
0.043816 0.049
0.042477 0.047
0.041773 0.034
0.0405 0.037
0.039401 0.045
0.038005 0.029
0.037487 0.039
0.036225 0.036Or if we want to compare the average number of clicks per session, we can do this as follows.
n_clicks_sim = fit_measures[:click_stats_sim][:no_clicks_per_session]
n_clicks_data = fit_measures[:click_stats_data][:no_clicks_per_session]
hcat(n_clicks_sim, n_clicks_data) # compare simulated and data clicks per session1×2 Matrix{Float64}:
1.84396 1.81Estimating Shock Variances
To estimate the variances of the unobserved shocks, we can pass the estimation_shock_variances keyword argument to the estimate function. The following would estimate the variance of the unobserved shock dE. Note that trying to run this code will not produce valid results. This is because the variance is not identified in this example, which is not checked by the code.
See the Estimation section for more details on the options.
estimation_shock_variances = [:σ_dE]
m_hat, estimates, likelihood_at_estimates,
result_solver, std_errors = estimate(m, e, d; seed, estimation_shock_variances)Specifying Solver Options
Solver options can be passed to the estimator. By default, the SMLE estimator uses the LBFGS optimizer, with gradients and Hessians computed using automatic differentiation. Specifically, calling e=SMLE(100) is equivalent to calling the following:
using Optimization # need to load to specify AutoDiff
using OptimizationOptimJL # provides Optim algorithms
using LineSearches
e = SMLE(100;
options_optimization = (
algorithm = LBFGS(; linesearch = LineSearches.BackTracking(order = 2)),
differentiation = Optimization.AutoForwardDiff()
)
)It uses the LBFGS algorithm with the BackTracking line search algorithm, which I found to be a bit more stable for these models.
While this is the default, we can use any solver/algorithm accessible through the Optimization package, provided we load the necessary packages first. For example, we can use the NelderMead algorithm as follows:
using Optimization # need to load to specify AutoDiff
using OptimizationOptimJL # provides Optim algorithms (LBFGS, NelderMead etc.)
e = SMLE(100;
options_optimization = (
algorithm = NelderMead(),
differentiation = Optimization.AutoForwardDiff()
)
)
Note that the NelderMead algorithm does not use gradients, but the differentiation still needs to be specified as it has no default.
We can also pass additional options to the optimizer that are used by the Optimization.jl package as described in the estimation section. For example, we can print information on the solvers' iterations as follows:
e = SMLE(100;
options_solver = (show_trace = true, show_every = 10)
)This again uses the default LBFGS optimizer. These options may differ across optimizers. To see what's all possible, check the estimation section and the documentation of the Optimization.jl package.
Parameter Rescaling
When model parameters differ greatly in how they affect the likelihood, the optimizer can take longer to converge. The parameter_rescaling option lets you provide a vector of scaling factors so that the optimizer works in a normalized space. The package offers a convenient way to construct a suitable scaling vector via build_inverse_hessian_scaler. This function approximates the curvature of the likelihood at a given point (which usually would be the starting values)and returns 1 ./ sqrt.(diag_H), which is a scaling vector ready to pass to the parameter_rescaling option.
# Compute scaling vector from diagonal Hessian at starting values
x0 = vectorize_parameters(m)
e = SMLE(100)
scale = build_inverse_hessian_scaler(m, e, d, x0)
# Update estimator with scaling vector
e.parameter_rescaling = scale
# Estimate with the updated estimator
m_hat, estimates, likelihood_at_estimates, result_solver, std_errors = estimate(m, e, d)The scaling is applied transparently: the optimizer works with rescaled parameters θ_rescaled = θ ./ scale, while all outputs (estimates, standard errors) are returned on the original scale.
Conditioning on Search
Sometimes, it is necessary to condition on search, for example, when estimating the model on data that only includes sessions with at least one click/search. This can be accommodated by setting the conditional_on_search option to true. The keyword applies to both the generate_data function and the SMLE estimator. Importantly, the keyword also needs to be set when calculating fit measures based on the estimated model. If the model is estimated conditional on search, while the fit measures are computed without it, the estimated model will likely not fit the data well.
Below is an example for an MC simulation for the WM model that conditions on search.
using StructuralSearchModels
using Distributions
# Define model
m = WM(
β = [-0.1, 3.0], # preference parameters
ρ = [-0.1], # parameter governing decrease in discovery value across positions
ξ = 1.5, # search value
dE = Normal(), # Distribution of εᵢⱼ (shocks revealed from search)
dV = Normal(), # Distribution of νᵢⱼ (shocks revealed on list)
dU0 = Uniform(), # distribution of outside option shock η
zsfun = "linear" # functional form specification
)
n_consumers = 1000
n_sessions_per_consumer = 1
seed = 1
d = generate_data(m, n_consumers, n_sessions_per_consumer; seed, conditional_on_search = true)
e = SMLE(200; conditional_on_search = true) # Simulated MLE estimator with 200 simulation draws
m_hat, estimates, likelihood_at_estimates, result_solver, std_errors = estimate(m, e, d; seed) ;
n_sim = 1000
fit_measures = calculate_fit_measures(m_hat, d, n_sim; seed, conditional_on_search = true)
m_hatWM{Float64}
β: Array{Float64}((2,)) [-0.0915665137081704, 2.614781940107058]
ξ: Float64 1.4196795462369074
ρ: Array{Float64}((1,)) [-0.09842784675927765]
dE: Distributions.Normal{Float64}
dV: Distributions.Normal{Float64}
dU0: Distributions.Uniform{Float64}
zsfun: String "linear"
information_structure: InformationStructureSpecification{Float64}
heterogeneity: HeterogeneitySpecification{Float64}
cs: Nothing nothing
Generating Products
The generate_data function generates products using the generate_products function with default options. By default, it generates many different products, so that the products for each session are different. For each session, 30 products are sampled from the total number of products generated. For each product, product characteristics are drawn from a multivariate normal distribution, and the last product characteristic is assumed to be a dummy for the outside option.
The following generates products with a single product characteristic using the default options.
using Distributions, StructuralSearchModels
n_sessions = 2
product_ids, product_characteristics = generate_products(n_sessions, Normal())
product_characteristics[1] # show the product characteristics for the first session31×2 Matrix{Float64}:
0.0 1.0
0.741352 0.0
-0.857603 0.0
-1.25114 0.0
0.20667 0.0
0.000541248 0.0
-1.58682 0.0
0.488452 0.0
1.35961 0.0
-1.72073 0.0
⋮
-0.749141 0.0
-0.513452 0.0
-0.49227 0.0
-0.6667 0.0
-0.280047 0.0
-1.73937 0.0
2.06415 0.0
-2.54269 0.0
-0.0935741 0.0We can update the defaults by specifying keyword arguments. The following generates 100 products, of which, 10 are sampled for each session. There are now five product characteristics, and the last one is a dummy for the outside option.
using LinearAlgebra # provides convenient identify matrix constructor I(n)
n_sessions = 2
product_ids, product_characteristics = generate_products(n_sessions, MvNormal(I(5));
n_products = 100, n_products_per_session = 10)
product_characteristics[1] # show the product characteristics for the first session11×6 Matrix{Float64}:
0.0 0.0 0.0 0.0 0.0 1.0
0.614065 1.57919 -1.51982 -0.704462 -1.50624 0.0
0.14925 0.0179227 0.041416 1.49609 -0.803541 0.0
-1.06478 0.254082 0.118937 0.256643 1.23644 0.0
0.20066 1.8167 0.410969 -0.601564 0.252858 0.0
-1.42424 -0.452639 -1.39084 0.0210551 0.782351 0.0
1.10252 -0.409713 -1.02907 0.608061 -0.590686 0.0
-1.35606 0.664026 0.558733 -0.640729 0.801281 0.0
-1.72077 -0.474351 0.177091 0.28531 0.611775 0.0
0.594312 -1.42816 1.0473 1.51966 -0.871822 0.0
-0.183216 -0.414896 1.08171 -0.315755 -0.256567 0.0We can also suppress the outside option indicator by adding outside_option = false.
product_ids, product_characteristics = generate_products(n_sessions, MvNormal(I(5));
n_products = 100, n_products_per_session = 10, outside_option = false)
product_characteristics[1] # show the product characteristics for the first session10×5 Matrix{Float64}:
0.00639063 1.07431 -0.380136 -0.195797 1.34874
-0.455624 -1.09156 0.532969 2.14716 -0.401033
-0.821375 0.0584874 1.45876 0.0906896 -0.170223
-1.42273 0.800884 -0.612195 1.25007 1.87401
0.840979 -0.944985 -1.42842 1.47113 0.931302
0.676918 0.827983 0.842577 -0.422245 -0.938152
-0.178433 -0.154088 0.29332 -0.0629499 0.125209
0.746974 1.5712 0.519627 -0.636784 -1.14582
0.224812 -1.46657 -0.660883 0.461805 -0.52447
-1.30006 -1.00155 -0.494485 0.488309 -0.306919Preparing Data
To estimate the model on real data, the data needs to be prepared. The following example shows how to do this when the data is in a DataFrame. The code loops over the sessions (which is fast in Julia) and creates the necessary structure for the DataSD type.
using DataFrames, StructuralSearchModels
# Generate example data
df = DataFrame(
session_id = [1, 1, 1, 2, 2, 2],
position = [0, 1, 2, 0, 1, 2],
product_id = [1, 2, 3, 1, 2, 3],
price = [1.0, 2.0, 3.0, 1.0, 2.0, 3.0],
quality = [0.5, 0.6, 0.7, 0.5, 0.6, 0.7],
discovered = [1, 1, 0, 1, 1, 1],
clicked = [1, 0, 0, 0, 1, 0],
purchased = [1, 0, 0, 0, 0, 0],
)
# Extract consumer ids. Have only one session per id.
n_sessions = length(unique(df.session_id))
consumer_ids = collect(1:n_sessions)
# Create empty arrays to store the data
product_ids = []
product_characteristics = []
positions = []
indices_list_characteristics = []
session_characteristics = []
consideration_sets = []
purchase_indices = []
stop_indices = []
# loop over sessions and add to arrays. For each session, adds outside option
for subdf in groupby(df, :session_id)
# Product ids
product_ids_i = vcat(0, subdf.product_id)
push!(product_ids, product_ids_i)
# Product characteristics
pc_i = Matrix(subdf[:, [:price, :quality]]) # converts to matrix
pc_i = hcat(pc_i, zeros(size(pc_i, 1))) # add characteristic for outside option
pc_i = vcat(zeros(1, size(pc_i, 2)), pc_i) # add outside option as product
pc_i[1, end] = 1 # set outside option dummy to 1 for outside option
push!(product_characteristics, pc_i)
# Positions
positions_i = vcat(0, subdf.position)
push!(positions, positions_i)
# Consideration sets
cs_i = vcat(0, subdf.clicked)
push!(consideration_sets, cs_i)
# stop indices
si = findlast(x -> x == 1, subdf.discovered) # finds last 1 in discovered
push!(stop_indices, isnothing(si) ? length(subdf.discovered) + 1 : si + 1) # + 1
# purchase indices
pi = findfirst(x -> x == 1, subdf.purchased) # returns nothing if not found
push!(purchase_indices, isnothing(pi) ? 1 : pi + 1) # +1 for outside option
end
# Put together into data type
d = DataSD(; consumer_ids,
product_ids,
product_characteristics,
positions,
consideration_sets,
purchase_indices,
stop_indices)DataSD{Float64}
consumer_ids: Array{Int64}((2,)) [1, 2]
product_ids: Array{Vector{Int64}}((2,))
product_characteristics: Array{Matrix{Float64}}((2,))
positions: Array{Vector{Int64}}((2,))
session_characteristics: Nothing nothing
consideration_sets: Array{Vector{Bool}}((2,))
purchase_indices: Array{Int64}((2,)) [2, 1]
min_discover_indices: Nothing nothing
search_paths: Nothing nothing
stop_indices: Array{Int64}((2,)) [3, 4]
Changing the Number of Products in Each Position
The number of products in each position can be set by adjusting the positions in the DataSD type. The package provides a function to do this. The following updates the positions in the data object d to have 3 products in the first position (those revealed without cost), and then 2 products per position after that.
update_positions!(d, 3, 2)To generate data with a different number of positions, we can also pass the keyword arguments n_A0 and n_d to the generate_data function. The following generates data with 3 products in the first position (those in the initial awareness set), and then 2 products per position after that.
generate_data(model, n_consumers, n_sessions_per_consumer;
n_A0 = 3, n_d = 2)