floswald / smm.jl Goto Github PK

View Code? Open in Web Editor NEW

73.0 12.0 25.0 12.92 MB

Simulated Method of Moments for Julia

Home Page: https://floswald.github.io/SMM.jl/latest

License: Other

Julia 84.13% Shell 0.55% Jupyter Notebook 15.32%

julia mcmc estimation

smm.jl's Introduction

SMM.jl: Simulated Method of Moments for Julia

Documentation	Build Status

Notice: this package was previously called MomentOpt.jl.

This package provides a Julia infrastructure for Simulated Method of Moments estimation, or other problems where we want to optimize a non-differentiable objective function. The setup is suitable for all kinds of likelihood-free estimators - in general, those require evaluating the objective at many regions. The user can supply their own algorithms for generating successive new parameter guesses. We provide a set of MCMC template algorithms. The code can be run in serial or on a cluster.

Installation

In your julia REPL, type

] add SMM

Documentation

Examples

Please check out a fully worked example in src/Examples.jl.

Here is a working session comparing serial vs parallel performance on a test objective function. Notice that parallel performance hinges on the objective function being reasonably expensive to compute (at least 0.1 seconds per function evaluation) - otherwise the overhead from data transfer is just too high.

julia> using SMM
[ Info: Precompiling SMM [bc769cb7-f249-5bba-802a-79f18cb247ec]

julia> x = SMM.serialNormal(2,200,slow = true)
[ Info: Starting estimation loop.
Progress: 100%|██████████████████████████████████████████████| Time: 0:01:05
┌ Warning: could not find 'filename' in algo.opts - not saving anything
└ @ SMM ~/.julia/dev/SMM/src/mopt/AlgoAbstract.jl:67
[ Info: Done with estimation after 1.1 minutes
summary(MA) = 3×5 DataFrame
 Row │ id     acc_rate   perc_exchanged  exchanged_most_with  best_val
     │ Int64  Float64    Float64         Int64                Float64
─────┼──────────────────────────────────────────────────────────────────
   1 │     1  0.0793651             5.5                    2  0.0023224
   2 │     2  0.0819672             8.5                    1  0.0126754
   3 │     3  0.115183              4.5                    2  0.0145625
(
BGP Algorithm with 3 BGPChains
============================

Algorithm
---------
Current iteration: 200
Number of params to estimate: 2
Number of moments to match: 2

, Plot{Plots.GRBackend() n=2}, Plot{Plots.GRBackend() n=3})

julia> using Distributed

julia> addprocs(2, exeflags="--project=.")  # you don't need the `exeflag` if you `add`ed the package regularly!
2-element Array{Int64,1}:
 2
 3

julia> @everywhere using SMM

julia> x = SMM.serialNormal(2,200,slow = true)
[ Info: Starting estimation loop.
Progress: 100%|█████████████████████████████████████████████| Time: 0:00:49
┌ Warning: could not find 'filename' in algo.opts - not saving anything
└ @ SMM ~/.julia/dev/SMM/src/mopt/AlgoAbstract.jl:67
[ Info: Done with estimation after 0.8 minutes
summary(MA) = 3×5 DataFrame
 Row │ id     acc_rate   perc_exchanged  exchanged_most_with  best_val
     │ Int64  Float64    Float64         Int64                Float64
─────┼───────────────────────────────────────────────────────────────────
   1 │     1  0.117347              2.0                    2  0.00246371
   2 │     2  0.0899471             5.5                    3  0.103399
   3 │     3  0.161458              4.0                    2  0.139263
(
BGP Algorithm with 3 BGPChains
============================

Algorithm
---------
Current iteration: 200
Number of params to estimate: 2
Number of moments to match: 2

, Plot{Plots.GRBackend() n=2}, Plot{Plots.GRBackend() n=3})

Contributing

We encourage user contributions. Please submit a pull request for any improvements you would like to suggest, or a new algorithm you implemented.

New algorithms:

You can model your algo on the basis of src/AlgoBGP.jl -
you need to implement the function computeNextIteration!( algo ) for your algo

History

This package grew out of the R package mopt.

Thanks to all Contributors!

smm.jl's People

Contributors

Stargazers

Watchers

Forkers

priscillafialho alancrawford julienpascal edoardociscato abolrezghi maximilianjhuber theloudis gregobad tyleransom lnsongxf augustbruno hdccer pedrovergaramerino jcaste94 standardgalactic tradler oliverevans96 junbiaochen1998 santerivtj

smm.jl's Issues

don't do exchange moves with replacement?

if i have 3 chains, all potential pairs are (1,2), (1,3), (2,3). Suppose i exchange 1 with 2 and 1 with 3. it becomes a bit hard to follow what happened to chain 1. to sample from all potential pairs with replacement is the exact translation of the paper.

1 <=> 2
1 now 2
2 now 1

1 <=> 3 (i.e. 2 <=> 3)
1 now 3 (i.e. 2 now 3) 
3 now 1 (i.e. 3 now 2)

result:
1 now 2
2 now 1
3 now 2

it may be more desirable to do this without replacement, i.e. once I select pair (1,2) I cannot exchange with either 1 or 2 any more? this reduces the amount of mixing in the chain.
the problem right now is that it says on chain 3 "echanged with 1", but in fact chain 1 was chain 2 at that point already. maybe this is irrelevant.

Error tagging new release

The tag name "v1.3" is not of the appropriate SemVer form (vX.Y.Z).
cc: @floswald

Error tagging new release

The tag name "v1.4" is not of the appropriate SemVer form (vX.Y.Z).
cc: @floswald

benchmark performance

could be faster
allocates a lot of memory

Pause and restart

Hi,

Thank you for the awesome work.
Is there any way to "pause" the estimation and to restart from the same point?

Cheers,
Julien

Re-scaling parameters to unit range

Hi,

It seems like it's harder for the algorithm to find the minimum when working with parameters of different order of magnitude (let's say 1 versus 0.01, as your examples in the README).

Maybe an easy fix would be to automatically re-scale each parameter to [0,1] (or [-1,1])?

Of course, when calling the objective function, we need to make sure that we transform back [0,1]^n to
its original value (n being the number of parameters to estimate). But the algorithm itself would be working only with [0,1]^n.

Does it make sense?

If yes, I can implement this on a new branch and make a PR.

Cheers

should we only track Array{Eval} ?

we could add a couple of more fields to the Eval type (like accept, prob, and a couple of other things we store in the infos dataframe now) and track only those. Currently, we set up the Eval anyway to do the computation, then store in dataframe, then unpick the dataframe to HDF5 arrays, then reassemble the dataframe at home from the HDF5. if we stored an Array{Eval} on each chain we could save that and cut out one layer.

warn macro needs logger arg; problem with pmap; saving MAlgoBGP not working

Hi @floswald
it's me again. I've been using the latest version of MomentOpt (thanks a lot for this package) and possibly found a couple of typos. Before sending a pull request, I wanted to make sure I'm not wrong.

It seems to me that the macro @warn also needs the logger argument. By the way, I found the package MiniLogging much better than Logging. See for instance (there are a couple of other cases)

https://github.com/floswald/MomentOpt.jl/blob/d50ef47b77773f5addfc1e31d52756c64a1a998d/src/mopt/AlgoAbstract.jl#L49

I had troubles implementing the BGP algorithm because of the following line

https://github.com/floswald/MomentOpt.jl/blob/d50ef47b77773f5addfc1e31d52756c64a1a998d/src/mopt/AlgoBGP.jl#L482

the function pmap does not store the output by replacing its argument algo.chains. In other words, there is no such thing as pmap!. The same is true for the non-parallelized version with map. As a result, the algorithm stopped after a single iteration because algo.chains[ic].iter was not updated on any chain. In order for the algorithm to work I had to replace that line with:

algo.chains = pmap( x->next_eval!(x), algo.chains )

Is this due to my Julia version? I'm using 0.6.0

I did not succeed in saving objects of type MAlgoBGP and it seems to me that it is because the function saveBGPChainToHDF5 is missing. I guess you might be still working on this.

Trying to add example that estimates a basic linear regression model using SMM

Hi @floswald, I'm getting more and more familiar with your package and I think it would be great to have some more examples (and, in particular, a regression example).

I've create a gist that attempts to do this. I will happily submit it as a PR to add to the examples if you're interested (and once it actually works!).

My example specifies the following moments to match:

E[y'X] in the data and E[y'X] in the model
V[y] in data and V[Xβ + ε] in the model (I specified this moment so as to be able to estimate the variance of the error term, like in the MSE formula)

This seems pretty easy to implement, but I am getting a couple of cryptic errors that I can't seem to resolve and was hoping it would be easy for you to spot where I'm going wrong:

Undefined error for dataMomentWd:

┌ Warning: caught exception UndefVarError(:dataMomentWd) at param OrderedDict(:b0 => -0.302835267638089,:b1 => -0.3927639716972684,:b2 => -0.7839915331115677,:b3 => 0.3630821477120001,:s => 1.474363363563423)`
└ @ SMM ~/.julia/packages/SMM/MDWS3/src/mopt/mprob.jl:167

No matter what I set for smpl_iters (it could be 100, 1_000 or 10_000), I always only get back 12 accepted draws

I couldn't see in the example files where dataMomentWd is defined, so that's why I'm getting the undefined error on that in my example.

I also had one other question about SMM in general (since this is my first foray into these types of models):

Does one usually simulate the X's in addition to drawing ε and using (X, β, ε) to generate y? The reason I ask is because the X's might have a different dimension than the number of draws of ε, so I can't seem to figure out how to resolve this. For example, suppose I wanted to estimate a regression using the iris data in R and this SMM.jl package. How would I allow for ns [link to a line in the aforementioned gist] to be larger than 150 (the number of rows in iris)?

Thanks a ton for all the help you've given me so far!

Is there a way to pause and restart?

Hi, thank you for the package. I need to submit the job on the HPC which has a 24 hour maximum time. So, I wonder if I can set a time limit on the optimization and if I can stop and restart the optimization? Thank you.

add taxicab search

need a simple taxicab search https://en.wikipedia.org/wiki/Coordinate_descent

replace minilogging

i can't stand the mutilplication of output streams with each reload of the pacakge. going back to Base.logging

very high dim makes MVNormal fail

ok, what about this (not 6 params, but 18!)

julia> MomentOpt.snorm_6()
INFO: These moments are our targets
INFO: Parameter p_i corresponds to moment m_i
18×3 DataFrames.DataFrame
│ Row │ name │ value │ weight │
├─────┼──────┼───────┼────────┤
│ 1   │ mu1  │ -1.0  │ -1.0   │
│ 2   │ mu2  │ 0.0   │ 0.0    │
│ 3   │ mu3  │ 0.5   │ 0.5    │
│ 4   │ mu4  │ -0.5  │ -0.5   │
│ 5   │ mu5  │ 0.008 │ 0.008  │
│ 6   │ mu6  │ -0.7  │ -0.7   │
│ 7   │ mu7  │ -1.0  │ -1.0   │
│ 8   │ mu8  │ 0.0   │ 0.0    │
│ 9   │ mu9  │ 0.5   │ 0.5    │
│ 10  │ mu10 │ -0.5  │ -0.5   │
│ 11  │ mu11 │ 0.008 │ 0.008  │
│ 12  │ mu12 │ -0.7  │ -0.7   │
│ 13  │ mu13 │ -1.0  │ -1.0   │
│ 14  │ mu14 │ 0.0   │ 0.0    │
│ 15  │ mu15 │ 0.5   │ 0.5    │
│ 16  │ mu16 │ -0.5  │ -0.5   │
│ 17  │ mu17 │ 0.008 │ 0.008  │
│ 18  │ mu18 │ -0.7  │ -0.7   │
21:52:41:INFO:Main:Starting estimation loop.
21:52:41:INFO:Main:Starting estimation loop.
21:52:41:INFO:Main:Starting estimation loop.
ERROR: no draw in support after 1000 trials. increase either opts[smpl_iters] or opts[bound_prob].
Stacktrace:
 [1] mysample(::Distributions.MvNormal{Float64,PDMats.PDiagMat{Float64,Array{Float64,1}},Array{Float64,1}}, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64) at /Users/74097/.julia/v0.6/MomentOpt/src/mopt/AlgoBGP.jl:342

Model simulation issue

Hi, Florian. One silly stupid question. To apply this code, do I need to simulate my model under every parameter combination (to obtain the empirical moments, I need dozens of simulations runs with 1000 simulated individuals each run), or alternatively, just one simulation at the initial time? As you know, each round of simulation take a lot of time, it's a crucial matter.

Eval not defined

I'm getting the following error when I try to run a (script) version of the example ipynb:

ERROR: LoadError: UndefVarError: Eval not defined

This is strange, given that Eval is part of the SMM package, and I'm invoking using SMM at the top of my script (included below in its entirety). The error shows up in the function objfunc_normal(ev::Eval; verbose = false) line

using SMM
using DataStructures
using DataFrames
using Plots
using Random

Random.seed!(1234);

# Define true values of parameters
#---------------------------------
trueValues = OrderedDict("mu1" => [-1.0], "mu2" => [1.0])
#------------------------------------------------
# Options
#-------------------------------------------------
# Boolean: do you want to save the plots to disk?
savePlots = false

#------------------------
# initialize the problem:
#------------------------

# Specify the initial values for the parameters, and their support:
pb = OrderedDict("p1" => [0.2,-3,3] , "p2" => [-0.2,-2,2])

# Specify moments to be matched + subjective weights:
moms = DataFrame(name=["mu1","mu2"],value=[trueValues["mu1"][], trueValues["mu2"][]], weight=ones(2))

# GMM objective function to be minized.
# It returns a weigthed distance between empirical and simulated moments
function objfunc_normal(ev::Eval; verbose = false)

    start(ev)


    # when running in parallel, display worker's id:
    #-----------------------------------------------
    if verbose == true
        if nprocs() > 1
          println(myid())
        end
    end

    # extract parameters from ev:
    #----------------------------
    mu  = collect(values(ev.params))

    # compute simulated moments
    #--------------------------
    # Monte-Carlo:
    #-------------
    ns = 100000 #number of i.i.d draws from N([mu], sigma)
    #initialize a multivariate normal N([mu], sigma)
    #sigma is set to be the identity matrix
    sigma = [1.0; 1.0]
    # draw ns observations from N([mu], sigma):
    randMultiNormal = SMM.MvNormal(mu,SMM.PDiagMat(sigma))
    # calculate the mean of the simulated data
    simM            = mean(rand(randMultiNormal,ns),2)
    # store simulated moments in a dictionary
    simMoments = Dict(:mu1 => simM[1], :mu2 => simM[2])

    # Calculate the weighted distance between empirical moments
    # and simulated ones:
    #-----------------------------------------------------------
    v = Dict{Symbol,Float64}()
    for (k, mom) in dataMomentd(ev)
        # If weight for moment k exists:
        #-------------------------------
        if haskey(SMM.dataMomentWd(ev), k)
            # divide by weight associated to moment k:
            #----------------------------------------
            v[k] = ((simMoments[k] .- mom) ./ SMM.dataMomentW(ev,k)) .^2
        else
            v[k] = ((simMoments[k] .- mom) ) .^2
        end
    end

    # Set value of the objective function:
    #------------------------------------
    setValue(ev, mean(collect(values(v))))

    # also return the moments
    #-----------------------
    setMoment(ev, simMoments)

    # flag for success:
    #-------------------
    ev.status = 1

    # finish and return
    finish(ev)

    return ev
end



# Initialize an empty MProb() object:
#------------------------------------
mprob = MProb()

# Add structural parameters to MProb():
# specify starting values and support
#--------------------------------------
addSampledParam!(mprob,pb)

# Add moments to be matched to MProb():
#--------------------------------------
addMoment!(mprob,moms)

# Attach an objective function to MProb():
#----------------------------------------
addEvalFunc!(mprob, objfunc_normal)
# estimation options:
#--------------------
# number of iterations for each chain
niter = 1000
# number of chains
nchains = 3

opts = Dict("N"=>nchains,
        "maxiter"=>niter,
        "maxtemp"=> 5,
        "coverage"=>0.025,
        "sigma_update_steps"=>10,
        "sigma_adjust_by"=>0.01,
        "smpl_iters"=>1000,
        "parallel"=>true,
        "maxdists"=>[0.05 for i in 1:nchains],
        "mixprob"=>0.3,
        "acc_tuner"=>12.0,
        "animate"=>false)


#---------------------------------------
# Let's set-up and run the optimization
#---------------------------------------
# set-up BGP algorithm:
MA = MAlgoBGP(mprob,opts)

# run the estimation:
@time SMM.run!(MA)

# show a summary of the optimization:
@show SMM.summary(MA)
# # Plot histograms for the first chain, the one with which inference should be done.
# # Other chains are used to explore the space and avoid local minima
# #-------------------------------------------------------------------------------
# p1 = histogram(MA.chains[1])
# display(p1)
# 
# if savePlots
#     savefig(p1, joinpath(pwd(),"histogram.svg"))
# end
# 
# # Plot the realization of the first chain
# #----------------------------------------
# p2 = plot(MA.chains[1])
# if savePlots
#     savefig(p2, joinpath(pwd(),"lines.svg"))
# end
# display(p2)


# Realization of the first chain:
#-------------------------------
dat_chain1 = SMM.history(MA.chains[1])

# discard the first 10th of the observations ("burn-in" phase):
#--------------------------------------------------------------
dat_chain1[round(Int, niter/10):niter, :]

# keep only accepted draws:
#--------------------------
dat_chain1 = dat_chain1[dat_chain1[:accepted ].== true, : ]

# create a list with the parameters to be estimated
parameters = [Symbol(String("mu$(i)")) for i=1:2]
# list with the corresponding priors:
#------------------------------------
estimatedParameters = [Symbol(String("p$(i)")) for i=1:2]


# Quasi Posterior mean and quasi posterior median for each parameter:
#-------------------------------------------------------------------
for (estimatedParameter, param) in zip(estimatedParameters, parameters)

  println("Quasi posterior mean for $(String(estimatedParameter)) = $(mean(dat_chain1[estimatedParameter]))")
  println("Quasi posterior median for $(String(estimatedParameter)) = $(median(dat_chain1[estimatedParameter]))")
  println("True value = $(trueValues[String(param)][])")
  println("abs(True value - Quasi posterior median)  = $(abs(median(dat_chain1[estimatedParameter]) - trueValues[String(param)][])) \n")

end

API: pass opts in another way

the current setup with a dict is no good
too easy to forget an option
or to set something whcih has no counterpart in the code.

think about other ways of setting opts and passing them to the chain constructor

replace HDF5 by JLD2

upgrade to v0.6

Minor issue with field best_id in BGPChain

https://github.com/floswald/MomentOpt.jl/blob/0c90a769d42aad3c29d1f5aecf5a49cd1e63f406/src/mopt/AlgoBGP.jl#L179-L199

As you can see in line 197 c.best_id[c.iter] is set to be equal to c.iter even when the new value at iteration c.iter is not the best value in the chain. Once the algorithm is done, I get that c.best_id is simply 1:opts["maxiter"]. I guess c.best_id[c.iter] should probably be set to c.best_id[c.iter-1] (exactly in the same way as c.best_val[c.iter] in line 196).

Note that I am almost sure that this does not affect the functioning of the BGP algorithm, nor of any other function in this package. It only makes it a bit more complicated to retrieve the iteration at which the lowest value has been found.

It should be any easy fix, just wanted to check with you that this is indeed a minor issue.

change swaprows

is hard to understand what's going on if i change the entire row. should only exchange selected elements of, i.e. the state of the chain according to BGP. when exchanging i with j, exchange

the parameter vector theta:

parameters(chain[i],iteration)[par_names] = parameters(chain[j],iteration)[par_names]

the summary statistic S, i.e. the moments:

moments(chain[i],iteration)[mom_names] = moments(chain[j],iteration)[mom_names]

leave a reference in infos[iteration, :exhanged_with]
but don't overwrite the entire row: chain_id, tempring etc stay fixed with the current chain.

Typo in definition of dataMomentW (Eval.jl)

Line 197 of Eval.jl is currently:

dataMomentW(ev::Eval) = dataMomentW(ev,collect(keys(ev.dataMomentW)))

When called, this function produces the error

ERROR: type Eval has no field dataMomentW Stacktrace:
[1] getproperty(::Any, ::Symbol) at ./Base.jl:20
[2] dataMomentW(::SMM.Eval) at /home/gregorymartin/.julia/packages/SMM/MDWS3/src/mopt/Eval.jl:197
[3] top-level scope at none:0

The field is dataMomentsW rather than dataMomentW. Corrected version:

dataMomentW(ev::Eval) = dataMomentW(ev,collect(keys(ev.dataMomentsW)))

ERROR: unsupported or misplaced expression =>

test fails with

ERROR: unsupported or misplaced expression =>

not sure what that is about. i think julia 0.4 has a new dict? we should require julia 0.3 only for now i think.

should leverage TuringLang

https://github.com/TuringLang/MCMCChains.jl

change exchange moves

the exchange moves should be independent of the distance between chains i and j. it should only depend on whether the proposed chain j is acceptable on chain i:

# this
if algo.dist_fun(evj.value - evi.value)  < algo["maxdists"][i]
            swap_ev_ij!(algo,i,j)
        end

# becomes
if evj.value < algo["maxdists"][i]
            swap_ev_ij!(algo,i,j)
end

need maybe to tighten up MProb to require only maximization or minimization

Example not producing correct results

@floswald I'm running the example Jupyter notebook, but it's taking almost 80 minutes and I'm getting the following as output:

SMM.summary(MA) = 3×5 DataFrame
│ Row │ id    │ acc_rate │ perc_exchanged │ exchanged_most_with │ best_val │
│     │ Int64 │ Float64  │ Float64        │ Int64               │ Float64  │
├─────┼───────┼──────────┼────────────────┼─────────────────────┼──────────┤
│ 1   │ 1     │ 0.001    │ 0.0            │ 0                   │ -1.0     │
│ 2   │ 2     │ 0.001    │ 0.0            │ 0                   │ -1.0     │
│ 3   │ 3     │ 0.001    │ 0.0            │ 0                   │ -1.0     │
Quasi posterior mean for p1 = 0.2
Quasi posterior median for p1 = 0.2
True value = -1.0
abs(True value - Quasi posterior median)  = 1.2

Quasi posterior mean for p2 = -0.2
Quasi posterior median for p2 = -0.2
True value = 1.0
abs(True value - Quasi posterior median)  = 1.2

Do you happen to know what might be going on? I'm using Julia 1.5.0. I would like to use your package to teach my PhD students (and myself!) how to do SMM. If you know of other "simple" examples (e.g. binary logit estimation via SMM or linear regression via SMM) I would love to use those as well!

The exact code I ran is below. The main thing I had to change was to use Random.seed() rather than srand() which no longer works.

# add three workers
addprocs(3)

@everywhere using SMM
using DataStructures
using DataFrames
using Plots
using Random

Random.seed!(1234);

# Define true values of parameters
#---------------------------------
trueValues = OrderedDict("mu1" => [-1.0], "mu2" => [1.0])
#------------------------------------------------
# Options
#-------------------------------------------------
# Boolean: do you want to save the plots to disk?
savePlots = false

#------------------------
# initialize the problem:
#------------------------

# Specify the initial values for the parameters, and their support:
pb = OrderedDict("p1" => [0.2,-3,3] , "p2" => [-0.2,-2,2])

# Specify moments to be matched + subjective weights:
moms = DataFrame(name=["mu1","mu2"],value=[trueValues["mu1"][], trueValues["mu2"][]], weight=ones(2))

# GMM objective function to be minized.
# It returns a weigthed distance between empirical and simulated moments
@everywhere function objfunc_normal(ev::Eval; verbose = false)

    start(ev)


    # when running in parallel, display worker's id:
    #-----------------------------------------------
    if verbose == true
        if nprocs() > 1
          println(myid())
        end
    end

    # extract parameters from ev:
    #----------------------------
    mu  = collect(values(ev.params))

    # compute simulated moments
    #--------------------------
    # Monte-Carlo:
    #-------------
    ns = 100000 #number of i.i.d draws from N([mu], sigma)
    #initialize a multivariate normal N([mu], sigma)
    #sigma is set to be the identity matrix
    sigma = [1.0; 1.0]
    # draw ns observations from N([mu], sigma):
    randMultiNormal = SMM.MvNormal(mu,SMM.PDiagMat(sigma))
    # calculate the mean of the simulated data
    simM            = mean(rand(randMultiNormal,ns),2)
    # store simulated moments in a dictionary
    simMoments = Dict(:mu1 => simM[1], :mu2 => simM[2])

    # Calculate the weighted distance between empirical moments
    # and simulated ones:
    #-----------------------------------------------------------
    v = Dict{Symbol,Float64}()
    for (k, mom) in dataMomentd(ev)
        # If weight for moment k exists:
        #-------------------------------
        if haskey(SMM.dataMomentWd(ev), k)
            # divide by weight associated to moment k:
            #----------------------------------------
            v[k] = ((simMoments[k] .- mom) ./ SMM.dataMomentW(ev,k)) .^2
        else
            v[k] = ((simMoments[k] .- mom) ) .^2
        end
    end

    # Set value of the objective function:
    #------------------------------------
    setValue(ev, mean(collect(values(v))))

    # also return the moments
    #-----------------------
    setMoment(ev, simMoments)

    # flag for success:
    #-------------------
    ev.status = 1

    # finish and return
    finish(ev)

    return ev
end



# Initialize an empty MProb() object:
#------------------------------------
mprob = MProb()

# Add structural parameters to MProb():
# specify starting values and support
#--------------------------------------
addSampledParam!(mprob,pb)

# Add moments to be matched to MProb():
#--------------------------------------
addMoment!(mprob,moms)

# Attach an objective function to MProb():
#----------------------------------------
addEvalFunc!(mprob, objfunc_normal)
# estimation options:
#--------------------
# number of iterations for each chain
niter = 1000
# number of chains
nchains = 3

opts = Dict("N"=>nchains,
        "maxiter"=>niter,
        "maxtemp"=> 5,
        "coverage"=>0.025,
        "sigma_update_steps"=>10,
        "sigma_adjust_by"=>0.01,
        "smpl_iters"=>1000,
        "parallel"=>true,
        "maxdists"=>[0.05 for i in 1:nchains],
        "mixprob"=>0.3,
        "acc_tuner"=>12.0,
        "animate"=>false)


#---------------------------------------
# Let's set-up and run the optimization
#---------------------------------------
# set-up BGP algorithm:
MA = MAlgoBGP(mprob,opts)

# run the estimation:
@time SMM.run!(MA);

# show a summary of the optimization:
@show SMM.summary(MA)
# # Plot histograms for the first chain, the one with which inference should be done.
# # Other chains are used to explore the space and avoid local minima
# #-------------------------------------------------------------------------------
# p1 = histogram(MA.chains[1])
# display(p1)
# 
# if savePlots
#     savefig(p1, joinpath(pwd(),"histogram.svg"))
# end
# 
# # Plot the realization of the first chain
# #----------------------------------------
# p2 = plot(MA.chains[1])
# if savePlots
#     savefig(p2, joinpath(pwd(),"lines.svg"))
# end
# display(p2)


# Realization of the first chain:
#-------------------------------
dat_chain1 = SMM.history(MA.chains[1])

# discard the first 10th of the observations ("burn-in" phase):
#--------------------------------------------------------------
dat_chain1[round(Int, niter/10):niter, :]

# keep only accepted draws:
#--------------------------
dat_chain1 = dat_chain1[dat_chain1[:accepted ].== true, : ]

# create a list with the parameters to be estimated
parameters = [Symbol(String("mu$(i)")) for i=1:2]
# list with the corresponding priors:
#------------------------------------
estimatedParameters = [Symbol(String("p$(i)")) for i=1:2]


# Quasi Posterior mean and quasi posterior median for each parameter:
#-------------------------------------------------------------------
for (estimatedParameter, param) in zip(estimatedParameters, parameters)

  println("Quasi posterior mean for $(String(estimatedParameter)) = $(mean(dat_chain1[estimatedParameter]))")
  println("Quasi posterior median for $(String(estimatedParameter)) = $(median(dat_chain1[estimatedParameter]))")
  println("True value = $(trueValues[String(param)][])")
  println("abs(True value - Quasi posterior median)  = $(abs(median(dat_chain1[estimatedParameter]) - trueValues[String(param)][])) \n")

end

Hope some econ examples can be provided

Hi, Florian:
This is a really cool stuff to make SMM estimation fesible to econ phD like me.
However, after reading the document, I still cannot easily map the code to the steps described in your paper. In your JMP and “Consumer Bankruptcy and Mortgage Default”, a parameter are estimated by matching one certain simulated moment to real moments from data, or even via multiple moments,which part of code correspond to this step? (I believe it's in the objective function). The arguments for addMoment are just two DataArrays, I don't quite understand what they represent.
Moreover, you also add some auxiliary models in these papers, can we now implement this in MOpt.jl?
Anyway, you guys did a great job.

Error tagging new release

The tag name "v0.6" is not of the appropriate SemVer form (vX.Y.Z).
cc: @floswald

Issue with parameters "labels"?

Hi,

With Edoardo Ciscato, we noticed something strange. It seems like the library "exchanges" the parameters to be estimated at some point (at least their names).

Here is an example, where I want to estimate the mean of 4-dimensional Gaussian process (I fix the variance-covariance matrix to be the identity matrix):

using MomentOpt
using GLM
using DataStructures
using DataFrames
using Plots
#plotlyjs()
pyplot()

#------------------------------------------------
# Options
#-------------------------------------------------
# Boolean: do you want to save the plots to disk?
savePlots = true

# initialize the problem:
# 4 means
#------------------------
# initial values:
#----------------
pb = OrderedDict("p1" => [0.2,-3,3] , "p2" => [-0.2,-2,2], "p3" => [0.1,0,10], "p4" => [-0.1,-10,0] )
# moments to be matched:
#-----------------------
moms = DataFrame(name=["mu1","mu2","mu3", "mu4"],value=[-1.0,1.0, 5.0, -4.0], weight=ones(4))



"""
    objfunc_normal(ev::Eval)

    GMM objective function to be minized.
    It returns a weigthed distance between empirical and simulated moments

    copy-paste of the function objfunc_norm(ev::Eval)
    I only made minor modifications to the original fuction

"""

function objfunc_normal(ev::Eval)

    start(ev)


    # extract parameters from ev:
    #----------------------------
    mu  = collect(values(ev.params))

    # compute simulated moments
    #--------------------------
    # Monte-Carlo:
    #-------------
    ns = 10000 #number of i.i.d draws from N([a,b], sigma)
    #initialize a multivariate normal N([a,b], sigma)
    #using a = mu[1], b=mu[2]
    sigma = [1.0 ;1.0; 1.0; 1.0]
    randMultiNormal = MomentOpt.MvNormal(mu,MomentOpt.PDiagMat(sigma))
    simM            = mean(rand(randMultiNormal,ns),2) #mean of simulated data
    simMoments = Dict(:mu1 => simM[1], :mu2 => simM[2], :mu3 => simM[3], :mu4 => simM[4])#store simulated moments in a dictionary


    # Calculate the weighted distance between empirical moments
    # and simulated ones:
    #-----------------------------------------------------------
    v = Dict{Symbol,Float64}()
    for (k, mom) in dataMomentd(ev)
        # If weight for moment k exists:
        #-------------------------------
        if haskey(MomentOpt.dataMomentWd(ev), k)
            # divide by weight associated to moment k:
            #----------------------------------------
            v[k] = ((simMoments[k] .- mom) ./ MomentOpt.dataMomentW(ev,k)) .^2
        else
            v[k] = ((simMoments[k] .- mom) ) .^2
        end
    end

    # Set value of the objective function:
    #------------------------------------
    setValue(ev, mean(collect(values(v))))

    # also return the moments
    #-----------------------
    setMoment(ev, simMoments)

    # flag for success:
    #-------------------
    ev.status = 1

    # finish and return
    finish(ev)

    return ev
end



# Initialize an empty MProb() object:
#------------------------------------
mprob = MProb()

# Add structural parameters to MProb():
# specify starting values and support
#--------------------------------------
addSampledParam!(mprob,pb)

# Add moments to be matched to MProb():
#--------------------------------------
addMoment!(mprob,moms)

# Attach an objective function to MProb():
#----------------------------------------
addEvalFunc!(mprob, objfunc_normal)


# estimation options:
#--------------------
# number of iterations for each chain
niter = 1000
# number of chains
# nchains = nprocs()
nchains = 4

opts = Dict("N"=>nchains,
        "maxiter"=>niter,
        "maxtemp"=> 5,
        "coverage"=>0.025,
        "sigma_update_steps"=>10,
        "sigma_adjust_by"=>0.01,
        "smpl_iters"=>1000,
        "parallel"=>true,
        "maxdists"=>[0.05 for i in 1:nchains],
        "mixprob"=>0.3,
        "acc_tuner"=>12.0,
        "animate"=>false)


# plot slices of objective function
#---------------------------------
s = doSlices(mprob,10)

# plot objective function over param values:
#-------------------------------------------
p1 = MomentOpt.plot(s,:value)

if savePlots == true
    Plots.savefig(p1, joinpath(pwd(),"slices_Normal1.svg"))
end


#---------------------------------------
# Let's set-up and run the optimization
#---------------------------------------
# set-up BGP algorithm:
MA = MAlgoBGP(mprob,opts)

# run the estimation:
@time MomentOpt.runMOpt!(MA)

# show a summary of the optimization:
@show MomentOpt.summary(MA)

# Plot histograms for the first chain:
#------------------------------------
p4 = histogram(MA.chains[1])

if savePlots == true
    savefig(p4, joinpath(pwd(),"histogram_chain1.svg"))
end

# Plot the "history" of one chain:
#--------------------------------
p10 = plot(MA.chains[1])
if savePlots == true
    savefig(p10, joinpath(pwd(),"history_chain_1.svg"))
end

The algorithm recovers the true parameters ([-1.0,1.0, 5.0, -4.0]) but not in the right order. It finds p1=5.0, p2 = -1, p3 = -4, p4 = 1.0, while it should be p1=-1.0, p2=1.0, p3=5.0, p4 = -4.0.

Is it an issue with the plotting functions, or something more serious?

Edit: on the other hand, slicing works fine

after this has been merged: JuliaStats/GLM.jl#242
can take on #24 only after update to 0.7

proposal() is difficult to handle with high-dim p

suppose you have 10 parameters in p
suppose they all have widely different supports
currently the proposal often fails to find a new candidate. it's a MVNormal that needs to be censored into the bounds of each param dimension, centered on the currentl value of p
particularly if the current draw is close to the boundary, that makes it very unlikely to find all draws inside the boundaries.

plotting

need to update plotSlices function

floswald / smm.jl Goto Github PK

smm.jl's Introduction

SMM.jl: Simulated Method of Moments for Julia

Installation

Documentation

Examples

Contributing

History

Thanks to all Contributors!

smm.jl's People

Contributors

Stargazers

Watchers

Forkers

smm.jl's Issues

Recommend Projects

Recommend Topics

Recommend Org