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:
Where W is the weighting matrix for the heteroskedastic case and W=I for i.i.d
Any help would be appreciated,
EDIT: Here is what I have so far, not sure how correct it is:
data=CSV.read(“data Ass_2.csv”, datarow=2)
y=docvis #x and z variables
β0=ones(4412) #collection of instruments
Z=[private chronic female age black hispanic]
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
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.
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
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.
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 θ̂, σ̂
You can test it out by comparing it to OLS using this moment function:
y = @view z[:,1]
X = @view z[:,2:end]
ε̂ = y - X*θ
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.
The below is a code for my case
names(df) #Multiple linear Model
reg4=Array(df[:reg4]); #Instrument variables
Z=[pcturban faminc reg2 reg4];
X=[B0 hsngval pcturban];