mindv_monte1.jl

using MinimumDivergence
using ModelsGenerators
using Ipopt
using Divergences

HellingerDistance()  = CressieRead(0.5)
ContinuousUpdating() = CressieRead(2.0)

ModifiedHellingerDistance(ϑ::Real)  = ModifiedCressieRead(0.5, ϑ)

Divergences.KullbackLeibler(::Int64) = KullbackLeibler()

stddiv = [ :KullbackLeibler, :ReverseKullbackLeibler, :HellingerDistance,
          :ContinuousUpdating ];

moddiv = [ :ModifiedKullbackLeibler, :ModifiedReverseKullbackLeibler, :ModifiedHellingerDistance]

stdprob = [:KL,  :RKL,  :HD,  :CUE]

modprob = [:MKL, :MRKL, :MHD]

g(θ) = z.*(y-x*θ)

lb = [-1.0, -20.0]
ub = [1.0, 20.0]
θ₀ = [0.0, 1.0]

dgp(;args...) = randiv_ts(;args...)
y, x, z = dgp()

for (fname, fdiv) in zip(stdprob, stddiv)
    eval( quote
        ($fname) = MinDivProb(
                              MomentFunction(g, :dual, IdentityKernel(), nobs = size(x,1), nmom = size(z, 2), npar = size(x,2)),
                              ($fdiv)(), θ₀, lb, ub, solver = IpoptSolver(print_level = 0))
    end)
end


for (fname, fdiv) in zip(modprob, moddiv)
    eval( quote
        ($fname) = MinDivProb(
                              MomentFunction(g, :dual, IdentityKernel(), nobs = size(x,1), nmom = size(z, 2), npar = size(x,2)),
                              ($fdiv)(.05), θ₀, lb, ub, solver = IpoptSolver(print_level = 0))
    end)
end


coefficients = [i => Array(Float64, 1)  for i = stdprob]
stderrors    = [i => Array(Float64, 1)  for i = stdprob]

coefficients = [i => Float64[]  for i = stdprob]
stderrors    = [i => Float64[]  for i = stdprob]



finalize(x::MinDivProb) = [coef(x), MinimumDivergence.stderr(x)]

plain_status(mdp::MinDivProb) = mdp.model.inner.status

ms_thetas = [[rand(Distributions.Uniform(-.1,.1)), rand(Distributions.Uniform(-5,5))] for i = 1:5]

function sim(sdata)
    y[:] = sdata[1]
    x[:] = sdata[2]
    z[:] = sdata[3]
    probnames = [stdprob, modprob]
    for fname in probnames
        qex = Expr(:quote, fname)
        eval(quote
            MinimumDivergence.multistart($fname, ms_thetas)
        end)
    end

    uu =  map(probnames) do x
        :(coef($x))
    end

    out = float([eval(uu[i])[2] for i = 1:7])

    uu =  map(probnames) do x
        :(MinimumDivergence.stderr($x))
    end

    out = [out, float([eval(uu[i])[2] for i = 1:7])]

    uu =  map(probnames) do x
        :(plain_status($x))
    end

    out = [out, float([eval(uu[i]) for i = 1:7])]

    return out'

end



srand(1)
jj = vcat([sim(dgp()) for i = 1:30]...)

srand(1)
jj = vcat([sim(dgp()) for i = 1:100]...)

df = convert(DataFrame, coefficients)
describe(df)

df = convert(DataFrame, stderrors)
describe(df)