Thank you all for your careful attention to my query.

I think my “minimal” working example was perhaps too minimal still. I will overshoot in the other direction and provide below a working copy of my actual code at the moment.

The key function for my query is `sim_data_like`

. It is the one that is called `n_sim`

times per call of the likelihood function. As of now, it has “too much” memory allocated within it so that each call of the likelihood function demands more than `n_sim`

times that amount of memory. (The function `CompParams`

currently sets the amount of simulations. Eventually I will want to simulate at least `n_sim=1000`

times for `n_firms=5000`

– that will cause memory issues for sure.

I’ve put an `@time`

call both at the bottom on the line that calls the likelihood function and also one one line within `sim_data_like`

. The purpose of the latter is so that you can easily see that it is the allocation of the vector of svectors that is the offending line. The “math” parts of the code are currently lean - they, for the most part, do not ask for much memory. It is just the ~10 instance of allocating vectors of svectors within `sim_data_like`

that adds up to the memory demands of that entire function. If I could reduce those only, I would trim much of the issue perhaps.

Some other notes: I know some other functions (e.g. `SimRawData`

) are also flabby at the moment. Don’t worry about those unless they have a direct bearing on the memory demand of `sim_data_like`

.

```
#setting up the linear program
# Required Packages
using Distributions, Random, LinearAlgebra, BenchmarkTools, SparseArrays
using StaticArrays, GLPK, Cbc, Ipopt, Optim, StatsBase, KernelDensitySJ
using Assignment, JuMP
function CompParams(;n_firms=50, n_sim=100, trim_percent=5, hmethod=1)
p=(
n_firms = n_firms,
# number of simulations for each iteration, i.e. n_sim 500 by 500 markets
n_sim=n_sim,
trim_percent = trim_percent, #which percentil as cutoff for trim fn
hmethod = hmethod) #hmethod=1 --> Silverman's RoT. =2 --> Sheather and Jones (1991) bandwidth
return p
end
function StrucParams(b::Vector{Float64})
Σ_m = exp(b[1])
Σ_f = exp(b[2])
Σ_u = [Σ_m 0.0;
0.0 Σ_f]
#Set up thetas
theta_0 = [b[3] b[4];
b[5] b[6]]
theta_1 = [b[7:8]';
b[9:10]']
theta_2 = [b[11:12]';
b[13:14]']
#set up Cov matrix for shocks to measures
sigma_1 = exp(b[15])
sigma_2 = exp(b[17])
cov12 = -sqrt(sigma_1)sqrt(sigma_2) + 2.0*sqrt(sigma_1)sqrt(sigma_2)*(exp(b[16]))/(1.0+exp(b[16])) #ensuring pos semi def cov matrix
Σ = [sigma_1 cov12;
cov12 sigma_2]
mean_price = b[18]
rm = b[19:21]
rf = b[22:24]
risk = exp(b[25]) #exp to keep risk coefficient positive
p = (
Σ_m = Σ_m,
Σ_f = Σ_f,
Σ_u = Σ_u,
theta_0 = theta_0,
theta_1 = theta_1,
theta_2 = theta_2,
Σ = Σ,
mean_price = mean_price,
rm = rm,
rf = rf,
risk = risk) #exp to keep risk coefficient positive
return p
end
function DrawUnobsData(n,Σ)
svectorscopy(rand(MvNormal([0.0, 0.0],Σ),n), Val{2}())
end
function svectorscopy(x::Array{T}, ::Val{N}) where {T,N}
size(x,1) == N || error("sizes mismatch")
isbitstype(T) || error("use for bitstypes only")
copy(reinterpret(SVector{N,T}, vec(x)))
end
function bstarCalc(ats::SVector{2, Float64},rΣ)
SVector{2,Float64}((hcat(ats,ats)+rΣ)\ats)
end
function SCalc(ats::SVector{2, Float64},bs::SVector{2, Float64})
SVector{1,Float64}(0.5*ats'*bs)
end
function eltorows(v::Vector{SVector{dim, T}}) where {dim, T}
copy(transpose(reshape(reinterpret(Float64, v), dim, :)))
end
function SumSVec(v::SVector{2, Float64})
SVector{1,Float64}(sum(v))
end
function MultSVecs(v1::SVector{dim, T}, v2::SVector{dim, T}) where {dim, T}
SVector{dim,Float64}(v1.*v2)
end
function ConcatSVecs(v1::Float64, v2::Float64)
SVector{2,Float64}(StaticArrays.SA[v1,v2])
end
function a0abeps(a0::SVector{1, Float64},ab::SVector{2, Float64},eps::SVector{2, Float64})
SVector{2,Float64}(a0 .+ ab .+ eps)
end
function absum(a::SVector{1, Float64},b::SVector{1, Float64},mp::Float64)
SVector{1,Float64}(a .+ b .+ mp)
end
function BMCalc(b::SVector{2, Float64},m::SVector{2, Float64})
SVector{1,Float64}(b'*m)
end
struct SimOutputs
down_match :: Vector{SVector{3,Float64}}
wages_match :: Vector{SVector{1,Float64}}
measures_match :: Vector{SVector{2,Float64}}
end
function SimRawData(b::Vector{Float64})
CR = CompParams();
Raw = StrucParams(b);
# Types
# Upstream observed: x_1, x_2
# Upstream unobs: delta^m
# Downstream Observed: y_1, y_2
# Downstream Unobs: delta^f
up_data = zeros(2, CR.n_firms);
up_data[1,:] = rand(Uniform(0.0,1.0), CR.n_firms);
up_data[2,:] = rand(Uniform(0.0,1.0), CR.n_firms);
down_data = zeros(2, CR.n_firms);
down_data[1,:] = rand(Uniform(0.0,1.0), CR.n_firms);
down_data[2,:] = rand(Uniform(0.0,1.0), CR.n_firms);
up_all = zeros(3, CR.n_firms);
down_all = zeros(3, CR.n_firms);
up_all[1:2,:] = up_data;
down_all[1:2,:] = down_data;
up_all[3,:] = rand(Normal(0.0,Raw.Σ_m[1,1]), CR.n_firms);
down_all[3,:] = rand(Normal(0.0, Raw.Σ_f[1,1]), CR.n_firms);
eps_measure = rand(MvNormal([0.0, 0.0],Raw.Σ),CR.n_firms);
alpha_0 = up_all[1:2,:]'*Raw.theta_0*down_all[1:2,:];
alpha_1 = up_all[1:2,:]'*Raw.theta_1*down_all[1:2,:];
alpha_2 = up_all[1:2,:]'*Raw.theta_2*down_all[1:2,:];
alpha_tilde = zeros(2,CR.n_firms,CR.n_firms);
alpha_tilde[1,:,:] = alpha_0+alpha_1;
alpha_tilde[2,:,:] = alpha_0+alpha_2;
bstar = zeros(2,CR.n_firms,CR.n_firms);
for n1 in 1:CR.n_firms;
for n2 in 1:CR.n_firms;
bstar[:,n1,n2] = (alpha_tilde[:,n1,n2]*ones(1,2).+(Raw.risk*Raw.Σ))\alpha_tilde[:,n1,n2]
end
end
S = zeros(CR.n_firms,CR.n_firms);
for n1 in 1:CR.n_firms;
for n2 in 1:CR.n_firms;
S[n1,n2] = 0.5*alpha_tilde[:,n1,n2]'*bstar[:,n1,n2];
end
end
#set up unobserved heterogeneity rhos:
#assume rho_m = y(1)*rm(1)*delta_m + y(2)*rm(2)*delta_m + rm(3)*delta_m
#assume rho_f = x(1)*rf(1)*delta_f + x(2)*rf(2)*delta_f + rf(3)*delta_f
rho_m = ([down_all[1:2,:]' ones(CR.n_firms,1)]* Raw.rm * up_all[3,:]')';
rho_f = [up_all[1:2,:]' ones(CR.n_firms,1)]* Raw.rf * down_all[3,:]';
C = -S .- rho_m .- rho_f; #negative of pairwise quadratic surplus
match, up_profit_lp, down_profit_lp = find_best_assignment(C)
down_match_lp= Array{Float64, 2}(undef, 3, CR.n_firms);
for i in 1:CR.n_firms
down_match_lp[:,i] = down_all[:, match[i][2]];
end
down_match_profit_lp = Array{Float64, 1}(undef, CR.n_firms);
for i in 1:CR.n_firms
down_match_profit_lp[i] = down_profit_lp[match[i][2]];
end
# minimum downstream profit =0
profit_diff = findmin(down_match_profit_lp)[1];
down_match_profit_lp .= down_match_profit_lp .+ profit_diff;
up_profit_lp .= up_profit_lp .- profit_diff;
alpha_0_match = diag(up_all[1:2,:]'*Raw.theta_0*down_match_lp[1:2,:]);
alpha_1_match = diag(up_all[1:2,:]'*Raw.theta_1*down_match_lp[1:2,:]);
alpha_2_match = diag(up_all[1:2,:]'*Raw.theta_2*down_match_lp[1:2,:]);
alpha_tilde_match = zeros(2,CR.n_firms);
alpha_tilde_match[1,:] = alpha_0_match.+alpha_1_match;
alpha_tilde_match[2,:] = alpha_0_match.+alpha_2_match;
bstar_match = zeros(2,CR.n_firms);
measures_match = zeros(2,CR.n_firms);
for n1 in 1:CR.n_firms;
bstar_match[:,n1] = (alpha_tilde_match[:,n1]*ones(1,2)+Raw.risk*Raw.Σ)\alpha_tilde_match[:,n1];
#simulate matched measures
measures_match[1,n1] = sum(alpha_0_match[n1]*bstar_match[:,n1]) + alpha_1_match[n1]*bstar_match[1,n1] + eps_measure[1,n1];
measures_match[2,n1] = sum(alpha_0_match[n1]*bstar_match[:,n1]) + alpha_2_match[n1]*bstar_match[2,n1] + eps_measure[2,n1];
end
S_match = zeros(CR.n_firms);
for n1 in 1:CR.n_firms;
S_match[n1] = 0.5*alpha_tilde_match[:,n1]'*bstar_match[:,n1]
end
rho_m_match = diag(([down_match_lp[1:2,:]' ones(CR.n_firms,1)]* Raw.rm * up_all[3,:]')');
rho_f_match = diag([up_all[1:2,:]' ones(CR.n_firms,1)]* Raw.rf * down_match_lp[3,:]');
up_valuation = -S_match .+ rho_m_match; #manager gets -S_match + unobs value, net of transfers
# eq transfers
up_prices_lp = up_profit_lp .- up_valuation;
#down_prices_lp= γ'*up_prices_lp;
#simulate wages:
wages_match = zeros(CR.n_firms)
for n1 in 1:CR.n_firms;
wages_match[n1] = Raw.mean_price + up_prices_lp[n1] + bstar_match[:,n1]'*measures_match[:,n1]
end
down_match_ret = down_match_lp[1:2,:]
out = (
up_data = up_data,
down_match_ret = down_match_ret,
wages_match = wages_match,
measures_match = measures_match)
return out
end
function sim_data_like(up_data_obs_in,down_data_obs_in,p_in::NamedTuple,i)
CR = CompParams();
# Types
# Upstream x_1, x_2, delta^m
# Downstream y_1, y_2, delta^f
# Generate simulated x and y data
Random.seed!(1234+i)
up_data_sim = Matrix{Float64}(undef, 3, CR.n_firms);
up_data_sim[1:2,:] = up_data_obs_in;
down_data_sim = similar(up_data_sim);
down_data_sim[1:2,:] = down_data_obs_in;
up_data_sim[3,:] = rand(Normal(0.0,p_in.Σ_u[1,1]), CR.n_firms);
down_data_sim[3,:] = rand(Normal(0.0, p_in.Σ_u[2,2]), CR.n_firms);
# noise to measures
eps_measure_sim = Vector{SVector{2,Float64}}(undef,CR.n_firms)
eps_measure_sim = DrawUnobsData(CR.n_firms , p_in.Σ)
alpha_0_sim = SMatrix{CR.n_firms,CR.n_firms,Float64}
alpha_1_sim = similar_type(alpha_0_sim)
alpha_2_sim = similar_type(alpha_0_sim)
alpha_0_sim = SMatrix{CR.n_firms,2,Float64}(up_data_obs_in')*SMatrix{2,2,Float64}(p_in.theta_0)*SMatrix{2,CR.n_firms,Float64}(down_data_obs_in);
alpha_1_sim = SMatrix{CR.n_firms,2,Float64}(up_data_obs_in')*SMatrix{2,2,Float64}(p_in.theta_1)*SMatrix{2,CR.n_firms,Float64}(down_data_obs_in);
alpha_2_sim = SMatrix{CR.n_firms,2,Float64}(up_data_obs_in')*SMatrix{2,2,Float64}(p_in.theta_2)*SMatrix{2,CR.n_firms,Float64}(down_data_obs_in);
alpha01 = similar_type(alpha_0_sim)
alpha02 = similar_type(alpha_0_sim)
alpha01 = alpha_0_sim + alpha_1_sim
alpha02 = alpha_0_sim + alpha_2_sim
alpha_tilde_simVec = Vector{SVector{2,Float64}}(undef,CR.n_firms*CR.n_firms)
alpha_tilde_simVec .= ConcatSVecs.(vec(alpha01),vec(alpha02))
bstar_sim = similar(alpha_tilde_simVec)
rΣ = Matrix{Float64}(undef,2,2)
rΣ = p_in.risk * p_in.Σ
bstar_sim .= bstarCalc.(alpha_tilde_simVec, Ref(rΣ))
@time S_simvec = Vector{SVector{1,Float64}}(undef, CR.n_firms*CR.n_firms)
S_simvec .= SCalc.(alpha_tilde_simVec,bstar_sim)
S_sim = SMatrix{CR.n_firms,CR.n_firms,Float64}
S_sim = SMatrix{CR.n_firms,CR.n_firms,Float64}(reshape(eltorows(S_simvec), (CR.n_firms,CR.n_firms)))
rho_m_sim = SMatrix{CR.n_firms,CR.n_firms,Float64}
rho_f_sim = similar_type(rho_m_sim)
rho_m_sim = SMatrix{CR.n_firms,1,Float64}(up_data_sim[3,:]')*SMatrix{1,3,Float64}(p_in.rm)*SMatrix{3,CR.n_firms,Float64}([down_data_obs_in; ones(1, CR.n_firms)]);
rho_f_sim = SMatrix{CR.n_firms,3,Float64}([up_data_obs_in' ones(CR.n_firms,1)])*SMatrix{3,1,Float64}(p_in.rf')*SMatrix{1,CR.n_firms,Float64}(down_data_sim[3,:]');
C_sim = SMatrix{CR.n_firms,CR.n_firms,Float64}
C_sim = -S_sim - rho_m_sim - rho_f_sim; #negative of pairwise quadratic surplus
match_sim, up_profit_sim, down_profit_sim = find_best_assignment(C_sim)
down_match_sim = similar(up_data_sim);
for i in 1:CR.n_firms
down_match_sim[:,i] = down_data_sim[:, match_sim[i][2]];
end
down_match_profit_sim = Array{Float64, 1}(undef, CR.n_firms);
for i in 1:CR.n_firms
down_match_profit_sim[i] = down_profit_sim[match_sim[i][2]];
end
# minimum downstream profit =0
profit_diff_sim = findmin(down_match_profit_sim)[1];
down_match_profit_sim .= down_match_profit_sim .+ profit_diff_sim;
up_profit_sim .= up_profit_sim .- profit_diff_sim;
alpha_0_match_sim = SVector{CR.n_firms,Float64}
alpha_1_match_sim = similar_type(alpha_0_match_sim)
alpha_2_match_sim = similar_type(alpha_0_match_sim)
alpha_0_match_sim = diag(SMatrix{CR.n_firms,2,Float64}(up_data_obs_in')*SMatrix{2,2,Float64}(p_in.theta_0)*SMatrix{2,CR.n_firms,Float64}(down_match_sim[1:2,:]));
alpha_1_match_sim = diag(SMatrix{CR.n_firms,2,Float64}(up_data_obs_in')*SMatrix{2,2,Float64}(p_in.theta_1)*SMatrix{2,CR.n_firms,Float64}(down_match_sim[1:2,:]));
alpha_2_match_sim = diag(SMatrix{CR.n_firms,2,Float64}(up_data_obs_in')*SMatrix{2,2,Float64}(p_in.theta_2)*SMatrix{2,CR.n_firms,Float64}(down_match_sim[1:2,:]));
alpha01_match = similar_type(alpha_0_match_sim)
alpha02_match = similar_type(alpha_0_match_sim)
alpha01_match = alpha_0_match_sim + alpha_1_match_sim
alpha02_match = alpha_0_match_sim + alpha_2_match_sim
alpha_tilde_match_simVec = Vector{SVector{2,Float64}}(undef,CR.n_firms)
alpha_tilde_match_simVec .= ConcatSVecs.(vec(alpha01_match),vec(alpha02_match))
alpha_0_match_simVec = copy(reinterpret(SVector{1,Float64}, vec(alpha_0_match_sim)))
bstar_match_sim = similar(alpha_tilde_match_simVec)
bstar_match_sim .= bstarCalc.(alpha_tilde_match_simVec, Ref(rΣ))
sumBstar = Vector{SVector{1,Float64}}(undef, CR.n_firms)
sumBstar .= SumSVec.(bstar_match_sim)
a0sumBstar = Vector{SVector{1,Float64}}(undef, CR.n_firms)
a0sumBstar .= MultSVecs.(alpha_0_match_simVec,sumBstar)
measures_match_sim = similar(bstar_match_sim)
a12 = similar(bstar_match_sim)
a12 .= ConcatSVecs.(vec(alpha_1_match_sim),vec(alpha_2_match_sim))
ab = similar(bstar_match_sim)
ab .= MultSVecs.(a12,bstar_match_sim)
measures_match_sim .= a0abeps.(a0sumBstar,ab,eps_measure_sim)
S_match_simvec = Vector{SVector{1,Float64}}(undef, CR.n_firms)
S_match_simvec .= SCalc.(alpha_tilde_match_simVec,bstar_match_sim)
rho_m_match_sim = diag(SMatrix{CR.n_firms,1,Float64}(up_data_sim[3,:]')*SMatrix{1,3,Float64}(p_in.rm)*SMatrix{3,CR.n_firms,Float64}([down_match_sim[1:2,:]; ones(1, CR.n_firms)]));
rho_m_match_simSvec = copy(reinterpret(SVector{1,Float64}, rho_m_match_sim))
#rho_f_match_sim = diag([up_data_sim[1:2,:]' ones(CR.n_firms,1)]* p_in.rf * down_match_sim[3,:]');
up_valuation_sim = -S_match_simvec .+ rho_m_match_simSvec; #manager gets -S_match + unobs value, net of transfers
up_profit_simSvec = copy(reinterpret(SVector{1,Float64}, up_profit_sim))
# eq transfers
up_prices_sim = up_profit_simSvec .- up_valuation_sim;
#down_prices_sim= γ_sim'*up_prices_sim;
#simulate wages:
wages_match_sim = Vector{SVector{1,Float64}}(undef, CR.n_firms)
bm = Vector{SVector{1,Float64}}(undef, CR.n_firms)
bm .= BMCalc.(bstar_match_sim,measures_match_sim)
wages_match_sim .= absum.(up_prices_sim, bm,p_in.mean_price);
down_match_sim_out = Vector{SVector{3,Float64}}(undef, CR.n_firms)
down_match_sim_out = svectorscopy(down_match_sim, Val{3}())
down_match = down_match_sim_out
wages_match = wages_match_sim
measures_match = measures_match_sim
Simout = SimOutputs(down_match, wages_match, measures_match)
return Simout
end
#liklihood function
function loglikepr(b_est,b_cal,i_cal,up_data_obs,down_data_obs,wages_obs,measures_obs)
# inputs are:
# b_est = parameter values of parameters to be estimated
# b_cal = parameters values of parameters to be calibrated
# i_cal = index in b for parameters to be calibrated
# up_data_obs - the observable part of up_data
# down_data_obs - the obs part of down_data
# wages_obs - obs wages
# measures_obs = obs measures
LLC = CompParams();
ic = 1;
ie = 1;
nb = length(b_est)+length(b_cal);
b = zeros(nb);
for ib in 1:nb
if ic<=length(i_cal)
if ib == i_cal[ic];
b[ib] = b_cal[ic];
ic += 1;
else
b[ib] = b_est[ie];
ie += 1;
end
else
b[ib] = b_est[ie];
ie += 1;
end
end
LLp = StrucParams(b);
@time sim_dat = Vector{SimOutputs}(undef,LLC.n_sim)
sim_dat = [sim_data_like(up_data_obs,down_data_obs,LLp,314159+i) for i in eachindex(sim_dat)]
ll=0.0;
# Bandwidths for y, wages, measures, cov of measures
data = zeros(5,LLC.n_firms*LLC.n_sim);
for i = 1:LLC.n_sim
data[1,1+(i-1)*LLC.n_firms:i*LLC.n_firms] = view(down_data_obs,1,:) .- getindex.(sim_dat[i].down_match,1);
data[2,1+(i-1)*LLC.n_firms:i*LLC.n_firms] = view(down_data_obs,2,:) .- getindex.(sim_dat[i].down_match,2);
data[3,1+(i-1)*LLC.n_firms:i*LLC.n_firms] = (wages_obs) .- getindex.(sim_dat[i].wages_match,1);
data[4,1+(i-1)*LLC.n_firms:i*LLC.n_firms] = view(measures_obs,1,:) .- getindex.(sim_dat[i].measures_match,1);
data[5,1+(i-1)*LLC.n_firms:i*LLC.n_firms] = view(measures_obs,2,:) .- getindex.(sim_dat[i].measures_match,2);
end
if (LLC.hmethod==1)
#Silverman's Rule of thumb:
Iqrnorm = 1.34;
A = zeros(5);
A[1] = min(std(data[1,:]),StatsBase.iqr(data[1,:])/Iqrnorm);
A[2] = min(std(data[2,:]),StatsBase.iqr(data[2,:])/Iqrnorm);
A[3] = min(std(data[3,:]),StatsBase.iqr(data[3,:])/Iqrnorm);
A[4] = min(std(data[4,:]),StatsBase.iqr(data[4,:])/Iqrnorm);
A[5] = min(std(data[5,:]),StatsBase.iqr(data[5,:])/Iqrnorm);
h=0.9 * A * LLC.n_firms^(-0.2);
else
#Using Sheather and Jones (1991) bandwidth
h = zeros(5);
h[1] = KernelDensitySJ.bwsj(data[1,:]);
h[2] = KernelDensitySJ.bwsj(data[2,:]);
h[3] = KernelDensitySJ.bwsj(data[3,:]);
h[4] = KernelDensitySJ.bwsj(data[4,:]);
h[5] = KernelDensitySJ.bwsj(data[5,:]);
end
like = zeros(LLC.n_firms);
for i in 1:LLC.n_firms
for j in 1:LLC.n_sim
like[i]+=(
pdf(Normal(),(down_data_obs[1,i] - getindex(sim_dat[j].down_match[i],1))/h[1])
*pdf(Normal(),(down_data_obs[2,i] - getindex(sim_dat[j].down_match[i],2))/h[2])
*pdf(Normal(),(wages_obs[i] - getindex(sim_dat[j].wages_match[i],1))/h[3])
*pdf(Normal(),(measures_obs[1,i] - getindex(sim_dat[j].measures_match[i],1))/h[4])
*pdf(Normal(),(measures_obs[2,i] - getindex(sim_dat[j].measures_match[i],2))/h[5])
)
end
end
if LLC.trim_percent > 0
#set up trimming function following Fermanian and Salanie
hall = prod(h);
delta_trim = log(percentile(like,LLC.trim_percent))/log(hall);
if isinf(abs(delta_trim)) #if above trim fn yields computational zeros, then we need less smooth trim fn
hd = percentile(like,LLC.trim_percent)
tau_trim =zeros(LLC.n_firms)
tau_trim .= (like.-hd)./(2.0*hd)
tau_trim[like .< hd] .=0.0
tau_trim[like .> 2.0*hd] .=1.0
maxtau = findmax(tau_trim)
else
hd = hall^delta_trim; #ensure hd, hd3, hd4 are not computationally zero
hd3 = hd^3.0;
hd4 = hd^4.0;
tau_trim = zeros(LLC.n_firms)
tau_trim .= 4 .*(like.-hd).^3 ./ hd3 .- 3.0 .*(like.-hd).^4.0 ./ hd4;
tau_trim[like .< hd] .=0.0
tau_trim[like .> 2.0*hd] .=1.0
maxtau = findmax(tau_trim)
end
else
tau_trim = zeros(LLC.n_firms)
tau_trim.=1.0
maxtau = findmax(tau_trim)
end
ll = 0.0
if tau_trim==0.0 || isnan(maxtau[1])
ll = NaN
else
for i in 1:LLC.n_firms
if tau_trim[i]>0.0
ll+=tau_trim[i]*log(like[i]/(LLC.n_sim*h[1]*h[2]*h[3]*h[4]*h[5]))
end
end
end
println("parameter: ", round.(b_est, digits=5), " likelihood: ", ll/LLC.n_firms)
return ll/LLC.n_firms
end
function SetTrueParams()
# Set up true parameters
b_true = zeros(25);
#unobs heterogeneity
Σ_m_true = exp(0.0);
b_true[1] = log(Σ_m_true);
Σ_f_true = exp(0.0);
b_true[2] = log(Σ_f_true);
#Set up thetas
theta_0_true = [0.5 0.7;
0.3 0.1];
b_true[3]= theta_0_true[1,1];
b_true[4]= theta_0_true[1,2];
b_true[5]= theta_0_true[2,1];
b_true[6]= theta_0_true[2,2];
theta_1_true = [0.2 0.05;
0.2 0.4];
b_true[7]= theta_1_true[1,1];
b_true[8]= theta_1_true[1,2];
b_true[9]= theta_1_true[2,1];
b_true[10]= theta_1_true[2,2];
theta_2_true = [0.1 0.0;
0.0 0.2];
b_true[11]= theta_2_true[1,1];
b_true[12]= theta_2_true[1,2];
b_true[13]= theta_2_true[2,1];
b_true[14]= theta_2_true[2,2];
#set up Cov matrix for shocks to measures
sigma_1_true = 0.3;
b_true[15] = log(sigma_1_true);
sigma_2_true = 0.1;
b_true[17] = log(sigma_2_true);
#cov btw 1 &2:
b_true[16]= 1.02;
cov12_true = -sqrt(sigma_1_true)sqrt(sigma_2_true) + 2*sqrt(sigma_1_true)sqrt(sigma_2_true)*(exp(b_true[16]))/(1+exp(b_true[16])); #ensuring pos semi def cov matrix
mean_price_true = 0.0;
b_true[18] = mean_price_true;
rm_true = [0.1 -0.3 0.2];
b_true[19:21] = rm_true;
rf_true = [0.0 0.2 0.1];
b_true[22:24] = rf_true;
risk_true = 0.05;
b_true[25] = log(risk_true);
return b_true
end
b_true = SetTrueParams();
Rawout = SimRawData(b_true);
igrid = Vector{Integer}(1:length(b_true));
ical = deleteat!(igrid,1);
b_cal = b_true[ical];
b_est = [b_true[1]+.01;
@time loglikepr(b_est,b_cal,ical,Rawout.up_data,Rawout.down_match_ret, Rawout.wages_match, Rawout.measures_match)
```