How to correct standard errors standard errors in a multinomial logit using IV

Question

I am trying to estimate a multinomial logit model using an instrumental variable. I didn't find any preexisting package, so I tried to estimate using a two-stage approach. First estimating the first stage as an OLS with the IV tsls1<-lm(d~x+z) Then I used d.hat<-fitted.values(tsls1) With that, I used multinom function from the nnet. tsls2<-multinom(y~x+d.hat) The problem is that the standard errors are wrong. I was wondering how I could correct them. Or if there is an easier way.