# Julia Manual GMM

Hey
I am having a bit of trouble understanding how to manually do GMM and in particular how to do GMM in Julia.
Here is the Stata code + data:

use http://www.stata-press.com/data/r13/docvisits
gmm (docvis - exp({xb:private chronic female income}+{b0})), instruments(private chronic female age black hispanic) twostep

As I understand it, I need to first find the residuals and then find estimates for theta but am not sure on how to do either. I am supposed to first assume i.i.d errors and then assume heteroskedasticity but no autocorrelation in the errors.

Is it something of the form:
inv(n)Z’(y-Xtheta))'W(inv(n)Z’(y-Xtheta))

Where W is the weighting matrix for the heteroskedastic case and W=I for i.i.d
Any help would be appreciated,
Thanks

EDIT: Here is what I have so far, not sure how correct it is:

using Optim
using CSV

docvis=Array(data[:docvis])
y=docvis
#x and z variables
private=Array(data[:private])
chronic=Array(data[:chronic])
female=Array(data[:female])
income=Array(data[:income])
age=Array(data[:age])
black=Array(data[:black])
hispanic=Array(data[:hispanic])
β0=ones(4412)
#collection of instruments
Z=[private chronic female age black hispanic]

#collection of regressors

X=[β0 private chronic female income]

“”"
L>K so X=Z inv(Z’Z)Z’X
“”"
X̂=(Z*inv(Z’*Z)*Z’*X)
W=eye(6)
θ=(inv(X̂’*X̂)*X̂’y)
#The solution for Optimal β
“”"
β=inv(X’ZW
Z’*X)X’ZWZ’y
eX=exp(X)
βe=inv(eX’ZW
Z’*eX)eX’ZWZ’y
“”"
u=y-X̂
θ
σ= u.^2

f(σ)=y-X̂*θ
optimize(f, [0, 0, 0, 0, 0])

You can include code on Discourse using code blocks,

using Optim
using CSV

docvis=Array(data[:docvis])
y=docvis
#x and z variables
private=Array(data[:private])
chronic=Array(data[:chronic])
female=Array(data[:female])
income=Array(data[:income])
age=Array(data[:age])
black=Array(data[:black])
hispanic=Array(data[:hispanic])
β0=ones(4412)
#collection of instruments
Z=[private chronic female age black hispanic]

#collection of regressors

X=[β0 private chronic female income]

"""
L>K so X=Z inv(Z’Z)Z’X
"""
X̂=(Z*inv(Z’*Z)*Z’*X)
W=eye(6)
θ=(inv(X̂’*X̂)*X̂’y)
#The solution for Optimal β
"""
β=inv(X’ZWZ’*X)X’ZWZ’y
eX=exp(X)
βe=inv(eX’ZWZ’*eX)eX’ZWZ’y
"""
u=y-X̂θ
σ= u.^2

f(σ)=y-X̂*θ
optimize(f, [0, 0, 0, 0, 0])


You define a function f(σ)=y-X̂*θ that does not depend on the input variable σ. So your function is constant. I think you may benefit from a better understanding of how to define methods and functions in Julia, see the resources at https://julialang.org/learning/

Optim works well to minimise least squares problems, but you can also consider using LsqFit.jl

I have a function which will do this, pretty easily, at https://github.com/mcreel/Econometrics.jl/blob/master/src/GMM/gmmresults.jl

It gives results as follows:

julia> gmmresults();
************************************************************
GMM example, two step
GMM Estimation Results    BFGS convergence: Normal
Observations: 100
Hansen-Sargan statistic: 1.84211
Hansen-Sargan p-value: 0.1747

estimate     st. err      t-stat     p-value
1    -0.13380     0.10507    -1.27341     0.20588
2     1.01924     0.15808     6.44768     0.00000

************************************************************
************************************************************
GMM example, CUE
GMM Estimation Results    BFGS convergence: Normal
Observations: 100
Hansen-Sargan statistic: 1.56835
Hansen-Sargan p-value: 0.21045

estimate     st. err      t-stat     p-value
1    -0.10846     0.10572    -1.02592     0.30745
2     0.98518     0.14669     6.71591     0.00000

************************************************************

julia>



For your problem, you would just need to define the appropriate moments function. An example which uses IV to estimate a linear model similar to yours is here: https://github.com/mcreel/Econometrics/blob/master/Examples/GMM/MeasurementErrorIV.jl

By the way, the gmmresults() function is general purpose, for possibly nonlinear in the parameters moment conditions. As such, it’s using Optim.jl, which is not necessary in the case of linear in parameters moment conditions. But, it works.

3 Likes

From the Stata Code I provided earlier, am I understanding correctly that the moment conditions would be something like:

[y-thetaexp(X) y-thetaexp(X) … y-thetaexp(X)]

Where there are 5 moment conditions?

Just found out that we are required to use

Optim
Optimize

Any thoughts on this?

β0=ones(n)
#collection of instruments
Z=[private chronic female age black hispanic]

X=[β0 private chronic female income]
W1=eye(6)

Q(θ)=(inv(n)(Z)(y-Xθ’))'W1(inv(n)(Z)(y-Xθ’))

optimize(Q)

Any julia GMM package now?

I’d say GMM is something that doesn’t really need a package. The estimator is defined by two things, the moment conditions, and the weight matrix (or just the moment conditions, if the continuously updating version is used). The following code will compute the ordinary two-step GMM estimator and everything needed to compute the estimated covariance matrix, given a function defining the moment conditions and the weight matrix (fminunc is just a convenience frontend to Optim). There’s not much there to need a dedicated package. I’d say it would be more appropriate to place GMM into a more general econometrics package.

# ordinary GMM with provided weight matrix
# weight should be gXg positive definite
function gmm(moments, theta, weight)
# average moments
m = theta -> vec(mean(moments(theta),dims=1)) # 1Xg
# moment contributions
momentcontrib = theta -> moments(theta) # nXg
# GMM criterion
obj = theta -> ((m(theta))'weight*m(theta))
# do minimization
thetahat, objvalue, converged = fminunc(obj, theta)
# derivative of average moments
D = (Calculus.jacobian(m, vec(thetahat), :central))'
# moment contributions at estimate
ms = momentcontrib(thetahat)
return thetahat, objvalue, D, ms, converged
end


4 Likes

@Yifan_Liu
This is the best one I found so far:

Also see:

By @floswald

I believe the example by @mcreel is very good. In case you want more of something close that, you may consider my lecture notes/notebook at GMM forth and back

3 Likes

These are all the Julia codes I’m aware of which can do GMM (none is a registered pkg):
@Paul_Schrimpf: https://github.com/schrimpf/GMMInference.jl

magnus dahlquist: https://github.com/MagnusDahlquist/PhD403

Chris Hedtrich: https://github.com/chrished/NonlinearGMM.jl

1 Like

thanks for your reply. I wrote my GMM functions in a style similar to the ones Paul_Soderlind mentioned in this thread.

The reason I am asking this is that I wan to see if there is another package that can export the results of a GMM function to LaTex tables.

1 Like

The function gmmresults.jl, mentioned a while ago upthread, uses the gmm code I posted recently, above, to compute the estimator, then it prints out nicely formatted results. To do this, it uses the PrettyTables.jl package. This package allows export to LaTeX. gmmresults.jl doesn’t use this feature, but it should be quite easy to modify the code to add this option.

By the way, thanks to all for the links, there are lots of good resources out there!

2 Likes

Just curious, does this PrettyTables package have the potential of translating the results of any econometrics functions into LaTex tables?

As far as I know, it will just nicely format tables, and has the option to export LaTeX. But the user must supply the basic table input.

I’m working on my own implementation of GMM for a project and thought I’d contribute to the many excellent implementations in this thread with a minimal example that makes good use of unicode characters to correspond to familiar mathematical notation. This example assumes you have already programmed a function g(z,θ) that returns an L \times T matrix of L moment contributions computed at T observations (i.e. the t-th column of the output of g is equal to g(z_t,\theta).

using Optim, ForwardDiff
function gmm(g, z, θ₀, W)
# Parameter estimates
f(θ) = (ḡ = mean(g(z,θ), dims=2); dot(ḡ, W, ḡ)) # GMM criterion
result = Optim.optimize(f, θ₀, LBFGS(); autodiff=:forward) # run optimization
θ̂ = Optim.minimizer(result) # access result

# Standard errors: (G'WG)⁻¹G'WΩ̂WG(G'WG)⁻¹ →ᵖ Var(√T̅(θ̂-θ))
T = size(z,1)
G = ForwardDiff.jacobian(θ -> mean(g(z,θ), dims=2), θ̂) # G →ᵖ E[∂g/∂θ′]
D = inv(G'*W*G)
Ω̂ = (ĝ = g(z,θ̂); ĝ*ĝ'/T) # Ω̂ = T⁻¹Σₜ[g(zₜ,θ̂)g(zₜ,θ̂)′] →ᵖ Ω
V̂ = D*G'*W*Ω̂*W*G*D # (G'WG)⁻¹G'WΩ̂WG(G'WG)⁻¹
σ̂ = sqrt.(diag(V̂/T))

return θ̂, σ̂
end


You can test it out by comparing it to OLS using this moment function:

function g(z,θ)
y = @view z[:,1]
X = @view z[:,2:end]
ε̂ = y - X*θ
transpose(X.*ε̂)
end


Understanding the mechanics of actually getting GMM to work was a bit of a hurdle for me coming from Matlab, so I hope this is helpful to others making the switch.

5 Likes

Hello. I am trying the same but my codes are not working. How to introduce the weight W=eye(6) to my model?

The below is a code for my case
using Pkg
import Pkg;
using DataFrames;
using CSV;
using Lathe;
using GLM;
using Statistics;
using StatsPlots;
using MLBase;
df=DataFrame(CSV.File(“C:Users/GABBYMOSH/OneDrive/Desktop/Work/hsng.csv”)
first(df,5);
names(df)
#Multiple linear Model
fm=@formula(rent~hsngval+pcturban)
linearRegressor=lm(fm,df)
##GMM
hsngval=Array(df[:hsngval]);
pcturban=Array(df[:pcturban]);
faminc=Array(df[:faminc]);
reg2=Array(df[:reg2]);
reg4=Array(df[:reg4]);
#Instrument variables
Z=[pcturban faminc reg2 reg4];
B0=ones(50);
X=[B0 hsngval pcturban];
rent=Array(df[:rent]);
y=rent;

I can’t answer everything, but this might help a bit: Assuming eye(6) is supposed to return a 6x6 floating point identity matrix, then this will work
using LinearAlgebra: I