I need to estimate parameters (a,b) and their corresponding standard errors using simulated MLE. I can successfully estimate the parameters, but I fail to get the SEs. In particular, when I try to obtain the hessian from the optimization using automatic differentiation, it gives me a zero matrix. Since the SE is the square root of the inverse of the hessian, I am unable to compute the SEs using this method. Can you please let me know if there is something wrong with how I am calling the hessian?
I first present the algorithm in words, I then present the implementation in Julia.
The algorithm for the SMLE is as follows:
-
We are given data for 10 types of individuals i for some share s^{data}_i and some variable Y_i. We need the relationship between s^{data}_i and Y_i. To do so we execute the following:
-
For each type of individual i, simulate Y_{ij}, u_{ij} for 10,000 individuals j as follows. For each j, draw (i) Y_{ij} according to some pre-determined rule, and (ii) make a uniform draw u_{i,j} between 0 and 1.
-
Construct logit probabilities for each i,j as follows p_{i,j} = \frac{\exp (a + b\cdot Y_{i,j})}{1+\exp (a + b\cdot Y_{i,j})}
-
Assign s_{i,j}=1 if p_{i,j} \geq u_{i,j}, set s_{i,j} = 0 otherwise
-
Construct the simulated likelihood contribution for type i as l_i = \frac{1}{10000}\sum_j s_{i,j}
-
Construct the overall log likelihood as the weighted sum L = \sum_i \ s^{data}_i \times \ln l_i + (1-s^{data}_i )\ln(1-l_i)
-
Search for the parameters (a,b) which minimize the simulated likelihood and obtain their standard errors. I can estimate the parameters in Julia. However the SEs are undefined because Julia returns a zero matrix for the hessian which cannot be inverted.
Below is the implementation in Julia.
using Random, LinearAlgebra, Distributions, Optim, StatsBase
using Optim: converged, maximum, maximizer, minimizer, iterations
using DelimitedFiles
using Statistics
using SharedArrays
using Optim, NLSolversBase, Random
using LinearAlgebra: diag
Random.seed!(1234)
## define fake data and make some simulations for the estimation
draws_uniform = rand(Uniform(0,1), 10, 10000)
draws_y = cumsum(ones(Float64, 10)) .+ rand(Uniform(1,2), 10, 10000)
share_hasgood = [.1 .2 .3 .4 .5 .6 .7 .8 .9 .9]'
## construct share by (i) comparing logit probability p to draws_uniform, and (ii) summing over all columns with p > draws_uniform
function model_derived_share(t, draws_uniform, draws_y)
a = t[1]
b = t[2]
p = exp.(a .+ b .* draws_y) ./ (1 .+ exp.(a .+ b .* draws_y)) # probability indiv has good
s = p .> draws_uniform # assign indiv as having or not having good by comparing model probability with uniform draw
return mean(s, dims=2)
end
# construct log likelihood by taking weighted sum of model-derived shares
loglike(model_prediction, actual_share) = -sum(log.(model_prediction) .* actual_share .+ (1 .- actual_share) .* log.(1 .- model_prediction))
# optimize by minimizing negative of the log likelihood
using Optim
params = Optim.minimizer(optimize((t) -> loglike(model_derived_share(t, draws_uniform, draws_income), share_hasgood), [.5; .5]))
display(params)
# evaluate hessian of loglikelihood at optimum params = [a*;b*]
functwicediff = TwiceDifferentiable(t -> momentcond(model_derived_share(t, draws_uniform, draws_y), share_hasgood), [.5; .5]; autodiff=:forward)
numerical_hessian = hessian!(functwicediff, params)
display(numerical_hessian)
# compute standard error of parameters as square root of the inverse of the hessian
std_err = sqrt.(diag(inv(numerical_hessian)))