Function reference

ParameterShock

A public struct encoding a perturbation to a GTAP model parameter.

Supports both commodity/endowment-specific shocks and region-specific shocks.

Fields

• parameter::String: The name of the parameter in model.data to shock
• indices::Union{Nothing, String, Vector{String}}: Row indices to apply the shock to (typically commodities or endowments).
  • nothing: Apply to all rows
  • String: Apply to a single row
  • Vector{String}: Apply to multiple rows
• scale::Union{Float64, Vector{Float64}, NamedArray}: The scaling factor $> 0.0$ to apply. Pass a vector of the same length as indices to scale each index separately. Alternatively, pass a NamedArray with rows representing indices (sectors, endowments) and columns representing regions to shock specific regions.

Constructors

# Basic constructor: commodity-specific or uniform regional scaling
ParameterShock(parameter::String, indices, scale::Float64)

# Region-specific: apply scaling to a commodity or endowment in specific regions
ParameterShock(parameter::String, indices, scale::NamedArray)

Examples

# Shock skilled and unskilled labor uniformly by 50%
ParameterShock("qe", ["skilled labour", "unskilled labour"], 0.5)

# Shock all rows of a parameter uniformly
ParameterShock("qe", nothing, 0.5)

# Apply shocks to specific commodities or endowments in specific regions
mag = reshape([0.1, 0.2], 2, 1)
mag = NamedArray(mag, (["skilled labour", "unskilled labour"], ["chn"]))
ParameterShock("qpa", ["skilled labour", "unskilled labour"], mag)

_apply_shocks!(data, base_data, shocks::Vector{ParameterShock})

Internal helper that applies a vector of parameter shocks to model data.

Mutates data in place, scaling values from base_data according to each shock's scale and region-specific factors.

Arguments

• data: The mutable model data dictionary to be modified
• base_data: The calibrated baseline data (used as reference for scaling)
• shocks::Vector{ParameterShock}: Vector of shocks to apply (may include region-specific shocks)

Details

Applies shocks in the following order:

1. Commodity/endowment-specific scaling (via indices and scale)
2. Region-specific scaling (via regions and region_scale, if provided)

For region-specific shocks, scaling is applied to the last dimension(s) of the parameter array (assumed to be regions).

compare_gtaps(initial::GTAP.model_container_struct, current::GTAP.model_container_struct)

Compare baseline and shocked GTAP models to extract economic outcomes.

Computes changes in economic indicators between two model states by calculating:

1. Absolute income/GDP by region in the current (shocked) model
2. Equivalent variation (welfare change relative to baseline) by region
3. Proportional change in GDP by region between baseline and current equilibrium

Arguments

• initial::GTAP.model_container_struct: Baseline model (typically the original or baseline equilibrium)
• current::GTAP.model_container_struct: Shocked or alternative model (typically the result of shock_gtap)

Returns

A named tuple (y, ev, delta_gdp) containing:

• y::NamedArray: Absolute income/GDP by region in the current model
• ev::NamedArray: Equivalent variation (change in expenditure function) by region, measuring welfare change relative to the initial equilibrium
• delta_gdp::NamedArray: Proportional change in GDP by region, calculated as (GDP_current / GDP_initial) - 1.0

Details

• delta_gdp is a proportional measure and can have negative values.
• All results are indexed by region names from the model's region set.

Example

# Create a shocked model
labour_shock = ParameterShock("qe", ["skilled labour", "unskilled labour"], 0.9)
shocked_model = shock_gtap(baseline_model, [labour_shock])

# Compare to extract outcomes
outcomes = compare_gtaps(baseline_model, shocked_model)
println("GDP change by region: ", outcomes.delta_gdp)
println("Welfare change by region: ", outcomes.ev)

shock_gtap(model::GTAP.model_container_struct, shocks::Vector{ParameterShock})

Apply parameter shocks to a GTAP economic model and compute a new equilibrium.

This function applies a vector of parameter shocks to a deep copy of the input model, runs the GTAP model to compute the new equilibrium response, and returns the modified model. The original model is left unchanged.

Arguments

• model::GTAP.model_container_struct: A GTAP model container with initialized data, sets, and parameters
• shocks::Vector{ParameterShock}: Vector of parameter shocks to apply. Each shock scales a specified parameter value by a given scale factor.

Returns

GTAP.model_container_struct: A modified copy of the input model containing:

• Updated data arrays reflecting the applied shocks
• Solved equilibrium values for endogenous variables (e.g., model.data["y"] for income)
• Original sets and parameters

To extract economic outcomes from the result, use compare_gtaps(original_model, result_model) to get GDP changes, equivalent variation, and other comparisons.

Example

# Apply a 50% reduction to skilled and unskilled labor
shocks = [ParameterShock("qe", ["skilled labor", "unskilled labor"], 0.5)]
shocked_model = shock_gtap(model, shocks)

# Compare with original to get economic outcomes
outcomes = compare_gtaps(model, shocked_model)

initial_gtap_model(datadir::String;
    roi::Union{String, Vector{String}} = ROI)::GTAP.model_container_struct

Initialize and calibrate a GTAP (Global Trade Analysis Project) economic model.

This function orchestrates the complete setup workflow for a GTAP model: loading raw trade and economic data from HAR files, aggregating regions/commodities/endowments according to regional interest (roi), building the model structure, and running calibration to ensure consistency.

Arguments

• datadir::String: Path to directory containing the three GTAP data files:
  • gsdfset.har: Set definitions (regions, commodities, endowments)
  • gsdfdat.har: Economic data (bilateral trade flows, domestic transactions, value added)
  • gsdfpar.har: Parameters (elasticities, behavioral coefficients)
• roi::Union{String, Vector{String}}: Regions of interest to preserve in aggregation. Defaults to ROI (see Helpers.ROI). All other regions are aggregated to "row" (rest of world) or "eur" for EU member states.

Returns

• model_container_struct: A calibrated GTAP model ready for scenario analysis.

Workflow

1. Loads three HAR files and converts to NamedArray dictionaries
2. Creates aggregation maps for regions, commodities, and endowments (grouping non-roi into defaults)
3. Aggregates all data/parameters/sets using these maps via aggregate_data()
4. Generates model structure from aggregated data via generate_initial_model()
5. Prepares calibration equations and data via generate_calibration_inputs()
6. Runs model in-place to enforce consistency and calibration constraints
7. Returns the calibrated model container

ROI::Vector{String}

A string vector of regions of interest. Includes the EU and some EU member states separately.

cluster_commodities(commodities::NamedArray)::NamedArray

Aggregate commodities in a NamedArray by clustering into major commodity groups.

Arguments

• commodities::NamedArray: Input array of commodity codes.

Returns

• NamedArray: New array with commodities clustered into groups: crop products, animal-derived products, extractive-industry products, processed food, manuf (manufacturing), transport, hospitality, and leisure, and other svces (services).

Example

commodities = NamedArray([1:65;], [:commodity])
result = cluster_commodities(commodities)

cluster_endowments(endowm::NamedArray)::NamedArray

Aggregate endowments in a NamedArray by clustering into major factor groups.

Arguments

• endowm::NamedArray: Input array of endowment codes.

Returns

• NamedArray: New array with endowments clustered into groups: land, skilled labor, unskilled labor, capital, and other.

Example

endowments = NamedArray([:land, :skilled, :unskilled, :capital, :other], [:endowment])
result = cluster_endowments(endowments)

cluster_named_array(x::NamedArray, preserve::Union{String, Vector{String}};
    fill_value::String = "XYZ",
    groups::Union{Nothing, Dict{String, <:AbstractVector{String}}, NamedTuple} = nothing)

Aggregate elements in a NamedArray by preserving specified indices and filling the rest with a specified value. Optionally group non-preserved elements into named categories.

Arguments

• x::NamedArray: Input NamedArray to cluster.
• preserve::Union{String, Vector{String}}: Index or indices to keep unchanged.
• fill_value::String: Value to assign to non-preserved elements. Defaults to "XYZ".
• groups::Union{Nothing, Dict, NamedTuple}: Optional mapping of group names to member lists. If provided, members of each group are set to the group name instead of fill_value. Can be a Dict{String, Vector{String}} or a NamedTuple mapping group names to their members. Defaults to nothing (no grouping, all non-preserved elements get fill_value).

Returns

• NamedArray: New array with preserved elements unchanged. Non-preserved elements are assigned to their group name (if in groups) or to fill_value otherwise.

Examples

# Basic clustering without groups
x = NamedArray(["a", "b", "c", "d"], [:item])
result = cluster_named_array(x, ["a", "c"]; fill_value = "other")
# result: ["a", "other", "c", "other"]

# Clustering with Dict groups
groups = Dict("group1" => ["b"], "group2" => ["d"])
result = cluster_named_array(x, ["a", "c"]; fill_value = "other", groups = groups)
# result: ["a", "group1", "c", "group2"]

# Clustering with NamedTuple groups
groups = (group1 = ["b"], group2 = ["d"])
result = cluster_named_array(x, ["a", "c"]; fill_value = "other", groups = groups)
# result: ["a", "group1", "c", "group2"]

cluster_regions(regions::NamedArray, preserve::Union{String, Vector{String}}
     = ["gbr", "usa", "chn", "eur"]; fill_value::String = "row")::NamedArray

Aggregate regions in a NamedArray by clustering non-preserved regions into a "row" (rest of world) category.

Arguments

• regions::NamedArray: Input array of region codes.
• preserve::Union{String, Vector{String}}: Region(s) to keep separate from aggregation. Defaults to ["gbr", "usa", "chn", "eur"].
  • If "eur" is specified, it is expanded to all EU member states (via get_eu_states()) unless it is on the preserve list, then uses the groups mechanism in cluster_named_array to map individual EU member states to "eur".
• fill_value::String: Value to assign to non-preserved, non-grouped regions. Defaults to "row" (rest of world).

Returns

• NamedArray: New array with preserved regions unchanged and all other regions aggregated into fill_value.
  • Individual EU member states are mapped to "eur" if "eur" was in the preserve list.

Example

regions = NamedArray(["gbr", "usa", "fra", "ind"], [:region])
result = cluster_regions(regions, ["gbr", "usa"])
# result: ["gbr", "usa", "row", "row"]

get_example_model()

Get a simple example model using example data from GlobalTradeAnalysisProject.jl.

shock_gtap_example(model::GTAP.model_container_struct, shocks::Vector{ParameterShock})

Run an example of passing a (labour supply) shock to a GTAP model, using the example dataset provided in GlobalTradeAnalysisProjectModelV7.jl.

This is a convenience function that delegates to shock_gtap but using example data. This function exists because the full GTAP data is used under license and cannot be shared online.

calc_consumption_avail(deaths::AbstractVector, phi::Real)

Calculate daily consumption availability from infection-avoidance behaviour.

Models the reduction in consumption (e.g., reduced shopping, dining out, travel) as exponential decay driven by daily new deaths. On days with high mortality, consumers avoid going out for discretionary purchases.

Arguments

• deaths::AbstractVector: Cumulative death count time series (length n_days)
• phi::Real: Infection-avoidance scaling parameter. Larger values mean stronger consumption reduction in response to deaths (dimensionally, 1/deaths)

Returns

Vector{Float64}: Daily consumption availability fractions, calculated as exp(-phi * new_deaths) where new_deaths = [0; diff(deaths)]

Example

cumulative_deaths = [0.0, 1.0, 3.0, 5.0, 5.0, 6.0]
phi = 0.01
avail = calc_consumption_avail(cumulative_deaths, phi)
# new_deaths: [0, 1, 2, 2, 0, 1]
# avail: [exp(0), exp(-0.01), exp(-0.02), exp(-0.02), exp(0), exp(-0.01)]
#      ≈ [1.0, 0.99, 0.98, 0.98, 1.0, 0.99]

calc_indivs(data::DataFrame, scaling::Dict)

Calculate person-equivalents over time from epidemiological compartment data.

This function converts epidemiological model outputs (such as compartmental prevalence counts from disease states like infected, hospitalized, or deceased) into weighted person-equivalents by applying scaling factors and summing across compartments.

Arguments

• data::DataFrame: Input DataFrame containing epidemiological compartment data with rows representing time points and columns representing disease compartments.
• scaling::Dict: Dictionary mapping column names (as String or Symbol) to scaling factors. Scaling factors are applied to weight the relative impact of each compartment. For example, Dict("infected" => 0.1, "hospitalized" => 1.0) represents that hospitalized individuals have 10x the impact of infected individuals.

Returns

Vector: A vector of person-equivalents at each time point, calculated as the row-wise sum of scaled compartment values.

Example

df = DataFrame(
    infected = [100, 150, 200],
    hospitalized = [10, 15, 20],
    deceased = [1, 2, 3]
)
scaling = Dict("infected" => 0.1, "hospitalized" => 1.0, "deceased" => 5.0)
person_equiv = calc_indivs(df, scaling)
# Returns: [101.1, 151.7, 203.3]

Details

The function performs the following steps:

1. Deep copies the input DataFrame to avoid modifying the original
2. Scales specified columns in-place by their corresponding scaling factors
3. Selects only the scaled columns (drops unscaled columns)
4. Sums across columns for each row to obtain person-equivalents
5. Returns the vector of row sums

Only columns specified in the scaling dictionary are included in the final calculation.

calc_labour_avail(epi_data::DataFrame, scaling::Dict, N_work::Real)

Calculate daily labour availability from epidemiological compartment data.

Computes the fraction of the workforce still available at each time point by:

1. Calculating weighted person-equivalents (absent workers) using calc_indivs
2. Normalizing by total workforce
3. Returning availability as 1.0 - person_equiv / N_work

Arguments

• epi_data::DataFrame: Epidemiological time series with compartment columns (e.g., "mild", "severe", "hospitalized", "deceased")
• scaling::Dict: Mapping of compartment column names to their productivity impact weights. For example, mildly symptomatic workers might have weight 0.36 (i.e., 1 - 0.64 productivity)
• N_work::Real: Total workforce size

Returns

Vector{Float64}: Daily labour availability fractions (typically in [0.0, 1.0], though can exceed 1.0 or be negative if data is inconsistent)

Example

df = DataFrame(
    mild = [100, 150, 200],
    severe = [10, 15, 20],
    hospitalized = [5, 7, 10],
    deceased = [1, 2, 3]
)
scaling = Dict("mild" => 0.36, "severe" => 1.0, "hospitalized" => 1.0, "deceased" => 1.0)
N_work = 10000.0
avail = calc_labour_avail(df, scaling, N_work)
# Person-equivalents: [136.6, 207.8, 283]
# Availability: [1 - 136.6/10000, 1 - 207.8/10000, ...] ≈ [0.9863, 0.9792, 0.9717]

integrate_shock(t::AbstractVector, avail::AbstractVector)

Compute time-weighted average of a daily availability time series using trapezoidal integration.

Integrates the availability curve over time and normalizes by the total time span to produce a single scalar shock value (typically in [0.0, 1.0]) suitable for use as a ParameterShock scale factor.

Arguments

• t::AbstractVector: Time vector (numeric, typically day-of-year or date-as-integer). Must be strictly increasing.
• avail::AbstractVector: Daily availability fractions, same length as t

Returns

Float64: Time-weighted average availability, calculated as:

area_under_curve / (t[end] - (t[1] - 1.0))

where area_under_curve is computed using the trapezoidal rule.

The normalization t[end] - (t[1] - 1.0) represents the total time span including the first day.

Example

t = [1.0, 2.0, 3.0, 4.0]
avail = [1.0, 0.95, 0.90, 0.90]
shock = integrate_shock(t, avail)
# Trapezoidal area ≈ (1.0 + 0.95)/2 + (0.95 + 0.90)/2 + (0.90 + 0.90)/2 = 2.8
# shock ≈ 2.8 / (4.0 - 0.0) = 0.7