# packages
using DataFrames, CSV, Statistics, LinearAlgebra, Plots, StatsPlots, GLM, Random, Optim

cd("C:/Users/Eric/Documents/GitHub/Econ-543/PS 2/data")  # adjust path as needed
consumer_census = DataFrame(CSV.File("consumer_census.csv"))
product_data = DataFrame(CSV.File("product_data_ps2.csv"))

2500×9 DataFrame

2475 rows omitted

Row	product_ID	price	sugar	caffeine	corn_syrup_price	caffeine_extract_price	quantity	t	market_share
	Int64	Float64	Float64	Float64	Float64	Float64	Int64	Int64	Float64
1	1	2.02133	1.58086	6.07446	0.199022	0.264899	223	1	0.0223
2	2	1.98025	3.03315	5.2606	0.199022	0.264899	39	1	0.0039
3	3	1.65556	1.71888	4.68573	0.199022	0.264899	83	1	0.0083
4	4	1.56647	2.63786	4.25037	0.199022	0.264899	16	1	0.0016
5	5	1.51763	1.47621	4.1378	0.199022	0.264899	141	1	0.0141
6	6	1.93513	1.01706	7.04804	0.199022	0.264899	92	1	0.0092
7	7	2.01592	1.83851	5.93992	0.199022	0.264899	151	1	0.0151
8	8	2.39419	2.46793	7.00869	0.199022	0.264899	322	1	0.0322
9	9	1.75108	2.56961	4.57863	0.199022	0.264899	47	1	0.0047
10	10	1.8627	2.83581	4.64917	0.199022	0.264899	71	1	0.0071
11	11	1.31229	1.10588	3.97262	0.199022	0.264899	69	1	0.0069
12	12	1.7945	3.01562	4.74948	0.199022	0.264899	21	1	0.0021
13	13	1.75079	1.64843	5.47186	0.199022	0.264899	44	1	0.0044
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
2489	39	2.64074	7.42992	4.67266	0.210564	0.255062	272	50	0.0272
2490	40	2.35862	3.32281	6.43918	0.210564	0.255062	101	50	0.0101
2491	41	3.00301	6.36023	6.13807	0.210564	0.255062	1073	50	0.1073
2492	42	2.5826	5.73171	5.31855	0.210564	0.255062	166	50	0.0166
2493	43	2.36855	5.79571	4.55279	0.210564	0.255062	112	50	0.0112
2494	44	1.96981	3.88361	5.00891	0.210564	0.255062	14	50	0.0014
2495	45	2.79746	5.64465	6.53525	0.210564	0.255062	162	50	0.0162
2496	46	2.36863	3.69903	5.59846	0.210564	0.255062	252	50	0.0252
2497	47	1.71946	4.03176	3.96157	0.210564	0.255062	12	50	0.0012
2498	48	2.52919	7.77044	3.60625	0.210564	0.255062	1158	50	0.1158
2499	49	1.67455	3.75178	3.64218	0.210564	0.255062	30	50	0.003
2500	50	2.01727	3.82539	4.59095	0.210564	0.255062	63	50	0.0063

Introduction

This homework assignment will give you practice working with the basics of the random coefficients multinomial logit demand model (Berry 1994). Compared to Problem set 1, we will now introduce three new sources of consumer heterogeneity:

1. Consumer demographic characteristics (income)
2. Unobserved preference for sugar
3. Unobserved preference for caffeine

The file product_data_ps2.csv contains information on 50 different soft drinks, for 50 different markets. This dataset contains the following variables: product_ID (the product ID), type (a variable indicating if the product is a diet or a regular soda), price (the product's price during period $t$, in USD), quantity (the number of cans sold in market $t$), sugar (the grams of sugar contained in product $j$ at time $t$), caffeine (the milligrams of caffeine contained in product $j$ at time $t$), caffeine_extract_price (the price of caffeine at time $t$), corn_syrup_price (the price of high fructose corn syrup at time $t$), and $t$ (the observation's time period). Throughout this problem set, assume that the market size is $M_t = \bar{M} = 10,000$ for all markets.

Let $\mathcal{J}_t = \{1,...,J_t\}$ denote the set of all products available in market $t$. As usual, let $j = 0$ denote the outside good (not consuming any product).

Consumer $i$'s indirect utility from consuming soda $j$ in market $t$ is given by the following expression:

$$u_{ijt} = \gamma + \alpha_i p_{jt} + \beta_1 sug_{jt} + \beta_2 caf_{jt} + \xi_{jt} + \varepsilon_{ijt}$$

with $$\alpha_i = p_{jt} (\hat{\alpha} + \alpha^y y_{it}),$$ $$\beta_{1it} = \bar{\beta}_1 + \nu_{1it} \sigma_1, \quad \nu_{1it} \sim N(0,1)$$ $$\beta_{2it} = \bar{\beta}_2 + \nu_{2it} \sigma_2, \quad \nu_{2it} \sim N(0,1)$$

Problem 1 (a)

The file consumer_census.csv contains micro-level data on households' income for each market $t$ in our sample. Let $\bar{y}_t$ denote the average income in market $t$, let $\sigma_{yt}$ denote the variance of $y_{it}$ in market $t$, and let $\bar{p}_t = \sum_j p_{jt} s_{jt}$ denote the average price paid in market $t$.

1. Compute $\bar{y}_t$, $\sigma_{yt}$, and $\bar{p}_t$ for each market.
2. Compute $\text{cor}(\bar{y}_t, \bar{p}_t)$. Do consumers in wealthier markets buy more expensive products on average?

# Problem 1(a) solution
mean_income = combine(groupby(consumer_census, :t), 
    :income => mean => :mean_income)

var_income = combine(groupby(consumer_census, :t), 
    :income => var => :var_income)

product_data[!, :market_share] .= product_data.quantity ./ 10000

grouped_product = groupby(product_data, :t)
p_bar = combine(grouped_product) do sdf
    DataFrame(p_bar = mean(sdf.price .* sdf.market_share))
end

market_data = innerjoin(mean_income, p_bar, on = :t)

correlation = cor(market_data.mean_income, market_data.p_bar)
println(correlation)

0.15430274953743087

Problem 1(b)

Write a function that takes the observed distribution of income of each market, the number of consumers to be simulated ($S$), and a random seed (for replicability) as inputs, and uses pseudo-Monte Carlo draws to return a simulated sample of $S$ consumers for each market $t$. For each consumer $i$ and each market $t$, the output of this function must contain $y_{it}$, $\nu_{1it}$, and $\nu_{2it}$. Use this function to draw a sample of 500 consumers for each market $t$.

function simulate_consumers(consumer_census::DataFrame, S::Int, seed::Int)
    Random.seed!(seed)

    # Split the data by the 't' column
    grouped_census = groupby(consumer_census, :t)

    simulated_consumers = DataFrame()
    
    tt = 1
    for sub_df in grouped_census
        income_sample = rand(sub_df.income, S)
        
        simulated_consumers_current_t = DataFrame(
            consumer_ID = 1:S,
            income = income_sample,
            t = fill(tt, S),
            nu_1 = randn(S),
            nu_2 = randn(S)
        )
        
        append!(simulated_consumers, simulated_consumers_current_t)
        
        tt += 1
    end

    return simulated_consumers
end

sample = simulate_consumers(consumer_census, 500, 123)

25000×5 DataFrame

24975 rows omitted

Row	consumer_ID	income	t	nu_1	nu_2
	Int64	Float64	Int64	Float64	Float64
1	1	5.68779	1	-0.563549	-0.264958
2	2	6.19198	1	0.227326	0.230161
3	3	5.85358	1	0.0102916	0.249837
4	4	5.56637	1	-0.760516	-0.283196
5	5	5.65982	1	1.39799	0.501889
6	6	6.10602	1	-0.343564	0.385318
7	7	6.89251	1	-1.54224	-2.85127
8	8	6.50529	1	-1.3	-0.259393
9	9	5.55876	1	-1.05767	1.12154
10	10	6.19681	1	-0.549124	1.62369
11	11	6.44758	1	1.12477	-0.372927
12	12	5.62771	1	1.23201	1.18332
13	13	6.26405	1	0.0875616	1.65032
⋮	⋮	⋮	⋮	⋮	⋮
24989	489	6.01369	50	0.0896844	0.11497
24990	490	5.62916	50	2.30871	-0.448778
24991	491	6.04779	50	0.991243	-0.460348
24992	492	5.32704	50	-1.5506	-0.856875
24993	493	5.89291	50	-0.18464	-1.16149
24994	494	5.62934	50	0.840827	0.8664
24995	495	6.07691	50	0.00597359	-1.69253
24996	496	7.33057	50	0.816041	1.5106
24997	497	5.82752	50	0.078145	0.506398
24998	498	6.89597	50	0.334505	-1.20558
24999	499	6.34592	50	-0.1372	0.303473
25000	500	6.02839	50	-0.392797	-1.46894

Problem 1(c)

Given our distributional assumptions over $\varepsilon_{ijt}$, the probability that consumer $i$ in market $t$, with heterogeneity given by $y_{it}$, $\nu_{1it}$, and $\nu_{2it}$ purchases product $j$ is given by:

$s_{ijt}(\delta_t, \theta_2) = \text{Pr}(u_{ijt} \geq u_{ikt} \forall k \in \mathcal{J}_t) = \frac{\exp(\delta_{jt} + \mu_{ijt}(\theta_2))}{1 + \sum_k \exp(\delta_{kt} + \mu_{ikt}(\theta_2))}$

Note that when computing $s_{ijt}$, it is important to only consider the goods available in $is market ($\mathcal{J}_t$).

Let $\mathcal{S}$ denote the set of simulated consumers you obtained in Problem 1(b). Write a function that, given the parameters $\delta_t = \{\delta_{1t}, \ldots, \delta_{J_t,t}\}$ and $\theta_2 = \{\alpha, \sigma_1, \sigma_2\}$ computes the vector of individual choice probabilities $s_{ijt} = \text{Pr}(u_{ijt} \geq u_{ikt} \forall k \in \mathcal{J}_t)$ for all $i \in \mathcal{S}$. You will use this function for part (d).

function s_ijt(simulated_consumers::DataFrame, delta_table::DataFrame, product_data::DataFrame, nonlinear_param::Dict{Symbol, Float64})
    
    alpha = nonlinear_param[:alpha]
    sigma_1 = nonlinear_param[:sigma_1]
    sigma_2 = nonlinear_param[:sigma_2]
    
    product_variables = innerjoin(delta_table, product_data, on = [:product_ID, :t])
    
    grouped_simulated = groupby(simulated_consumers, :t)
    grouped_product = groupby(product_variables, :t)

    individual_shares = DataFrame()
    
    for tt in 1:length(grouped_product)
        
        data = sort(grouped_product[tt], :product_ID)
        n_j = nrow(data)
        
        consumer_census_t = sort(grouped_simulated[tt], :consumer_ID)
        
        income_data = transpose(repeat(consumer_census_t.income, 1, n_j))
        sugar_taste = transpose(repeat(consumer_census_t.nu_1, 1, n_j))
        caffeine_taste = transpose(repeat(consumer_census_t.nu_2, 1, n_j))
        
        mu_ijt = alpha .* data.price .* income_data .+ sigma_1 .* sugar_taste .* data.sugar .+ sigma_2 .* caffeine_taste .* data.caffeine
        utility = mu_ijt .+ data.delta
        
        # Calculate utility from outside option
        outside_option = zeros(size(utility, 2))     
       
         # Calculate inclusive values (log-sum-exp over rows)
        inclusive_values = log.(sum(exp.(vcat(utility, outside_option')), dims=1))
        
         current_shares = exp.(utility .- inclusive_values)
        
        current_shares_df = DataFrame(current_shares, :auto)
        rename!(current_shares_df, string.(consumer_census_t.consumer_ID))
        current_shares_df.product_ID = data.product_ID
        current_shares_df.t = fill(tt, nrow(current_shares_df))
        
        append!(individual_shares, current_shares_df)
        
    end
    
    return individual_shares
    
end

s_ijt (generic function with 1 method)

Problem 1(d)

Suppose the linear parameters are $\hat{\alpha} = -1$, $\bar{\beta}_1 = 0.5$, $\bar{\beta}_2 = 0.5$, and $\gamma = -1.5$. Compute $\delta_{jt}$ for each product and each market, assuming $\xi_{jt} = 0 \forall j, t$. Call the vector of $\delta$'s you get $\delta^0$. Let $\theta_2^0 = (\alpha^y, \sigma_1^y, \sigma_2^y) = (0, 0, 0)$.

Compute the conditional shares $s_{ijt}(\delta_t^0, \theta_2^0)$ for each consumer. Let $\theta_2^1 = (\alpha^y, \sigma_1^1, \sigma_2^1) = (0.05, 1, 1)$. Compute the conditional shares $s_{ijt}(\delta_t^0, \theta_2^1)$. Comment on how the conditional probabilities change if we use $\theta_2^1$ or $\theta_2^0$.

function get_delta(product_data::DataFrame, linear_parameters::Dict{Symbol, Float64}; xi::Bool=false, xi_hat::DataFrame=DataFrame())
    
    product_variables = copy(product_data)
    
    alpha_bar = linear_parameters[:alpha_bar]
    beta_1 = linear_parameters[:beta_1]
    beta_2 = linear_parameters[:beta_2]
    intercept = linear_parameters[:intercept]
    
    if xi
        product_variables = innerjoin(product_variables, xi_hat, on = [:product_ID, :t])
        product_variables.delta = intercept .+ alpha_bar .* product_variables.price .+ beta_1 .* product_variables.sugar .+ beta_2 .* product_variables.caffeine .+ product_variables.xi
    else
        product_variables.delta = intercept .+ alpha_bar .* product_variables.price .+ beta_1 .* product_variables.sugar .+ beta_2 .* product_variables.caffeine
    end
    
    delta_table = select(product_variables, :product_ID, :delta, :t)
    
    return delta_table
    
end

get_delta (generic function with 1 method)

linear_parameters = Dict(
    :alpha_bar => -1.0,
    :beta_1 => 0.5,
    :beta_2 => 0.5,
    :intercept => -1.5
)

delta_table = get_delta(product_data, linear_parameters)

theta_2_prime = Dict(
    :alpha => 0.0,
    :sigma_1 => 0.0,
    :sigma_2 => 0.0
)

theta_2_double_prime = Dict(
    :alpha => 0.0,
    :sigma_1 => 1.0,
    :sigma_2 => 1.0
)

individual_shares_prime = s_ijt(sample, delta_table, product_data, theta_2_prime)

individual_shares_prime_prime = s_ijt(sample, delta_table, product_data, theta_2_double_prime)

2500×502 DataFrame

402 columns and 2475 rows omitted

Row	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	⋯
	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	Float64	⋯
1	0.0209467	0.0069844	0.0129582	0.0217894	7.40465e-5	0.0275348	3.56659e-9	0.0193389	0.0435752	0.0284815	0.000100953	0.000253598	0.0116217	0.000156998	0.0156112	0.00984782	0.000237556	2.85719e-6	0.00241107	0.0403777	0.000913044	0.000279508	0.000880872	0.00368382	0.0377294	3.58705e-5	5.02292e-5	0.00206006	0.00593158	0.00188098	0.0126545	0.0111695	0.00136373	0.0144161	0.013461	0.00536722	0.0264436	0.0283992	0.0136363	5.46991e-7	0.0223877	0.000572919	0.026502	3.87828e-5	0.0346481	0.0396781	0.0288478	0.00121783	0.0001299	0.000127429	0.000207222	0.009626	0.00128878	0.00256755	0.000921898	2.03344e-5	0.00740714	0.0205036	5.51174e-5	0.0072492	0.00112714	0.00775113	0.0188332	0.00274489	0.00272472	0.00202818	0.000179995	0.00901417	0.0242006	0.0515012	0.0349647	0.00127058	0.0217814	0.00413354	2.86526e-6	0.0501407	0.0182311	0.0482817	9.56796e-5	0.0237622	8.94103e-6	1.85517e-5	0.000190029	0.0258612	0.0178465	0.000182225	0.00616494	8.27957e-6	0.00114678	0.00704081	6.71712e-5	0.0380566	0.00518514	0.000741895	0.000159233	0.0096115	0.00336346	0.000701299	0.0140449	0.00313568	⋯
2	0.0164361	0.0115514	0.0153888	0.013036	0.00053744	0.0175175	5.54373e-9	0.00518387	0.00539767	0.00490674	0.00100422	0.00083067	0.00493922	0.00131936	0.0145978	0.00454311	0.00206724	8.134e-6	0.0023997	0.00568685	0.00130783	0.00198622	0.000966051	0.00153081	0.00197377	0.000603647	0.000562636	0.000580499	0.01095	0.00548046	0.0174062	0.00388591	0.00554653	0.0166019	0.0173209	0.012101	0.0184275	0.0033145	0.0178845	3.9493e-6	0.0104578	0.00162297	0.0112698	6.4653e-5	0.0115368	0.00925235	0.0154645	0.00266935	4.5028e-5	0.00143087	0.00160388	0.00612674	0.0047186	0.00800042	0.00427422	0.000424441	0.0123523	0.0172276	0.000364437	0.00647448	0.0053509	0.00870207	0.0178129	0.00726301	0.002891	0.000444594	0.00160626	0.0155156	0.0150807	0.00503333	0.00828992	0.00429762	0.0192788	0.000768173	6.11529e-5	0.00390757	0.0183774	0.0103912	0.001065	0.0010992	0.000139427	0.000198487	0.00122628	0.00816708	0.0179635	0.00158721	0.0117477	1.88694e-5	0.00385818	0.01262	6.99616e-6	0.0116784	0.00567533	0.00382754	0.00165847	0.010574	0.00985757	0.00332335	0.0138405	0.0067502	⋯
3	0.0215974	0.00403837	0.00707529	0.0224246	3.4504e-5	0.0118621	1.16608e-7	0.0178739	0.0061191	0.00213618	0.00015266	4.48308e-5	0.000917159	0.000211994	0.0192883	0.0127583	0.000651858	3.17753e-5	0.000211378	0.022417	7.74705e-5	0.00145906	0.00246826	0.000122873	0.0164103	0.000127721	5.68726e-5	0.00327054	0.00380498	0.00121951	0.0140277	0.0129819	0.00178264	0.00864473	0.0104558	0.00602222	0.0215813	0.00161141	0.0128246	1.09813e-5	0.00387254	0.00013182	0.00498475	0.000224487	0.00715276	0.024621	0.0242007	0.000282122	0.000373703	0.000481324	0.000283708	0.00115492	0.00109573	0.00300734	0.00116708	0.000210399	0.00473299	0.00965261	0.000442422	0.00119378	0.00227071	0.0130966	0.0105004	0.00201476	0.000288585	0.00299799	0.0010405	0.0113056	0.00765983	0.0144076	0.024086	0.00419631	0.0146773	0.00503418	1.99834e-6	0.0170317	0.011754	0.0162143	0.000207036	0.000727647	6.30015e-6	8.29419e-6	0.00113838	0.00315391	0.0108118	0.00106325	0.00458496	7.04705e-5	0.000685388	0.0119647	0.000152639	0.0251092	0.000932872	0.00106039	0.00062817	0.00302777	0.00607695	0.000696674	0.00538925	0.00151788	⋯
4	0.0201044	0.00626794	0.00891913	0.017557	0.000139519	0.0101835	1.3615e-7	0.0084359	0.00197797	0.000885503	0.00070283	0.000115679	0.000674632	0.000882404	0.019863	0.00834383	0.0026575	5.82962e-5	0.000266006	0.00723542	0.000123252	0.00502626	0.00267769	9.3877e-5	0.00285621	0.000780918	0.000287464	0.00154213	0.0062953	0.0026985	0.0186049	0.007137	0.00467744	0.0106251	0.0135794	0.0109425	0.0188896	0.000529818	0.0166654	3.52136e-5	0.0028849	0.000305097	0.00347953	0.000304477	0.00424679	0.0108913	0.0178855	0.000554091	0.000193995	0.00226476	0.00111953	0.00107183	0.00275954	0.00670102	0.00333588	0.00138351	0.0073403	0.00984421	0.001418	0.00134926	0.00640591	0.0148717	0.0114317	0.00416074	0.000373987	0.00120941	0.00411992	0.0171737	0.00661215	0.00383901	0.0107002	0.00925213	0.0151411	0.00184919	1.57004e-5	0.00388136	0.0132107	0.00707348	0.00100206	0.00013755	4.04848e-5	4.30764e-5	0.00365677	0.00186913	0.0121863	0.004143	0.00767181	0.000114003	0.00167075	0.0183102	3.72992e-5	0.0131868	0.00119414	0.00322112	0.00280971	0.00375836	0.012683	0.00204927	0.00616126	0.00281802	⋯
5	0.0221346	0.0026044	0.00475813	0.0243496	1.44323e-5	0.0080706	6.25144e-7	0.0218361	0.00330749	0.000775108	0.000110194	1.34395e-5	0.000281016	0.000149616	0.020914	0.0165256	0.000660893	8.51457e-5	5.68987e-5	0.0241351	1.89458e-5	0.00216157	0.00396276	2.41437e-5	0.0192793	0.000128585	3.44059e-5	0.00524092	0.00252949	0.000738478	0.0132205	0.0168913	0.00144762	0.00614559	0.00830493	0.00513174	0.0200892	0.000550016	0.0112086	3.27296e-5	0.00175038	4.67886e-5	0.00240209	0.000480523	0.00381522	0.0251845	0.0241898	0.000106237	0.000780679	0.000545516	0.000205306	0.000403386	0.000727375	0.00244329	0.000905958	0.00035263	0.00320994	0.00647609	0.000833242	0.000457026	0.00223443	0.016069	0.00749975	0.00132908	8.45502e-5	0.00488733	0.00154743	0.0108249	0.00445747	0.0117437	0.0259236	0.00575084	0.0117394	0.00768847	8.16564e-7	0.0161337	0.00893187	0.0122522	0.000175817	0.000220023	2.7703e-6	3.12145e-6	0.00185027	0.00129318	0.00795602	0.00159734	0.00326984	0.00017236	0.000385586	0.0133049	0.000365868	0.0249875	0.000357886	0.000856855	0.000740231	0.00155394	0.00627211	0.000474316	0.00313966	0.000842731	⋯
6	0.0297506	0.0102851	0.0219828	0.0339737	7.34234e-5	0.0650652	7.09124e-10	0.0418299	0.315388	0.25233	4.98261e-5	0.000536092	0.0737981	8.68628e-5	0.0190755	0.0149467	9.10146e-5	7.93887e-7	0.0108698	0.140404	0.00369644	8.31089e-5	0.000676141	0.0363042	0.207755	9.63398e-6	2.75264e-5	0.00333979	0.00798679	0.00215662	0.0143165	0.019707	0.000981909	0.0235788	0.0186203	0.00510047	0.043971	0.330425	0.0170266	8.23224e-8	0.0928745	0.00111498	0.108847	1.79402e-5	0.147611	0.109898	0.0517321	0.00257759	0.000145661	3.81344e-5	0.000117093	0.0431771	0.0011933	0.0021388	0.000645109	2.97246e-6	0.0103451	0.0421533	1.3536e-5	0.0247186	0.000619792	0.00760667	0.0341646	0.00305589	0.0109321	0.0037123	4.70554e-5	0.00888089	0.0672195	0.286003	0.0913801	0.000611995	0.0367525	0.00881996	1.60746e-6	0.27287	0.0300887	0.186978	3.79743e-5	0.516246	5.5532e-6	1.64565e-5	5.44935e-5	0.146066	0.0304026	4.78092e-5	0.00762484	2.83778e-6	0.00130413	0.00587505	0.00012768	0.0930875	0.0157815	0.000471042	4.78339e-5	0.022051	0.00230198	0.000550073	0.0303022	0.00461448	⋯
7	0.0200735	0.00767724	0.0134333	0.0198965	0.000106096	0.0255864	3.76152e-9	0.0153181	0.03051	0.0212486	0.000151655	0.000317653	0.0101795	0.000229063	0.0154142	0.00857064	0.000346907	3.38362e-6	0.00245557	0.0286281	0.000992443	0.000391805	0.000889321	0.00323464	0.0224653	5.87969e-5	7.72365e-5	0.00163966	0.00664373	0.00228503	0.0133949	0.0092519	0.00174886	0.0148526	0.0141179	0.00620343	0.0248548	0.0198089	0.0143299	7.61187e-7	0.0198252	0.000697921	0.0230646	4.19578e-5	0.0288483	0.03075	0.0258714	0.00141721	0.000106799	0.000194265	0.000297911	0.0090319	0.00162732	0.00314263	0.00121044	3.43582e-5	0.0081478	0.0200073	7.603e-5	0.00720773	0.00148099	0.00788811	0.0187451	0.00327474	0.00280282	0.00154449	0.000262558	0.00992119	0.0224573	0.0343904	0.0271562	0.00156581	0.0213943	0.00306097	4.96002e-6	0.0321116	0.0183329	0.0370575	0.000146246	0.0141208	1.46351e-5	2.84927e-5	0.000261571	0.0214147	0.0179498	0.000264646	0.00693652	9.44308e-6	0.00143015	0.0077876	4.46573e-5	0.0309574	0.00534209	0.000991915	0.000239441	0.00986982	0.00405888	0.000926079	0.0141214	0.00361726	⋯
8	0.0169851	0.0181406	0.0282772	0.0145853	0.000700287	0.0498219	1.08348e-10	0.00820209	0.0832481	0.136569	0.000330898	0.00391247	0.100497	0.000516191	0.0117775	0.00425825	0.000369304	4.33971e-7	0.0341429	0.0188592	0.0173807	0.000186337	0.000304686	0.080247	0.0113351	6.72245e-5	0.000245768	0.000501686	0.0148791	0.00646189	0.0141596	0.00444444	0.00276581	0.0277827	0.0211727	0.00843641	0.0254252	0.139802	0.0176462	8.47842e-8	0.0879491	0.00597598	0.0888475	8.19656e-6	0.0881071	0.0235321	0.0223152	0.0104972	1.91707e-5	0.000167693	0.000632792	0.0703201	0.00387744	0.0048526	0.00212623	1.38737e-5	0.0173525	0.0410799	2.18206e-5	0.047747	0.0015923	0.00473537	0.034141	0.00768734	0.0328586	0.000446254	0.000126327	0.0103463	0.0604046	0.0391471	0.0187226	0.000821356	0.0306224	0.000998533	3.77123e-5	0.0262309	0.0295165	0.0528495	0.000229838	0.118379	9.29935e-5	0.000242209	0.000101754	0.112475	0.0308923	0.000124423	0.0134611	1.43639e-6	0.00477166	0.0059714	5.43113e-6	0.025684	0.0358574	0.00161771	0.000189382	0.0361969	0.00379116	0.00214063	0.0394203	0.0119367	⋯
9	0.0181344	0.00630208	0.00916049	0.0159552	0.00014159	0.011202	5.61736e-8	0.00801614	0.00290895	0.00148328	0.000545299	0.000148498	0.00109151	0.000705065	0.0171359	0.00725504	0.00181358	3.00538e-5	0.000411293	0.00782543	0.000190146	0.00299108	0.00191944	0.000183859	0.0033454	0.000494549	0.000237038	0.00129847	0.00616823	0.00260674	0.0162762	0.00640324	0.00383595	0.0106973	0.0128742	0.00942552	0.0179981	0.000953561	0.015138	1.55504e-5	0.0039143	0.000370876	0.00464576	0.000183894	0.00559005	0.0113436	0.0170538	0.000677034	0.000142452	0.0014306	0.000893707	0.00156421	0.00249574	0.00566961	0.00274585	0.000688179	0.00721738	0.0105419	0.000768679	0.00181983	0.00475508	0.0119151	0.0117562	0.00392024	0.000553653	0.0010413	0.0023829	0.0145345	0.00779081	0.00487769	0.0108853	0.00621106	0.0149577	0.00166897	1.41714e-5	0.00475215	0.0131343	0.00845622	0.000716137	0.000290644	3.67577e-5	4.36755e-5	0.00211615	0.00276787	0.0122843	0.00239282	0.0072691	6.27152e-5	0.00163711	0.0144744	2.9771e-5	0.013453	0.00157107	0.00257402	0.00176129	0.004368	0.00976555	0.0017787	0.00687462	0.00292682	⋯
10	0.0162137	0.00720213	0.00989483	0.0135192	0.000225254	0.011118	3.22505e-8	0.00589354	0.00251459	0.00152099	0.000758381	0.000237163	0.00132847	0.000966342	0.015392	0.0056163	0.0022226	2.33198e-5	0.000582701	0.00551221	0.000291247	0.00316104	0.00155158	0.000257353	0.0020153	0.000661664	0.00035432	0.000880644	0.00708095	0.00327142	0.0160778	0.00477195	0.00457454	0.0113053	0.0133199	0.010253	0.0163437	0.000948122	0.0152259	1.31555e-5	0.00425653	0.000550573	0.00488919	0.00014309	0.00551644	0.00867156	0.01463	0.000953909	9.11381e-5	0.00174842	0.00120219	0.00192345	0.00313241	0.0065308	0.0033792	0.000803791	0.00812456	0.0108958	0.000745079	0.00227031	0.00543924	0.0105324	0.0120932	0.004727	0.000769682	0.000679849	0.00259055	0.0147058	0.00810876	0.00358335	0.00820393	0.00611292	0.01471	0.00108916	2.61232e-5	0.0032304	0.0133544	0.00705113	0.000959183	0.000266494	6.35421e-5	7.55892e-5	0.00214383	0.00295194	0.0126173	0.00258476	0.00819211	4.86869e-5	0.00214331	0.0140312	1.55941e-5	0.0107512	0.00201098	0.00317076	0.00210518	0.00509538	0.0103341	0.00228968	0.00759036	0.00363169	⋯
11	0.0267676	0.00216535	0.00427319	0.0317701	7.43672e-6	0.00808006	1.66518e-6	0.0346546	0.00381984	0.000682487	7.25944e-5	6.58018e-6	0.000194604	0.000100092	0.0253412	0.0247703	0.000535557	0.000140742	3.19363e-5	0.0397323	9.50757e-6	0.00222658	0.00579344	1.31728e-5	0.0391521	9.31636e-5	2.01894e-5	0.00956426	0.00209995	0.000540675	0.0139309	0.0265828	0.00117748	0.00571759	0.00806002	0.00472031	0.023506	0.0004993	0.0114241	4.92183e-5	0.00147517	2.57224e-5	0.00212884	0.000767211	0.00371982	0.0371201	0.0307168	6.29078e-5	0.00159121	0.000450512	0.000141127	0.000280793	0.000539312	0.00207661	0.000703739	0.00032559	0.00273931	0.0061476	0.00098987	0.000316374	0.00196752	0.0202126	0.0072247	0.00103578	4.91811e-5	0.00938118	0.00153877	0.0110901	0.00411227	0.0176381	0.0391772	0.00641235	0.0121599	0.0146334	3.48106e-7	0.0272939	0.00880739	0.015514	0.000124242	0.000217605	1.29405e-6	1.42485e-6	0.00203659	0.00110198	0.00771076	0.00160354	0.00280589	0.000281145	0.000264885	0.014691	0.000978885	0.0346849	0.000239895	0.000668505	0.00063218	0.00121963	0.00613475	0.000341796	0.0026885	0.000618432	⋯
12	0.0175703	0.00945557	0.0125183	0.0141147	0.000375107	0.0133792	2.26115e-8	0.00559764	0.00286537	0.00199724	0.00110134	0.000410447	0.00196126	0.00139499	0.016622	0.00546034	0.00289187	2.09748e-5	0.000952286	0.00517866	0.000503086	0.00362364	0.00151083	0.000432227	0.00171382	0.000906238	0.000546807	0.000775184	0.00928986	0.00456427	0.0186223	0.00453664	0.00600251	0.0140388	0.016143	0.0126925	0.0180585	0.00123125	0.017967	1.22936e-5	0.00564496	0.000896906	0.00631858	0.0001321	0.00679424	0.0085878	0.015523	0.0014982	7.53641e-5	0.00224389	0.00171216	0.00281605	0.00432703	0.00839945	0.00453738	0.000963063	0.0105164	0.0134423	0.000791649	0.00331831	0.00682361	0.0111189	0.0147777	0.00638508	0.00123279	0.000584359	0.00301719	0.0172405	0.010201	0.00355519	0.00801397	0.00678317	0.0173011	0.000942614	4.74764e-5	0.00301507	0.0161006	0.0075877	0.00133915	0.000334674	0.0001104	0.000133666	0.00237435	0.00391255	0.0153547	0.00299608	0.0105529	4.4241e-5	0.00308876	0.0157893	1.13205e-5	0.0109562	0.00298242	0.00424125	0.00265603	0.00693974	0.0123313	0.00320056	0.00987813	0.00508172	⋯
13	0.0237265	0.0061927	0.0111886	0.0246238	6.03249e-5	0.0213931	1.79681e-8	0.0207721	0.0207016	0.0103475	0.000136785	0.000135498	0.00433766	0.000202407	0.0191344	0.0122707	0.00041572	9.07784e-6	0.000946094	0.0346979	0.000353985	0.000644481	0.00153753	0.000948775	0.0290259	7.04852e-5	6.01244e-5	0.00279402	0.00550316	0.00175762	0.0148346	0.0132609	0.00172674	0.0129579	0.013534	0.00634231	0.0270662	0.00912953	0.014889	2.25297e-6	0.0117217	0.000342622	0.014378	9.27368e-5	0.0195395	0.0358713	0.0298354	0.000729394	0.000227848	0.000256265	0.00026863	0.00431281	0.0013551	0.00309627	0.00115214	6.33324e-5	0.00685669	0.0165765	0.000152898	0.00372953	0.00172086	0.0108825	0.0163875	0.00270422	0.00116027	0.00266373	0.000434219	0.0111582	0.0164631	0.032878	0.0330681	0.00239549	0.0205564	0.00498976	2.79384e-6	0.0347558	0.0168938	0.0334933	0.000151366	0.00582671	8.74143e-6	1.48792e-5	0.000464674	0.0116229	0.0160999	0.000441305	0.00609685	2.34148e-5	0.0010355	0.00993561	0.000105726	0.0353809	0.00277438	0.000977512	0.000326136	0.00654593	0.00488673	0.000789653	0.010407	0.00257923	⋯
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋱
2489	0.186948	3.22556e-5	0.0124551	0.000506988	7.9098e-6	0.345574	0.0133567	6.48624e-5	0.230843	0.210642	3.60048e-7	0.230781	9.25411e-6	0.086963	4.03705e-5	0.000223891	0.143574	0.116007	0.198127	0.0522129	0.000515069	0.212466	9.31717e-5	0.213525	0.00169942	0.0265114	0.00031459	0.000378467	0.253499	0.129689	0.139735	0.306227	0.108643	0.0164464	0.0074966	5.9168e-5	0.275985	0.00184365	0.32945	0.000765645	3.25797e-5	0.00355114	0.170771	0.0716283	0.0430864	0.140025	0.000312149	0.211981	0.0588576	0.0222813	0.0185314	0.215941	2.74433e-5	0.199195	0.000413384	0.00120152	0.1265	0.233413	0.000463503	0.11701	7.92475e-6	0.192853	0.0243812	0.01729	0.0490552	0.298393	0.0787513	3.91366e-6	0.065819	0.230707	0.168202	0.000139498	0.00101464	0.00924736	3.59103e-6	0.0341167	0.282092	0.016366	0.401629	0.00280864	0.00655419	0.0104324	0.0402303	0.341495	0.0353692	9.52398e-9	0.0957917	0.00105948	0.251233	0.0300049	6.99025e-9	0.00875849	0.0367777	0.000165857	0.204798	0.00480167	0.330875	2.29279e-5	0.0151905	1.80199e-6	⋯
2490	0.0054276	0.0215971	0.0381749	0.00446343	2.62429e-5	0.0002363	0.0330393	0.0253778	0.00172214	0.00129768	3.47556e-7	0.00165191	3.33388e-6	0.00702471	0.00101342	0.0323338	0.00767186	5.76446e-5	0.0033424	0.0133629	0.0357483	0.00197002	0.00218882	0.00262451	0.0250134	0.0117098	0.000407833	0.00782078	0.00153907	0.0112463	9.10016e-6	0.000403685	0.00715009	0.000451283	0.0301059	0.000128739	0.00136635	1.96349e-5	0.000222013	0.0313901	1.32759e-5	0.0269246	2.92917e-7	0.0111933	0.018759	0.00601015	0.029701	0.00355944	0.0154726	0.0192494	0.0281955	0.00235878	0.014405	0.00476844	0.0305565	0.0380793	0.000316968	0.000455479	0.0245106	0.00988176	0.00353638	0.00399585	0.0326985	0.0174238	0.0240757	0.000715605	0.0191275	0.00833476	0.00923722	0.000325124	0.00301449	0.0291526	0.00246509	0.0142601	0.00134761	0.0293476	3.28188e-5	0.000443128	1.87429e-6	0.0258475	0.0228655	0.00164805	0.000205147	6.44099e-6	0.0298188	1.7407e-8	0.0136643	1.53648e-5	0.00210354	0.010775	3.06263e-7	0.00160099	0.00333795	0.0180979	0.00305981	0.0144034	0.000270632	0.00049835	0.0102375	0.00227582	⋯
2491	0.0808199	0.000794376	0.0425222	0.00062567	3.27847e-6	0.168644	0.0555419	0.0013089	0.173004	0.0337269	5.46804e-8	0.0497232	1.31159e-6	0.166066	5.7273e-5	0.00307234	0.0566901	0.00333656	0.0541962	0.129057	0.0048201	0.0440203	0.000145567	0.0550713	0.00458339	0.014575	0.000133115	0.000728373	0.154482	0.0900673	0.00309273	0.0608157	0.160703	0.00177686	0.0139485	2.58847e-5	0.10168	9.14718e-5	0.0615061	0.00412761	5.76758e-6	0.0386838	0.00250509	0.140785	0.0312519	0.15749	0.00260934	0.0736251	0.120993	0.0931349	0.02445	0.164826	0.000317529	0.0943584	0.00305073	0.00749966	0.00726843	0.0273683	0.0139772	0.132776	4.18361e-5	0.144993	0.0602828	0.0144519	0.0856723	0.0779595	0.0739297	7.0179e-5	0.151133	0.0235316	0.0346942	0.00208828	0.00070353	0.00813487	1.45283e-5	0.0471062	0.0209619	0.00175499	0.0303658	0.0355344	0.0557512	0.0025492	0.00219747	0.02448	0.0638736	1.24435e-9	0.0546999	5.69092e-5	0.115541	0.0150531	4.1856e-9	0.135166	0.177315	0.000927776	0.0552487	0.060146	0.0692488	2.7857e-5	0.00906307	1.29699e-5	⋯
2492	0.0390656	0.000380698	0.0166329	0.00116758	1.30548e-5	0.0139768	0.0159803	0.000616106	0.0254898	0.0243649	3.77735e-7	0.0277939	6.35799e-6	0.0250408	0.000145005	0.00142342	0.0394284	0.00518056	0.0339688	0.0243063	0.00245803	0.0287344	0.00032251	0.0320525	0.00458371	0.0177603	0.000344282	0.00121736	0.0260391	0.0419341	0.00263347	0.0177769	0.0292349	0.003768	0.0118507	8.09e-5	0.0270605	0.000298664	0.0143544	0.00304334	2.31796e-5	0.00638083	0.00068393	0.0274231	0.0280231	0.0320103	0.00172528	0.0355493	0.0281541	0.0168551	0.0197616	0.028107	0.000310436	0.0379782	0.00206852	0.00421721	0.0107803	0.0168014	0.00178087	0.0358884	8.92234e-5	0.033213	0.023344	0.0160132	0.0312194	0.0219652	0.038874	7.85112e-5	0.0237344	0.0145823	0.0303189	0.00103409	0.00140026	0.010263	3.83422e-5	0.0282744	0.00618153	0.00373021	0.00216668	0.00542219	0.00857603	0.0047819	0.00469133	0.00339923	0.0284081	1.33612e-8	0.0391065	0.000195018	0.0303936	0.0185103	3.44049e-8	0.00301892	0.0103065	0.00100647	0.033343	0.00569469	0.0155351	7.83625e-5	0.0121478	3.10911e-5	⋯
2493	0.0350845	0.000102597	0.00742558	0.00132454	2.78005e-5	0.00594413	0.00594422	0.000174883	0.012038	0.0345602	1.35353e-6	0.0325487	2.01577e-5	0.00904247	0.000178504	0.000450164	0.042129	0.0166589	0.0388103	0.00889308	0.000905312	0.0362419	0.000365582	0.0361687	0.00320533	0.0215679	0.000612224	0.00115767	0.0138059	0.0322253	0.00718311	0.0164189	0.0122847	0.0092758	0.00866359	0.000153752	0.0200309	0.00110256	0.0123309	0.00153209	6.37826e-5	0.00142467	0.00133201	0.0108868	0.0276944	0.0152583	0.000729693	0.0340556	0.0118495	0.00519196	0.0156804	0.0140041	0.000139595	0.0308304	0.0009171	0.00182216	0.0277107	0.0244799	0.000264806	0.0191418	6.93808e-5	0.0187247	0.011423	0.0165254	0.0165274	0.0184252	0.0286909	3.27307e-5	0.00824837	0.0222284	0.0430748	0.000327373	0.00193473	0.0109926	3.51631e-5	0.0196301	0.00767162	0.00921175	0.00150903	0.0011222	0.00200349	0.00879254	0.0141795	0.00314648	0.016137	5.98159e-8	0.0378681	0.000717015	0.0209658	0.0231826	8.79976e-8	0.000271266	0.00192433	0.000585656	0.0376732	0.000975044	0.0126959	0.000105394	0.0149625	2.3116e-5	⋯
2494	0.00730314	0.00136479	0.00991354	0.0046281	8.36574e-5	0.000144845	0.00654948	0.00185108	0.000957513	0.00433319	2.9396e-6	0.00398391	2.58399e-5	0.00194395	0.00106651	0.00313678	0.0127586	0.000952308	0.00704086	0.00332134	0.00480983	0.00520203	0.00203189	0.00561369	0.0112355	0.0189129	0.00107986	0.00554993	0.00105498	0.0104213	0.000144637	0.000802046	0.00258408	0.0031005	0.0167887	0.000358128	0.00164003	0.000287914	0.000422968	0.00740888	8.07881e-5	0.00207048	4.05854e-6	0.00335235	0.0217764	0.00279101	0.00482087	0.00564863	0.00480451	0.00316648	0.0199908	0.00137012	0.00213347	0.00557117	0.00550378	0.00723109	0.00304845	0.00180518	0.000788708	0.0051	0.00128857	0.00263716	0.0107933	0.0197505	0.0100904	0.001139	0.0148069	0.000934826	0.00225968	0.00143698	0.00895972	0.00271853	0.00405329	0.0163676	0.000652985	0.0178812	0.00013754	0.0030645	4.43305e-6	0.00172148	0.00206227	0.00582257	0.0023662	2.0616e-5	0.0131802	1.96258e-7	0.0173826	0.000216259	0.00212552	0.0187159	9.83663e-7	4.30443e-5	0.000315129	0.0048365	0.00641721	0.000800903	0.000465636	0.000614611	0.0161382	0.000686994	⋯
2495	0.0534164	0.00343085	0.0682771	0.0010763	4.41665e-6	0.0639547	0.0878684	0.00512082	0.0977932	0.0161459	5.59085e-8	0.0250009	1.14845e-6	0.146199	0.00011686	0.00999129	0.0407555	0.000937366	0.0316833	0.138379	0.0135628	0.022936	0.000296586	0.030531	0.00913453	0.0148472	0.000155881	0.00148279	0.0831051	0.0732523	0.000624469	0.0235105	0.134085	0.00102933	0.0221041	3.28755e-5	0.0507182	4.28587e-5	0.0213675	0.0102131	5.24781e-6	0.0789626	0.00028152	0.137286	0.032559	0.120223	0.0076237	0.0439343	0.128141	0.12575	0.0324543	0.0988564	0.00123599	0.0610673	0.00848514	0.0180828	0.00277284	0.0107121	0.0418857	0.112357	0.000148193	0.0956626	0.0830164	0.0173085	0.0982723	0.0338039	0.0721581	0.000350949	0.147022	0.00866049	0.0199611	0.00723681	0.000946528	0.0104369	4.89609e-5	0.0575985	0.0051457	0.00101394	0.00459268	0.0756888	0.0990069	0.00206855	0.000929117	0.00455059	0.0797448	1.38248e-9	0.047134	2.81356e-5	0.0636358	0.0147417	8.60394e-9	0.16211	0.172512	0.00267023	0.0316723	0.108058	0.0250247	5.48545e-5	0.0099127	5.52636e-5	⋯
2496	0.0053714	0.00325088	0.0138591	0.00352113	3.98377e-5	0.000157106	0.010311	0.00417651	0.00104822	0.00220399	9.36643e-7	0.00232516	8.66337e-6	0.002779	0.000794829	0.006348	0.00849758	0.000252449	0.00430427	0.00489437	0.00854435	0.00290997	0.00159947	0.00342414	0.0121536	0.0123932	0.000564036	0.00491104	0.00105869	0.0088617	4.00791e-5	0.000520568	0.00330355	0.00113382	0.0168467	0.000181849	0.00129776	7.70406e-5	0.000282086	0.0104159	2.99977e-5	0.0048192	1.2535e-6	0.00459106	0.016445	0.00322295	0.00786768	0.0038996	0.00643276	0.0054785	0.0183373	0.00146439	0.00353801	0.00436497	0.00861954	0.0111994	0.000997002	0.000866575	0.00259805	0.00557669	0.00145545	0.00264514	0.0137677	0.0148968	0.0117782	0.000806449	0.0134113	0.00171359	0.00338029	0.000661225	0.00472131	0.00556756	0.00257514	0.0122082	0.000654192	0.017641	6.69906e-5	0.00111771	2.88358e-6	0.00425057	0.00453043	0.00276556	0.000708923	1.15555e-5	0.014927	5.47986e-8	0.012815	5.86454e-5	0.00179795	0.0119209	4.50764e-7	0.000165464	0.000713488	0.00636109	0.00393696	0.00216965	0.000323566	0.000428749	0.0106243	0.000829364	⋯
2497	0.00648017	0.000208848	0.00322754	0.00513649	0.000218105	5.02761e-5	0.00167757	0.000306047	0.000368422	0.00723865	1.56234e-5	0.00514363	0.000118441	0.000504384	0.0012988	0.00060817	0.0141369	0.00492006	0.00868173	0.000861604	0.00115534	0.00738497	0.00221916	0.00684139	0.00649343	0.0240617	0.00224631	0.00480617	0.000475082	0.00738807	0.000619067	0.000774541	0.000823083	0.0105526	0.0105014	0.000807557	0.00115565	0.0017059	0.000374688	0.00272061	0.00030603	0.000264087	1.17245e-5	0.00097692	0.0210585	0.00105852	0.00138561	0.00550828	0.00150144	0.000641076	0.0141553	0.000562862	0.00064821	0.00433359	0.00169098	0.00217468	0.0114007	0.00318081	5.68144e-5	0.00222351	0.000814213	0.00126483	0.00402142	0.0199749	0.00423989	0.000958704	0.00979154	0.00025133	0.000552415	0.00270252	0.0147562	0.00052482	0.00596657	0.0172817	0.000515405	0.0106904	0.000204345	0.0104761	3.20437e-6	0.000197878	0.000282736	0.0130784	0.0108622	2.1309e-5	0.00603306	1.38963e-6	0.016641	0.00126841	0.00135248	0.0248885	3.14678e-6	1.72853e-6	3.38924e-5	0.00212204	0.00781065	7.27477e-5	0.000384826	0.000844023	0.0207878	0.000402443	⋯
2498	0.178456	3.63052e-6	0.0037017	0.000450471	1.63895e-5	0.164879	0.00320556	8.0585e-6	0.107834	0.392148	1.54324e-6	0.335162	3.62526e-5	0.0249949	3.71587e-5	3.39065e-5	0.163715	0.702686	0.264242	0.0141773	9.88042e-5	0.332938	7.77088e-5	0.286097	0.000799752	0.0310079	0.000552718	0.000255681	0.138585	0.0951012	0.78439	0.368362	0.0382116	0.0553538	0.0040493	0.000109283	0.231658	0.0109181	0.377829	0.000224173	0.000106142	0.000420197	0.793638	0.0222167	0.0390184	0.0592132	6.9731e-5	0.227061	0.0189965	0.00451939	0.0117903	0.105253	5.93723e-6	0.16972	9.95173e-5	0.000295582	0.523071	0.483338	2.92136e-5	0.0545252	3.42782e-6	0.105281	0.00856148	0.0156678	0.0202529	0.306492	0.0514376	7.04928e-7	0.0174136	0.515948	0.293198	2.06623e-5	0.00125081	0.00854338	1.92892e-6	0.0190558	0.572065	0.0553467	0.486315	0.000298667	0.00086201	0.0218411	0.194579	0.539405	0.0154042	5.17345e-8	0.090781	0.00612133	0.1868	0.0369511	1.54487e-8	0.000422317	0.00461218	5.46626e-5	0.271458	0.000432898	0.35224	2.37814e-5	0.0176045	6.8995e-7	⋯
2499	0.00416251	0.000142583	0.00196118	0.00590337	0.000338234	1.54751e-5	0.00089929	0.000206719	0.000145248	0.00541423	3.00222e-5	0.00347225	0.00019604	0.000203805	0.00169967	0.000405568	0.0105583	0.00493104	0.00617163	0.000383658	0.000797118	0.00531611	0.00275066	0.00471055	0.00582118	0.0228418	0.00290919	0.00531237	0.000197488	0.00470737	0.000480577	0.000398156	0.000358279	0.0123064	0.00874156	0.00113265	0.00057687	0.00235359	0.000175127	0.00217086	0.000469999	0.000113914	5.6187e-6	0.000439565	0.0176619	0.00048053	0.00106378	0.00352428	0.000728255	0.000283103	0.0114766	0.000235963	0.000589157	0.00258219	0.00130936	0.00155295	0.0113339	0.00221151	2.23586e-5	0.00113738	0.00101275	0.000599766	0.00242081	0.0183453	0.00247447	0.00048743	0.00659615	0.000243379	0.00022771	0.00188217	0.0120002	0.000360365	0.0069802	0.016771	0.00069787	0.00765859	0.00010379	0.0122274	8.9926e-7	8.28168e-5	0.000117982	0.0146793	0.0125654	7.89382e-6	0.00382909	3.20365e-6	0.0126468	0.00178801	0.000674228	0.0236819	6.64705e-6	3.24864e-7	9.3284e-6	0.00199124	0.00546402	2.49331e-5	0.000177239	0.00115142	0.020563	0.00052827	⋯
2500	0.00524883	0.000603337	0.00522672	0.00437745	0.000112842	6.08632e-5	0.00307936	0.000833593	0.000453228	0.00395031	5.26329e-6	0.0032488	4.2004e-5	0.000828433	0.00108588	0.00147232	0.0103244	0.00133715	0.00577254	0.00146583	0.00243147	0.00446239	0.00196582	0.00453436	0.00804585	0.0173681	0.00127866	0.0048178	0.000532679	0.00699896	0.000175398	0.000545082	0.00120028	0.00407559	0.0119758	0.000446732	0.00100682	0.00047201	0.00027091	0.00446228	0.000120973	0.000735784	3.69127e-6	0.00153587	0.0176661	0.0013843	0.00266768	0.00416607	0.00230777	0.00129697	0.0146404	0.000675736	0.00127502	0.00375146	0.00311442	0.00394957	0.00388731	0.00162776	0.00022988	0.00274942	0.00106564	0.00142201	0.00589685	0.016772	0.00564434	0.00074099	0.00994422	0.000551748	0.000953391	0.00132119	0.0083839	0.00128794	0.00418655	0.014356	0.000593788	0.0119258	0.00010541	0.00403462	2.22714e-6	0.000587208	0.000739423	0.0066944	0.00338141	1.26956e-5	0.00778956	4.07699e-7	0.0135692	0.000356252	0.00126635	0.0174232	1.57372e-6	8.11584e-6	9.06461e-5	0.00325151	0.00521126	0.000237469	0.000290768	0.000657678	0.0150958	0.000566752	⋯

Problem 1(e)

Let $\hat{s}_{jt}$ denote the predicted market share of product $j$ at time $t$, which is given by:

$\hat{s}_{jt} (\delta_t, \theta_2) = \int \frac{\exp(\delta_{jt} + \mu_{ijt}(\theta_2))}{1 + \sum_k \exp(\delta_{kt} + \mu_{ikt}(\theta_2))} dF_t(i)$

where $\theta_2$ denotes the non-linear parameters of this model, $\delta_t$ is the vector of mean utilities of all products, and $F_t(i)$ denotes the distribution of consumer characteristics in market $t$.

Given that we have a sample of 500 consumers drawn from $F_t(i)$, we can approximate the integral over $F_t(i)$ by summing over all consumers in market $t$:

$\hat{s}_{jt} (\delta_t, \theta_2) \approx \sum_{i \in \mathcal{S}_t} s_{ijt} = \sum_{i \in \mathcal{S}_t} \frac{\exp(\delta_{jt} + \mu_{ijt}(\theta_2))}{1 + \sum_{k \in \mathcal{J}_t} \exp(\delta_{kt} + \mu_{ikt}(\theta_2))}$

Where $\mathcal{S}_t$ denotes the set of simulated consumers.

Write a function that aggregates the individual choice probabilities of all consumers to compute the predicted market shares for each product, as a function of $\delta_t = (\delta_{1t}, \ldots, \delta_{J_t,t})$ and $\theta_2 = (\alpha, \sigma_1, \sigma_2)$. Given nonlinear parameters $\alpha^y = 0$, $\sigma_1^y = 1$, $\sigma_2^y = 1$, the vector $\delta_t^0$ from Problem 1(d), and the sample of simulated consumers you obtained in part 1 (b), compute $\hat{s}_{jt}(\delta_t^0, \theta_2^y)$ for all products and all markets. Compute $\bar{p}^{''' } = \sum_j p_{jt} \hat{s}_{jt}(\delta_t^0, \theta_2^y)$ and $\text{cor}(\bar{p}^{'''}, \bar{y}_t)$. Compare the correlation between the average prices paid by consumers predicted by the model and those observed in the data (from Problem 1(a)). What does the difference between these results suggest about the true value of $\alpha$?

function s_jt_hat(simulated_consumers::DataFrame, delta_table::DataFrame, product_data::DataFrame, nonlinear_param::Dict{Symbol, Float64})
    
    alpha = nonlinear_param[:alpha]
    sigma_1 = nonlinear_param[:sigma_1]
    sigma_2 = nonlinear_param[:sigma_2]
    
    product_variables = innerjoin(delta_table, product_data, on = [:product_ID, :t])
    
    grouped_simulated = groupby(simulated_consumers, :t)
    grouped_product = groupby(product_variables, :t)

    product_shares_all = DataFrame()
    
    for tt in 1:length(grouped_product)
        
        data = sort(grouped_product[tt], :product_ID)
        n_j = nrow(data)
        
        consumer_census_t = sort(grouped_simulated[tt], :consumer_ID)
        
        income_data = transpose(repeat(consumer_census_t.income, 1, n_j))
        sugar_taste = transpose(repeat(consumer_census_t.nu_1, 1, n_j))
        caffeine_taste = transpose(repeat(consumer_census_t.nu_2, 1, n_j))
        
        mu_ijt = alpha .* data.price .* income_data .+ sigma_1 .* sugar_taste .* data.sugar .+ sigma_2 .* caffeine_taste .* data.caffeine
        utility = mu_ijt .+ data.delta
        
        # Calculate utility from outside option
        outside_option = zeros(size(utility, 2))     
       
         # Calculate inclusive values (log-sum-exp over rows)
        inclusive_values = log.(sum(exp.(vcat(utility, outside_option')), dims=1))
        
        current_shares = exp.(utility .- inclusive_values)
        
        product_shares = sum(current_shares, dims=2) ./ nrow(consumer_census_t)
        
        current_shares_df = DataFrame(product_shares, [:s_jt_hat])

        current_shares_df.product_ID = data.product_ID
        current_shares_df.t = fill(tt, nrow(current_shares_df))
        
        append!(product_shares_all, current_shares_df)
        
    end
    
    return product_shares_all
    
end

s_jt_hat (generic function with 1 method)

predicted_shares = s_jt_hat(sample, delta_table, product_data, theta_2_double_prime)

product_comb = innerjoin(product_data, predicted_shares, on=[:product_ID, :t])

# Add market_share column to product_comb DataFrame
product_comb[!, :market_share] = product_comb[!, :s_jt_hat]

grouped_product = groupby(product_comb, :t)
p_bar = combine(grouped_product) do sdf
    DataFrame(p_bar = mean(sdf.price .* sdf.market_share))
end

market_data = innerjoin(mean_income, p_bar, on = :t)

correlation = cor(market_data.mean_income, market_data.p_bar)
println(correlation)

0.044607516356691856

Problem 2(a)

Let $s_{jt}$ be the market shares for product $j$ in market $t$ observed in the data, and $\hat{s}_{jt} (\delta_t, \theta_2)$ be the market shares predicted by the model. Berry (1994) shows that, given $\theta_2$, there exists a unique vector $\hat{\delta}_t (\theta_2)$ that equates the observed and the predicted market shares:

$$s_{jt} = \hat{s}_{jt} \left( \hat{\delta}_t (\theta_2) , \theta_2 \right)$$

In order to find the vector $\hat{\delta}_t (\theta_2)$ that equates the observed and the predicted market shares, we can use the following contraction mapping:

$$\delta_{jt}^{(\tau+1)} = \delta_{jt}^{(\tau)} + \log(s_{jt}) - \log\left( \hat{s}_{jt} \left( \delta_{jt}^{(\tau)}, \theta_2 \right) \right)$$

where $\delta_{jt}^{(\tau)}$ denotes the current-period guess for $\hat{\delta}_{jt} (\theta_2)$. Write a function that employs this contraction mapping (or one of the alternatives discussed in Conlon & Gortmaker (2020)) to find $\hat{\delta}_t (\theta_2)$ for any given $\theta_2$.

function delta_inversion(delta_initial, next_delta, tolerance)
    delta_t = delta_initial

    while norm((delta_t) .- next_delta(delta_t)) > tolerance
        delta_t = next_delta(delta_t)
    end

    return delta_t
end

function delta_hat(nonlinear_param, product_data, simulated_consumers, tolerance)
    
    each_period_product_data = groupby(product_data, :t)
    each_period_census = groupby(simulated_consumers, :t)
    
    T = length(each_period_product_data)
    
    final_delta = DataFrame()
    
    function processPeriod(tt)
        current_period_var = sort!(each_period_product_data[tt], :product_ID)
        consumer_census_tt = sort!(each_period_census[tt], :consumer_ID)
        
        current_period_var[!, :delta] = log.(current_period_var[!, :market_share]) .- log(1 - sum(current_period_var[!, :market_share]))
        
        delta_initial = select(current_period_var, :product_ID, :delta, :t)
        
        select!(current_period_var, Not(:delta))
        
        s_jt_obs = current_period_var[!, :market_share]

        current_period_var_df = DataFrame(current_period_var)
        consumer_census_tt_df = DataFrame(consumer_census_tt)

function next_delta(par)
            delta_current = DataFrame(product_ID = delta_initial[!, :product_ID], delta = par, t = delta_initial[!, :t])
            
            # Compute predicted shares using current delta
            s_jt_mod = s_jt_hat(consumer_census_tt_df, delta_current, current_period_var_df, nonlinear_param)
            sort!(s_jt_mod, :product_ID)
            s_jt_pred = s_jt_mod[!, :s_jt_hat]
            
            # Update delta
            delta_new = delta_current[!, :delta] .+ log.(s_jt_obs) .- log.(s_jt_pred)
            
            return delta_new
        end

        result = delta_inversion(delta_initial[!, :delta], next_delta, tolerance) 
        
        delta_current = DataFrame(product_ID = delta_initial[!, :product_ID], delta = result, t = delta_initial[!, :t])
        
        return delta_current
    end
    
    for tt in 1:T
        delta_result = processPeriod(tt)
        append!(final_delta, delta_result)
    end
    
    return final_delta
end

delta_hat (generic function with 1 method)

Problem 2(b)

In this exercise, we will study how the accuracy of $\hat{\delta}_t (\theta_2)$ changes depending on the size of our sample of simulated consumers. As before, let $\theta_2^y = (0, 1, 1)$. For market $t = 1$ and $n = 1, 2, \ldots, 100$, draw a sample $\mathcal{S}_n$ of consumers using the function you wrote for Problem 1(a), each of them containing $S_n = 100$ simulated consumers (100 different samples, each sample must include 100 consumers). For each sample $\mathcal{S}_n \in \{\mathcal{S}_1, \mathcal{S}_2, \ldots, \mathcal{S}_{100}\}$, let

$$\hat{\delta}_n(\theta_2^y) = \left( \hat{\delta}_{1t_n}(\theta_2^y), \ldots, \hat{\delta}_{J_t,t_n}(\theta_2^y) \right)$$

1. Plot the distribution of $\hat{\delta}_{1t_n}(\theta_2^y)$.

2. Report the sample mean:

$$\bar{\delta}_{1t} = \frac{1}{100} \sum_{n=1}^{100} \hat{\delta}_{1t_n}(\theta_2^y)$$

3. Report the sample variance

$$\frac{1}{100} \sum_{n=1}^{100} \left( \hat{\delta}_{1t_n}(\theta_2^y) - \bar{\delta}_{1t} \right)^2$$

4. Repeat 1, 2, and 3 increasing the size of each sample from 100 to 500. How do the results from parts 2(a) and 2(b) compare to each other?

function delta_hat_sim(nonlinear_param, product_data, n_sim, n_salmple)
    
    first_market_data = product_data[product_data.t .== 1, :]
    
    function compute_delta(i)
        
        sample_i = simulate_consumers(consumer_census, n_salmple, i)
        
        first_market_consumers = sample_i[sample_i.t .== 1, :]
        
        delta_hat_i = delta_hat(theta_2_double_prime, first_market_data, first_market_consumers, 1e-8)
        
        delta_i_1 = delta_hat_i[!, :delta][1]
        
        return(delta_i_1)
        
    end
    
    results = zeros(n_sim)  # create an empty array that can hold any type
    
    for j in 1:n_sim
        
        results[j] = compute_delta(j)
        
    end
    
    return(results)
    
end

delta_hat_sim (generic function with 1 method)

delta_tilde_1 = delta_hat_sim(theta_2_double_prime, product_data, 100, 100)

mean(delta_tilde_1)

3.005090488551078

var(delta_tilde_1)

0.8992948086119623

histogram(delta_tilde_1, bins=30, label="delta_tilde_1", alpha=0.6, color=:blue)

delta_tilde_2 = delta_hat_sim(theta_2_double_prime, product_data, 100, 500)

mean(delta_tilde_2)

2.872556452151993

var(delta_tilde_2)

0.14957281962634364

histogram(delta_tilde_2, bins=30, label="delta_tilde_2", alpha=0.6, color=:blue)

Problem 3(a)

In this problem, we will follow Gandhi and Houde (2023) to construct differentiation instruments $A_j (x_t, w_t)$ using the products' observed characteristics

$$x_{jt} = (1, sug_{jt}, caf_{jt})$$

and price instruments

$$w_{jt} = (\text{corn\_syrup\_price}_t \times sug_{jt}, \text{caffeine\_extract\_price}_t \times caf_{jt})$$

Construct the instruments $A_j (x_t, w_t)$

function build_GH_instruments(product_data, mean_income)
    # Copy product_data
    dd = copy(product_data)
    
    # Fit the model and predict
    model = lm(@formula(price ~ sugar * corn_syrup_price + caffeine * caffeine_extract_price), dd)
    dd.p_hat = predict(model)

    # Creating pairwise combinations
    dd_pairs = innerjoin(dd, dd, on = :t, makeunique=true, validate=(false, false))
    dd_pairs = filter(row -> row.product_ID ≠ row.product_ID_1, dd_pairs)

    # Calculate squared differences for sugar, caffeine, and p_hat
    dd_pairs[!, :sq_diff_sugar] = (dd_pairs[!, :sugar] .- dd_pairs[!, :sugar_1]).^2
    dd_pairs[!, :sq_diff_caffeine] = (dd_pairs[!, :caffeine] .- dd_pairs[!, :caffeine_1]).^2
    dd_pairs[!, :sq_diff_p_hat] = (dd_pairs[!, :p_hat] .- dd_pairs[!, :p_hat_1]).^2

    # Aggregate these differences by product_ID and t
    aggregated_diff_sugar = combine(groupby(dd_pairs, [:product_ID, :t]), 
                                    :sq_diff_sugar => sum => :GH_instrument_sugar )
    aggregated_diff_caffeine = combine(groupby(dd_pairs, [:product_ID, :t]), 
                                       :sq_diff_caffeine => sum => :GH_instrument_caffeine )
    aggregated_diff_p_hat = combine(groupby(dd_pairs, [:product_ID, :t]), 
                                    :sq_diff_p_hat => sum => :GH_instrument_p_hat )

    # Join aggregated results back to the original DataFrame
    dd = leftjoin(dd, aggregated_diff_sugar, on = [:product_ID, :t])
    dd = leftjoin(dd, aggregated_diff_caffeine, on = [:product_ID, :t])
    dd = leftjoin(dd, aggregated_diff_p_hat, on = [:product_ID, :t])

    # Merging with mean_income and other computations
    instrument_data = innerjoin(dd, mean_income, on=:t)
    instrument_data.avg_income_inst = instrument_data.mean_income .* instrument_data.p_hat

    # Adding a constant column
    instrument_data.constant = ones(nrow(instrument_data))

    return instrument_data
end

build_GH_instruments (generic function with 1 method)

instrument_data = build_GH_instruments(product_data, mean_income)

2500×16 DataFrame

2475 rows omitted

Row	product_ID	price	sugar	caffeine	corn_syrup_price	caffeine_extract_price	quantity	t	market_share	p_hat	GH_instrument_sugar	GH_instrument_caffeine	GH_instrument_p_hat	mean_income	avg_income_inst	constant
	Int64	Float64	Float64	Float64	Float64	Float64	Int64	Int64	Float64	Float64	Float64?	Float64?	Float64?	Float64	Float64	Float64
1	1	2.02133	1.58086	6.07446	0.199022	0.264899	223	1	0.0223	1.92786	298.941	91.7637	7.51015	6.00259	11.5721	1.0
2	1	1.94121	2.40883	5.37299	0.189651	0.263386	106	2	0.0106	1.87266	183.9	53.1036	7.74265	4.5502	8.52098	1.0
3	1	1.65012	1.37439	4.79142	0.190519	0.283622	81	3	0.0081	1.62128	342.576	35.6997	18.6042	5.9295	9.61336	1.0
4	1	1.65104	1.829	5.40206	0.196738	0.234751	78	4	0.0078	1.63139	317.574	63.7647	18.9123	4.43252	7.23115	1.0
5	1	2.18907	1.15713	6.59757	0.31464	0.250255	712	5	0.0712	2.02285	381.115	157.937	16.5204	4.45157	9.00485	1.0
6	1	1.16253	1.94425	4.57279	0.174609	0.182037	43	6	0.0043	1.17332	256.924	58.7992	9.35902	5.5722	6.53798	1.0
7	1	1.70539	1.75074	7.43689	0.290227	0.161537	634	7	0.0634	1.72901	276.034	357.133	12.8515	4.85486	8.39411	1.0
8	1	1.68218	1.1832	4.27315	0.204226	0.314919	131	8	0.0131	1.58654	443.935	62.6004	35.5052	4.51054	7.15617	1.0
9	1	1.72647	2.36464	5.2934	0.26121	0.242201	8	9	0.0008	1.90328	211.132	72.4287	13.703	4.78392	9.10513	1.0
10	1	1.99782	1.55731	5.14066	0.220158	0.299709	243	10	0.0243	1.88404	302.391	38.3452	16.645	5.1851	9.76894	1.0
11	1	1.74007	2.00323	4.53769	0.205343	0.290355	48	11	0.0048	1.72775	253.601	57.7183	19.9383	4.8054	8.30254	1.0
12	1	2.06338	2.30017	5.29768	0.239624	0.282229	49	12	0.0049	2.04745	176.08	47.2596	12.0884	5.96622	12.2155	1.0
13	1	1.40313	1.6347	4.23523	0.348398	0.213371	15	13	0.0015	1.47531	284.63	112.836	51.5643	4.69344	6.92428	1.0
⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮	⋮
2489	50	3.13989	5.72339	6.28477	0.236956	0.278179	306	39	0.0306	3.10466	410.4	87.5781	40.5707	5.93635	18.4304	1.0
2490	50	1.92962	5.18098	4.0355	0.116454	0.308221	235	40	0.0235	1.83187	298.785	92.9716	5.94231	5.84791	10.7126	1.0
2491	50	2.09891	3.80218	4.8143	0.240711	0.219965	86	41	0.0086	1.97515	176.074	52.4533	10.6179	5.42932	10.7237	1.0
2492	50	3.52263	6.44061	6.66856	0.31482	0.247241	291	42	0.0291	3.68829	584.269	196.069	112.643	4.77021	17.5939	1.0
2493	50	3.47946	8.20551	3.96725	0.288353	0.202105	3913	43	0.3913	3.16926	1363.41	90.9647	88.1933	5.74445	18.2056	1.0
2494	50	2.36507	4.21225	4.26267	0.274055	0.279371	52	44	0.0052	2.34486	151.358	78.0126	13.597	5.54058	12.9919	1.0
2495	50	2.35623	3.57038	5.4192	0.321409	0.222686	29	45	0.0029	2.36198	145.944	69.8106	18.8515	4.28395	10.1186	1.0
2496	50	2.1976	4.29359	4.13298	0.24389	0.265363	76	46	0.0076	2.14135	162.411	57.5581	11.9051	5.65435	12.1079	1.0
2497	50	2.22313	4.76312	5.1124	0.179816	0.284419	22	47	0.0022	2.30425	229.073	55.9681	11.5902	5.55554	12.8014	1.0
2498	50	2.86181	4.96189	6.09999	0.307171	0.227802	66	48	0.0066	2.92356	235.3	104.618	40.5804	5.56763	16.2773	1.0
2499	50	2.29245	4.22963	4.78552	0.251569	0.24075	110	49	0.011	2.21666	166.056	45.6461	10.154	5.02191	11.1319	1.0
2500	50	2.01727	3.82539	4.59095	0.210564	0.255062	63	50	0.0063	1.97377	171.148	77.847	12.7467	5.94398	11.7321	1.0

Problem 3(c)

Code the GMM objective function that you would use to estimate the nonlinear parameters $\theta_2$, using the instruments you derived in part (a). Evaluate the objective function at $\tilde{\theta} = (-0.5, 2, 2)$.

function GMM_objective(par,product_data,sample,mean_income)
    
    theta = Dict(
    :alpha => par[1],
    :sigma_1 => par[2],
    :sigma_2 => par[3]
)
    
    N = nrow(product_data)
    
    delta_table = delta_table = delta_hat(theta, product_data, sample, 1e-6)
    
    demand_instruments = build_GH_instruments(product_data, mean_income)
    
    sort!(product_data, [:t, :product_ID])
    sort!(delta_table, [:t, :product_ID])
    sort!(demand_instruments, [:t, :product_ID])
    
    X = hcat(product_data[!, :price], product_data[!, :sugar], product_data[!, :caffeine], ones(size(product_data, 1)))
    
    Z  = hcat(
    demand_instruments[!, :sugar], demand_instruments[!, :caffeine],demand_instruments[!, :p_hat],
    demand_instruments[!, :GH_instrument_sugar],demand_instruments[!, :GH_instrument_caffeine],
    demand_instruments[!, :avg_income_inst],demand_instruments[!, :constant]
)

    #fix type just in case
    Z = Matrix{Float64}(Z)
    
    delta = hcat(delta_table[!, :delta])

    W = inv((Z' * Z) / N)
    
    P_Z = Z * inv(Z' * Z) * Z'
    
    Beta_2sls = inv(X' * P_Z * X) * X' * P_Z * delta
    
    xi_hat = delta - X * Beta_2sls
    
    g = ( Z' * xi_hat ) / N
    
    objective_value = g' * W * g
    
    return objective_value[1]
    
end
    

GMM_objective (generic function with 1 method)

GMM_objective([-0.5,2,2],product_data,sample,mean_income)

1.1737387832864838

Problem 3(d)

Write a function that computes the analytic gradient as a function of the vector of nonlinear parameters $\theta_2$.

The gradient of the GMM objective function is

$$\nabla q (\theta_2) = 2G(\theta_2)^T W g (\theta_2)$$

with

$$G (\theta_2) = \frac{1}{N} \underbrace{Z^T}_{T \times 3} \underbrace{\frac{\partial \xi}{\partial \theta_2}}_{T \times N \times 3}$$

$$\frac{\partial \xi}{\partial \theta_2} = \frac{\partial \delta}{\partial \theta_2}$$

$\frac{\partial \delta}{\partial \theta_2}$ can be separated by market $t$, after using the implicit function theorem:

$$\frac{\partial \delta_t}{\partial \theta_2} (\theta_2) = - \left( \frac{\partial \hat{s}_t}{\partial \delta_t} (\theta_2) \right)^{-1} \left[ \frac{\partial \hat{s}_t}{\partial \theta_2} (\theta_2) \right]$$

$$\underbrace{}_{J_t \times 3} \qquad \underbrace{}_{J_t \times J_t} \qquad \underbrace{}_{J_t \times 3}$$

$$\frac{\partial \hat{s}_t}{\partial \delta_t} (\theta_2) = \begin{bmatrix} \frac{\partial s_{1t}}{\partial \delta_{1t}} & \cdots & \frac{\partial s_{1t}}{\partial \delta_{J_t,t}} \\ \vdots & \ddots & \vdots \\ \frac{\partial s_{J_t,t}}{\partial \delta_{1t}} & \cdots & \frac{\partial s_{J_t,t}}{\partial \delta_{J_t,t}} \end{bmatrix}$$

With

$$\frac{\partial \hat{s}_{jt}}{\partial \delta_{kt}} = \int \frac{\partial s_{ijt}}{\partial \delta_{kt}} dF(i)$$

$$= \int s_{ijt} (1 - s_{ijt}) dF(i) \text{ if } j = k$$

$$= \int -s_{ijt} s_{ikt} dF(i) \text{ if } j \neq k$$

function partial_s_delta(product_data::DataFrame, simulated_consumers::DataFrame, theta_2::Vector, time_period::Int)
    # Assuming theta_2 is a vector with at least 3 elements
    nonlinear_param = Dict(
        :alpha   => theta_2[1],
        :sigma_1 => theta_2[2],
        :sigma_2 => theta_2[3]
    )
    
    # Filter product_data for the given time_period
    product_data_t = product_data[product_data.t .== time_period, :]

    delta_table = delta_hat(nonlinear_param, product_data_t, simulated_consumers, 1e-7)

    # You need to define or import the s_ijt function in Julia
    individual_shares_t = s_ijt(simulated_consumers, delta_table, product_data_t, nonlinear_param)

    # Assuming `product_ID` is a column in `individual_shares_t`
    individual_shares_t = individual_shares_t[sort(individual_shares_t.product_ID), :]

    deriv_matrix_t = zeros(size(individual_shares_t, 1), size(individual_shares_t, 1))

    diag_t = zeros(50)
    
    for i in 1:size(individual_shares_t, 1) - 2
        # Assuming individual_shares_t[i, :] gets the i-th column representing a consumer
        consumer_i = DataFrame(product_ID = individual_shares_t.product_ID, market_share = individual_shares_t[:, i])
        consumer_i.market_share = convert.(Float64, consumer_i.market_share)
        deriv_matrix_t_i = -1 * consumer_i.market_share * transpose(consumer_i.market_share)
        deriv_matrix_t += deriv_matrix_t_i

        diag_i = consumer_i.market_share .* (1 .- consumer_i.market_share)
        diag_t += diag_i
    end

    for i in 1:size(deriv_matrix_t, 1)
        deriv_matrix_t[i, i] = diag_t[i]
    end

    return deriv_matrix_t / 500
end

partial_s_delta (generic function with 1 method)

$$\frac{\partial \hat{s}_t}{\partial \theta_2} (\theta_2) = \begin{bmatrix} \frac{\partial s_{1t}}{\partial \alpha} & \frac{\partial s_{1t}}{\partial \sigma_1} & \frac{\partial s_{1t}}{\partial \sigma_2} \\ \vdots & \vdots & \vdots \\ \frac{\partial s_{J_t,t}}{\partial \alpha} & \frac{\partial s_{J_t,t}}{\partial \sigma_2} & \frac{\partial s_{J_t,t}}{\partial \sigma_2} \end{bmatrix}$$

$$\frac{\partial \hat{s}_{jt}}{\partial \alpha^y} = \int y_{it} p_{jt} s_{ijt} (1 - s_{ijt}) - y_{it} \left( \sum_{k \neq j} p_{kt} s_{ikt} \right) dF(i)$$

$$\frac{\partial \hat{s}_{jt}}{\partial \sigma_1} = \int \nu_{1it} sug_{jt} s_{ijt} (1 - s_{ijt}) - \nu_{1it} \left( \sum_{k \neq j} sug_{kt} s_{ikt} \right) dF(i)$$

$$\frac{\partial \hat{s}_{jt}}{\partial \sigma_2} = \int \nu_{2it} caf_{jt} s_{ijt} (1 - s_{ijt}) - \nu_{2it} \left( \sum_{k \neq j} caf_{kt} s_{ikt} \right) dF(i)$$

function partial_s_theta(product_data::DataFrame, simulated_consumers::DataFrame, theta_2::Vector, time_period::Int)
    nonlinear_param = Dict(
        :alpha   => theta_2[1],
        :sigma_1 => theta_2[2],
        :sigma_2 => theta_2[3]
    )

    product_data_t = product_data[product_data.t .== time_period, :]

    delta_table = delta_hat(nonlinear_param, product_data_t, simulated_consumers, 1e-7)

    individual_shares_t = s_ijt(simulated_consumers, delta_table, product_data_t, nonlinear_param)
    individual_shares_t = individual_shares_t[sort(individual_shares_t.product_ID), :]

    deriv_matrix_theta_t = zeros(size(individual_shares_t, 1), 3)

    for i in 1:size(individual_shares_t, 1) - 2
        consumer_i = DataFrame(
            product_ID = individual_shares_t.product_ID, 
            market_share = individual_shares_t[:, i],
            income = simulated_consumers.income[i],
            nu_1 = simulated_consumers.nu_1[i],
            nu_2 = simulated_consumers.nu_2[i]
        )
        
        select!(product_data_t, :product_ID, :price, :sugar, :caffeine)

        consumer_i = innerjoin(consumer_i, product_data_t, on = :product_ID)
        
        consumer_i.market_share = convert.(Float64, consumer_i.market_share)

        consumer_i.own_partial_alpha = consumer_i.income .* consumer_i.price .* consumer_i.market_share .* (1 .- consumer_i.market_share)
        consumer_i.cross_partial_alpha = consumer_i.income .* (sum(consumer_i.price .* consumer_i.market_share) .- consumer_i.price .* consumer_i.market_share)
        consumer_i.partial_alpha = consumer_i.own_partial_alpha .- consumer_i.cross_partial_alpha .* consumer_i.market_share

        consumer_i.own_partial_sigma_1 = consumer_i.nu_1 .* consumer_i.sugar .* consumer_i.market_share .* (1 .- consumer_i.market_share)
        consumer_i.cross_partial_sigma_1 = consumer_i.nu_1 .* (sum(consumer_i.sugar .* consumer_i.market_share) .- consumer_i.sugar .* consumer_i.market_share)
        consumer_i.partial_sigma_1 = consumer_i.own_partial_sigma_1 .- consumer_i.cross_partial_sigma_1 .* consumer_i.market_share

        consumer_i.own_partial_sigma_2 = consumer_i.nu_2 .* consumer_i.caffeine .* consumer_i.market_share .* (1 .- consumer_i.market_share)
        consumer_i.cross_partial_sigma_2 = consumer_i.nu_2 .* (sum(consumer_i.caffeine .* consumer_i.market_share) .- consumer_i.caffeine .* consumer_i.market_share)
        consumer_i.partial_sigma_2 = consumer_i.own_partial_sigma_2 .- consumer_i.cross_partial_sigma_2 .* consumer_i.market_share

        deriv_matrix_theta_t[:, 1] .+= consumer_i.partial_alpha
        deriv_matrix_theta_t[:, 2] .+= consumer_i.partial_sigma_1
        deriv_matrix_theta_t[:, 3] .+= consumer_i.partial_sigma_2
    end

    return deriv_matrix_theta_t ./ 500
end


function evaluated_moments(par)
    
    theta = Dict(
    :alpha => par[1],
    :sigma_1 => par[2],
    :sigma_2 => par[3]
)
    
    N = nrow(product_data)
    
    delta_table = delta_table = delta_hat(theta, product_data, sample, 1e-6)
    
    demand_instruments = build_GH_instruments(product_data, mean_income)
    
    sort!(product_data, [:t, :product_ID])
    sort!(delta_table, [:t, :product_ID])
    sort!(demand_instruments, [:t, :product_ID])
    
    X = hcat(product_data[!, :price], product_data[!, :sugar], product_data[!, :caffeine], ones(size(product_data, 1)))
    
    Z  = hcat(
    demand_instruments[!, :sugar], demand_instruments[!, :caffeine],demand_instruments[!, :p_hat],
    demand_instruments[!, :GH_instrument_sugar],demand_instruments[!, :GH_instrument_caffeine],
    demand_instruments[!, :avg_income_inst],demand_instruments[!, :constant]
)
    
    delta = hcat(delta_table[!, :delta])

    Z = Matrix{Float64}(Z)
    W = inv((Z' * Z) / N)
    
    P_Z = Z * inv(Z' * Z) * Z'
    
    Beta_2sls = inv(X' * P_Z * X) * X' * P_Z * delta
    
    xi_hat = delta - X * Beta_2sls
    
    g = ( Z' * xi_hat ) / N
    
    return Dict(:g => g, :W => W)
    
end

evaluated_moments (generic function with 1 method)

function Analytic_gradient(product_data::DataFrame, simulated_consumers::DataFrame, theta::Vector)
  
  # Split data on times (assuming product_data has a column `t` for time periods)
  each_period_data = groupby(product_data, :t)
  
  partial_delta_theta = []

  TT = length(each_period_data)

  # Process each time period
  for tt in 1:TT
    current_period_data = each_period_data[tt]
    current_period_data = DataFrame(current_period_data)

    # Assuming order function sorts the DataFrame and product_ID is a column in current_period_data
    sort!(current_period_data, :product_ID)
     
    partial_s_delta_t = partial_s_delta(current_period_data, simulated_consumers, theta, tt)
     
    partial_s_theta_t = partial_s_theta(current_period_data, simulated_consumers, theta, tt)
     
    partial_delta_theta_t = (- inv(partial_s_delta_t) * partial_s_theta_t)

    push!(partial_delta_theta, partial_delta_theta_t)
  end

  # Concatenate all results vertically
  partial_delta_theta = vcat(partial_delta_theta...)
  
  instruments = build_GH_instruments(product_data, mean_income)

  Z_D = Matrix(instruments[:, [:sugar, :caffeine, :p_hat, :avg_income_inst, :GH_instrument_sugar, :GH_instrument_caffeine, :constant]])
  
  G = (transpose(Z_D) * partial_delta_theta) / nrow(product_data)
  
  moments = evaluated_moments(theta) 
  
  g = moments[:g]
  W = moments[:W]
  
  gradient = 2 * transpose(G) * W * g
  
  return gradient
end


Gradient = Analytic_gradient(product_data,sample, [-0.5,2,2])

3×1 Matrix{Float64}:
 -1913.6682930428492
   373.8989663611442
   417.549439473164

Problem 3(e)

Estimate the vector $\hat{\theta}_2$ such that

$$\hat{\theta}_2 = \arg \min_{\theta_2} q (\theta_2) = g_1 (\theta_2)^T W g_1 (\theta_2)$$

function SolveGMM(product_data, sample, initial_guess)
    
    GMM_problem = function(theta_hat)
        objective_value = GMM_objective(theta_hat,product_data,sample,mean_income)
        return objective_value
    end

    # Optimize using BFGS
    result = optimize(GMM_problem, initial_guess, LBFGS())
    
    return result
end

SolveGMM (generic function with 1 method)

result = SolveGMM(product_data, sample, [0.00,0.5,0.5])

result.minimizer

3-element Vector{Float64}:
 0.14069303403169534
 1.0535769259103636
 1.023221164523225